text:
function model = setup_fmincon_params(model)

monomtable = model.monomtable;
nonlinearindicies = model.nonlinearindicies;
linearindicies = model.linearindicies;

model.evalinobjective  = 0;

if ~isempty(model.evalMap)
    [i,j,k] = find(model.monomtable(:,model.evalVariables));
    evalInvolved = union(model.evalVariables,i);
    model.evalinobjective = any(model.c(evalInvolved)) | any(model.Q(evalInvolved,1));
    if size(model.Anonlineq,1)>0
        if nnz(model.Anonlineq(:,evalInvolved))>0
            model.evalinconstraint = 1;
        end
    end
    if size(model.Anonlinineq,1)>0
        if nnz(model.Anonlinineq(:,evalInvolved))>0
            model.evalinconstraint = 1;
        end
    end
    if any(model.K.q)
        if nnz(model.F_struc(1+model.K.f + model.K.l:end,1+evalInvolved)) > 0
            model.evalinconstraint = 1;
        end
    end
    % Speed up callbacks by removing last marker argument (only used in
    % expandmodel etc) and figuring out what R^m -> R^n operators computes
    for i = 1:length(model.evalMap)
        model.evalMap{i}.prearg = {model.evalMap{i}.fcn,model.evalMap{i}.arg{1:end-1}};
        model.evalMap{i}.prearg{1+model.evalMap{i}.argumentIndex} = [];
        if ~isfield(model.evalMap{i},'computes')
            model.evalMap{i}.computes = model.evalVariables(i);
        end
    end
end

% Figure out if YALMIP easily can compute the gradient of the objective
% This will done completely general later
model.SimpleLinearObjective = 0;    % Linear
model.SimpleQuadraticObjective = 0; % Quadratic
model.SimpleNonlinearObjective = 1; % Polynomial
if isempty(model.evalMap)
    if nnz(model.c(nonlinearindicies)) == 0
        if (nnz(model.Q)==0)
            model.SimpleLinearObjective = 1;
        else
            if nnz(model.Q(nonlinearindicies,nonlinearindicies))==0
                model.SimpleQuadraticObjective = 1;
            end
        end
    end
elseif ~model.evalinobjective
    if nnz(model.c(nonlinearindicies)) == 0
        if (nnz(model.Q)==0)
            model.SimpleLinearObjective = 1;
        else
            if nnz(model.Q(nonlinearindicies,nonlinearindicies))==0
                model.SimpleQuadraticObjective = 1;
            end
        end
    elseif nnz(model.Q)==0
        r = find(model.c);
        if all(model.variabletype(r)<=2)
            if isempty(intersect(r,model.evalVariables))
                % Simple quadratic function
                for i = r(:)'
                    if model.variabletype(i)==2
                        j = find(model.monomtable(i,:));
                        model.Q(j,j) = model.c(i);
                        model.c(i) = 0;
                    elseif model.variabletype(i)==1
                        j = find(model.monomtable(i,:));
                        model.Q(j(1),j(2)) = model.c(i)/2;
                        model.Q(j(2),j(1)) = model.c(i)/2;
                        model.c(i) = 0;
                    end
                end
                model.SimpleQuadraticObjective = 1;
            end
        end
    end
else
    model.SimpleNonlinearObjective = 0;
end

model.linearconstraints = isempty(model.Anonlinineq) & isempty(model.Anonlineq) & nnz(model.K.q)==0;
model.nonlinearinequalities = ~isempty(model.Anonlinineq);
model.nonlinearequalities = ~isempty(model.Anonlineq);
if any(model.K.q)
    if nnz(model.F_struc(1+model.K.f + model.K.l:end,1+model.nonlinearindicies)) > 0
        model.nonlinearcones = 1;
    else
        model.nonlinearcones = 0;
    end
else
    model.nonlinearcones = 0;
end

% Structure for fast evaluation of monomial terms in differentiation
 if isempty(model.evalMap) & (model.nonlinearinequalities | model.nonlinearequalities | model.nonlinearcones) & ~isfield(model,'fastdiff')
    allA = [model.Anonlineq;model.Anonlinineq];
    dgAll = [];
    n = length(model.c);
    linearindicies = model.linearindicies;
    mtNonlinear = model.monomtable(model.nonlinearindicies,:);

    allDerivemt = [];
    news = [];
    c = [];
    mtNonlinear = mtNonlinear(:,linearindicies);
    mtNonlinear = mtNonlinear';

    [jj,ii,val] = find(mtNonlinear');
    for k = 1:length(ii)
        i = ii(k);
        j = jj(k);
        s=mtNonlinear(:,j);
        c = [c;s((i))];
        s((i)) = s((i))-1;
        allDerivemt = [allDerivemt s(:)];
        news = [news;j i];
    end

    allDerivemt = allDerivemt';

    model.fastdiff.news = news;
    model.fastdiff.allDerivemt = allDerivemt;
    model.fastdiff.c = c;
    model.fastdiff.univariateDifferentiates = 0;

    a1 =  model.fastdiff.news(1:length(model.fastdiff.c),2);
    a2 =  model.nonlinearindicies(model.fastdiff.news(1:length(model.fastdiff.c),1))';
    zzz = [ones(length(a1),1)];
    nn = max(max(length(model.linearindicies)),max(news(:,2)));
    mm = max(max(linearindicies),max(model.nonlinearindicies(news(:,1))));
    a1f = [a1(:);(1:length(model.linearindicies))'];
    a2f = [a2(:);model.linearindicies(:)];
    zzzf = [zzz;ones(length(linearindicies),1)];

    %model.fastdiff.newdxx = sparse(a1f,a2f,zzzf,nn,mm);
    %model.fastdiff.linear_in_newdxx = sub2ind([nn mm],a1(:),a2(:));
    model.fastdiff.newdxx = sparse(a1f,a2f,zzzf,nn,mm)';
    model.fastdiff.linear_in_newdxx = sub2ind([mm nn],a2(:),a1(:));
    model.fastdiff.newdxx(model.fastdiff.linear_in_newdxx)=1;

    if all(sum(allDerivemt | allDerivemt,2)==1)
        model.fastdiff.univariateDifferentiates = 1;
        [i,j,k] = find(allDerivemt');
        if any(diff(j)<0)
            error;
        end
        model.fastdiff.univariateDiffMonom = i(:);
        model.fastdiff.univariateDiffPower = k(:);
    end
 else
      allA = [model.Anonlineq;model.Anonlinineq];
      if any(model.K.q)
          allA = [allA;model.F_struc(1+model.K.f + model.K.f:end,2:end)];
      end
      requested = any(allA',2);
      [i,j,k] = find((model.deppattern(find(requested),:)));
      requested(j) = 1;
      model.fastdiff.requested = requested;
 end

author: zarathustr

repo: LibQPEP


text:
function dx = dynamics(x,u,p)
% dx = dynamics(x,u,p)
%
% Computes the dynamics for the simple pendulum
%

q = x(1,:);
dq = x(2,:);

k = p.k;    c = p.c;
ddq = -c*dq - k*sin(q) + u;
dx = [dq;ddq];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function model = assignGPRA(model,gpraFile)
% Assigns each reaction a Gene-protein-reaction association
%
% USAGE:
%
%    model = assignGPRA(model, gpraFile)
%
% INPUTS:
%    model:       COBRA model structure
%    gpraFile:    SimPheny GPRA file
%
% OUTPUT:
%    model:       COBRA model with gene-protein-reaction assoction
%
% .. Author: -  Markus Herrgard 10/4/06

[rxnList,genes,rxnGeneMat,subSystems] = parseSimPhenyGPRA(gpraFile);

[tmp,gpraRxnInd,modelRxnInd] = intersect(rxnList,model.rxns);

model.genes = genes;

model.rxnGeneMat = sparse(length(model.rxns),length(genes));
model.rxnGeneMat(modelRxnInd,:) = rxnGeneMat(gpraRxnInd,:);
model.rxnGeneMat=sparse(model.rxnGeneMat);

author: opencobra

repo: cobratoolbox


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function f = RegCoeff(S, g, df, B, Nb, Nr)
% calculates regression coefficients to be used with longstaff-schwartz

v = g(:,end);   % start for backward induction

f = zeros(Nb, Nr-1);

% backward induction and regression from t_{Nr-1} up to t_1
for i = Nr-1:-1:1
        index = find(g(:,i) > 0); % all ITM paths
        s = S(index,i+1);         % values of S at given time point
        v = v * df(i+1);          % option value at t_i

        Acell = B(s);             % evaluate basis function in cell array B
        A = cell2mat(Acell{:,:}); % convert to matrix

        f(:,i) = (A'*A)\(A'*v(index)); % determine coefficients
        c = A*f(:,i);                   % continuation value
        exercise = g(index,i) >= c;    % early exercise
        v(index(exercise)) = g(index(exercise),i);
end

end

author: Sable

repo: mcbench-benchmarks


text:
function ostr = mlstripcommentsstr(istr)
%MLSTRIPCOMMENTSSTR Strip comments from a string with MATLAB code.
%
%   OSTR = MLSTRIPCOMMENTSSTR(ISTR) takes the input string ISTR, assuming it
%   contains MATLAB code, strips all MATLAB comments, and returns the result as
%   the output string OSTR.  The input string must be a row vector or a cell
%   array of row vectors.  Any input string may be a multi-line string, that is
%   a string containing newline characters, for example a string containing the
%   entire contents of an m-file.
%
%   Example:  Read all MATLAB code from an m-file opened for reading on file
%   identifier "fid" and strip all comments:
%
%       str = fread(fid);               % data as numerical column vector
%       str = char(str)';               % data as character row vector
%       str = mlstripcommentsstr(str);  % strip all MATLAB comments
%
%   See also MLSTRIPCOMMENTSFILE, MLSTRIPCOMMENTSFID.

%   Author:      Peter J. Acklam
%   Time-stamp:  2004-03-18 08:23:19 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % Regex for removing a comment on a line possibly containing MATLAB code.
   %
   % Use a persistent variable so that the variable is set up only once.
   %
   % Find the part of the line which is not a comment and then replace the
   % whole line with the part that was not a comment.
   %
   % A comment is started by a percent sign, but not all percent signs start a
   % comment.  A percent sign might be embedded in a string, in which case the
   % percent sign obviously does not start a comment.  So we must look for
   % percent signs which are not embedded in strings.  To do this we must find
   % all strings.  Strings are delimited by single quotes, but not all single
   % quotes start strings.  A single quote might also be a tranpose operator.
   % So we must look for single quotes which are not tranpose operators.
   %
   % The MATLAB grammar (syntax and semantics) is not available to the public,
   % but from my knowledge, a single quote is a transpose operator if and only
   % if one of the following criteria is satisfied
   %
   %   1) it follows a right closing delimiter
   %   2) it follows right after a variable name or function name
   %   3) it follows right after a period (as part of the ".'" operator)
   %   4) it follows right after another tranpose operator
   %
   % This can be immediately after a right closing bracket ("]"), right closing
   % parenthesis (")"), right closing brace ("}"), an English letter ("a" to
   % "z" and "A" to "Z", a digit ("0" to "9"), an underline ("_"), a period
   % ("."), or another single quote ("'").

   persistent mainregex

   if isempty(mainregex)

      mainregex = [ ...
         ' (                   ' ... % Grouping parenthesis (content goes to $1).
         '   ( ^ | \n )        ' ... % Beginning of string or beginning of line.
         '   (                 ' ... % Non-capturing grouping parenthesis.
         '                     ' ...
         '' ... % Match anything that is neither a comment nor a string...
         '       (             ' ... % Non-capturing grouping parenthesis.
         '           [\]\)}\w.]' ... % Either a character followed by
         '           ''+       ' ... %    one or more transpose operators
         '         |           ' ... % or else
         '           [^''%]    ' ... %   any character except single quote (which
         '                     ' ... %   starts a string) or a percent sign (which
         '                     ' ... %   starts a comment).
         '       )+            ' ... % Match one or more times.
         '                     ' ...
         '' ...  % ...or...
         '     |               ' ...
         '                     ' ...
         '' ...  % ...match a string.
         '       ''            ' ... % Opening single quote that starts the string.
         '         [^''\n]*    ' ... % Zero or more chars that are neither single
         '                     ' ... %   quotes (special) nor newlines (illegal).
         '         (           ' ... % Non-capturing grouping parenthesis.
         '           ''''      ' ... % An embedded (literal) single quote character.
         '           [^''\n]*  ' ... % Again, zero or more chars that are neither
         '                     ' ... %   single quotes nor newlines.
         '         )*          ' ... % Match zero or more times.
         '       ''            ' ... % Closing single quote that ends the string.
         '                     ' ...
         '   )*                ' ... % Match zero or more times.
         ' )                   ' ...
         ' [^\n]*              ' ... % What remains must be a comment.
                  ];

      % Remove all the blanks from the regex.
      mainregex = mainregex(~isspace(mainregex));

   end

   % Initialize output.
   ostr = istr;

   % If input is a character array, put it into a cell array.  We'll later make
   % sure output is a character array if input is a character array.
   if ischar(ostr)
      ostr = {ostr};
   end

   % Iterate over each element in the cell array.
   for i = 1 : numel(ostr)

      % Get the ith input string.
      str = ostr{i};

      % All character arrays must be row vectors.
      if (ndims(str) > 2) | (any(size(str) > 0) & (size(str, 1) ~= 1))
         error('All character arrays must be row vectors.');
      end

      % Note:  With MATLAB 6.5.1 it seems that the '^' and '$' anchors in the
      % regexes '(^|\n)' and '($|\n)' are not working (bug?), so use a
      % workaround:  prepend and append a newline and remove them later.

      use_workaround = 1;

      if use_workaround
         lf   = sprintf('\n');     % LF = Line feed character
         str = [lf, str, lf];      % prepend and append an LF
      end

      % Remove comment lines.
      str = regexprep(str, '(^|\n)[^\S\n]*%[^\n]*', '$1', 'tokenize');

      % Remove a comment on a line possibly containing MATLAB code.
      str = regexprep(str, mainregex, '$1', 'tokenize');

      % Remove trailing whitespace from each line.
      str = regexprep(str, '[^\S\n]+($|\n)', '$1', 'tokenize');

      % Compress multiple blank lines.
      %str = regexprep(str, '\n\n\n+', sprintf('\n\n'));

      if use_workaround
         str = str(2:end-1);       % remove leading and trailing LF
      end

      % Insert ith string into output cell array.
      ostr{i} = str;

   end

   % If the input was a character array, make sure output is so too.
   if ischar(istr)
      ostr = ostr{1};
   end


text:
%EarthView Image from Google maps
%
% A concrete subclass of ImageSource that acquires images from Google
% maps.
%
% Methods::
% grab    Grab a frame from Google maps
% size    Size of image
% close   Close the image source
% char    Convert the object parameters to human readable string
%
% Examples::
% Create an EarthView camera
%          ev = EarthView();
% Zoom into QUT campus in Brisbane
%          ev.grab(-27.475722,153.0285, 17);
% Show aerial view of Brisbane in satellite and map view
%          ev.grab('brisbane', 14)
%          ev.grab('brisbane', 14, 'map')
%
% Notes::
% - Google limit the number of map queries limit to 1000 unique (different)
%   image requests per viewer per day.  A 403 error is returned if the daily
%   quota is exceeded.
% - Maximum size is 640x640 for free access, business users can get more.
% - There are lots of conditions on what you can do with the images,
%   particularly with respect to publication.  See the Google web site for
%   details.
%
% Author::
%  Peter Corke, with some lines of code from from get_google_map by Val
%  Schmidt.
%
% See also ImageSource.

%TODO:
% map features/element control
% street view
% roadmap, satellite, terrain, hybrid

classdef EarthView < ImageSource

% e = EarthView() is an object that returns images of the Earth's surface
% obtained from Google Maps.
%
% e = EarthView(googlekey) as above but the google key string is passed in.
%
% e.grab(lat, lon, zoom, options)

    properties
        key
        type
        scale
    end

    methods
        function ev = EarthView(varargin)
        %EarthView.EarthView Create EarthView object
        %
        % EV = EarthView(OPTIONS)
        %
        % Options::
        % 'satellite'    Retrieve satellite image
        % 'map'          Retrieve map image
        % 'hybrid'       Retrieve satellite image with map overlay
        % 'scale'        Google map scale (default 18)
        % 'width',W      Set image width to W (default 640)
        % 'height',H     Set image height to H (default 640)
        % 'key',S        The Google maps key string
        %
        % see also options for ImageSource.
        %
        % Notes::
        % - A key is required before you can use the Google Static Maps API.
        %   The key is a long string that can be passed to the constructor or
        %   saved as an environment variable GOOGLE_KEY.  You need a Google
        %   account before you can register for a key.
        %
        % Notes::
        % - Scale is 1 for the whole world, ~20 is about as high a resolution
        %   as you can get.
        %
        % See also ImageSource, EarthView.grab.

        %TODO
        %  clone this for new StreetView API
        %  method to return lat/long or NE matrices corresp to pixels

            ev = ev@ImageSource(varargin{:});

            opt.type = {'satellite', 'map', 'hybrid'};
            opt.scale = 18;
            opt.key = [];

            [opt,args] = tb_optparse(opt, varargin);

            ev.type = opt.type;
            ev.scale = 1;

            % set default size params if not set
            if isempty(ev.width)
                ev.width = 640;
                ev.height = 640;
            end

            if isempty(opt.key)
                opt.key = getenv('GOOGLE_KEY');
            end
            ev.key = opt.key;

        end

        function [im,E,N] = grab(ev, varargin)
            % EarthView.grab Grab an aerial image
            %
            % IM = EV.grab(LAT, LONG, OPTIONS) is an image of the Earth
            % centred at the geographic coordinate (lat, long).
            %
            % IM = EarthView.grab(LAT, LONG, ZOOM, OPTIONS) as above with the
            % specified zoom.  ZOOM is an integer between 1 (zoom right out)
            % to a maximum of 18-20 depending on where in the world you
            % are looking.
            %
            % [IM,E,N] = EarthView.grab(LAT, LONG, OPTIONS) as above but also
            % returns the estimated easting E and northing N.  E and N are both
            % matrices, the same size as IM, whose corresponding elements are
            % the easting and northing are the coordinates of the pixel.
            %
            % [IM,E,N] = EarthView.grab(NAME, OPTIONS) as above but uses a
            % geocoding web site to resolve the name to a location.
            %
            % Options::
            % 'satellite'      Retrieve satellite image
            % 'map'            Retrieve map image
            % 'roadmap'        Retrieve map image (synonym for 'map')
            % 'terrain'        Retrieve terrain map
            % 'hybrid'         Retrieve satellite image with map overlay
            % 'roads'          Retrieve a binary image that shows only roads, no labels.
            %                  Roads are white, everything else is black.
            % 'noplacenames'   Don't show placenames.
            % 'noroadnames'    Don't show roadnames.
            %
            %
            % Examples::
            % Zoom into QUT campus in Brisbane
            %          ev.grab(-27.475722,153.0285, 17);
            % Show aerial view of Brisbane in satellite and map view
            %          ev.grab('brisbane', 14)
            %          ev.grab('brisbane', 14, 'map')
            %
            % Notes::
            % - If northing/easting outputs are requested the function
            %   deg2utm is required (from MATLAB Central)
            % - The easting/northing is somewhat approximate, see
            %   get_google_map on MATLAB Central.
            % - If no output argument is given the image is displayed using idisp.

            opt.type = {'satellite', 'map', 'hybrid', 'terrain', 'roadmap', 'roads'};
            opt.scale = ev.scale;
            opt.roadnames = true;
            opt.placenames = true;
            opt.landscape = false;
            opt.onlyroads = false;

            [opt,args] = tb_optparse(opt, varargin);

            switch opt.type
                case 'map'
                    opt.type = 'roadmap'
                case 'roads'
                    opt.type = 'roadmap';
                    opt.onlyroads = true;
            end

            % build the URL
            if ischar(args{1})
                % given a string name, do a geocode lookup
                place = args{1};
                zoom = args{2};

                % build the URL, and load the XML document
                url = sprintf('http://maps.googleapis.com/maps/api/geocode/xml?address=%s&sensor=false', place);
                doc = xmlread(url);

                % walk the XML document
                locations = doc.getElementsByTagName('location')

                if locations.getLength > 1
                    fprintf('%d places called %s found\n', locations.getLength, place);
                end

                location = locations.item(0);   % take the first return

                assert(~isempty(location), 'location not found');

                node_lat = location.getElementsByTagName('lat');
                el = node_lat.item(0);
                lat = str2num( el.getFirstChild.getData );

                node_lon = location.getElementsByTagName('lng');
                el = node_lon.item(0);
                lon = str2num( el.getFirstChild.getData );
            else
                lat = args{1};
                lon = args{2};
                if length(args) == 3
                    zoom = args{3};
                else
                    zoom = 18;
                end

            end
            % now read the map
            baseurl = sprintf('https://maps.googleapis.com/maps/api/staticmap?center=%.6f,%.6f&zoom=%d&size=%dx%d&scale=%d&format=png&maptype=%s&key=%s&sensor=false', lat, lon, zoom, ev.width, ev.height, opt.scale, opt.type, ev.key);

            opturl = '';

            if ~opt.roadnames
                opturl = strcat(opturl, '&style=feature:road|element:labels|visibility:off');
            end
            if ~opt.placenames
                opturl = strcat(opturl, '&style=feature:administrative|element:labels.text|visibility:off&style=feature:poi|visibility:off');
            end

            if opt.onlyroads
                opturl = strcat(opturl, ['&style=feature:landscape|element:geometry.fill|color:0x000000|visibility:on'...
                    '&style=feature:landscape|element:labels|visibility:off'...
                    '&style=feature:administrative|visibility:off'...
                    '&style=feature:road|element:geometry.fill|color:0xffffff|visibility:on'...
                    '&style=feature:road|element:labels|visibility:off'...
                    '&style=feature:poi|element:all|visibility:off'...
                    '&style=feature:transit|element:all|visibility:off'...
                    '&style=feature:water|element:all|visibility:off'...
                    ]);
            end

            %opturl = urlencode(opturl)
            opturl = strrep(opturl, '|', '%7C');

            %opturl

            try
            [idx,cmap] = imread([baseurl opturl], 'png');
            catch
                error('Failed to read the URL: check your Google API key');
            end
            cmap = iint(cmap);

            % apply the color map
            view = cmap(idx(:)+1,:);

            % knock it into shape
            view = shiftdim( reshape(view', [3 size(idx)]), 1);
            view = ev.convert(view);

            if nargout == 0
                idisp(view);
            else
                im = view;
            end

            if nargout > 1
                % compute the northing/easting at each pixel.
                %
                % the following lines of code from get_google_map by Val Schmidt
                % ESTIMATE BOUNDS OF IMAGE:
                %
                % Normally one must specify a center (lat,lon) and zoom level for a map.
                % Zoom Notes:
                % ZL: 15, hxw = 640x640, image dim: 2224.91 x 2224.21 (mean 2224.56)
                % ZL: 16, hxw = 640x640, image dim: 1128.01m x 1111.25m (mean 1119.63)
                % This gives an equation of roughly (ZL-15)*3.4759 m/pix * pixels
                % So for an image at ZL 16, the LHS bound is
                % LHS = centerlonineastings - (zl-15) * 3.4759 * 640/2;
                [lonutm latutm zone] = deg2utm(lat,lon);
                Hdim = (2^15/2^zoom) * 3.4759 * ev.width;
                Vdim = (2^15/2^zoom) * 3.4759 * ev.height;

                ell = lonutm - Hdim/2;
                nll = latutm - Vdim/2;
                eur = lonutm + Hdim/2;
                nur = latutm + Vdim/2;

                Nvec = linspace(nur,nll,ev.height); % lat is highest at image row 1
                Evec = linspace(ell,eur,ev.width);
                [E,N] = meshgrid(Evec, Nvec);
            end
        end

        function close()
        end

        function s = char(m)
        %EarthView.char Convert to string
        %
        % EV.char() is a string representing the state of the EarthView object in
        % human readable form.
        %
        % See also EarthView.display.

            s = sprintf('EarthView: %d x %d', ...
                m.width, m.height);
            s2 = char@ImageSource(m);
            if ~isempty(s2)
                s = strvcat(s, strcat(' - ', s2));
            end
        end

        function paramSet(varargin)
        end

    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% *************************************************************************
% Video Super-Resolution with Convolutional Neural Networks
%
% This function calculates the motion compensated, bicubically upsampled
% input frames
%
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         02/19/2016
%
% *************************************************************************
function [input_data,idx_gt] = preprocess_frames(im_LR,upscale,MOTIONCOMPENSATION,ADAPTIVEMOTIONCOMPENSATION)

%% parameters
DParam.prescaling = 1;
DParam.upscaleFactor = upscale;
DParam.patchSize = 36;
DParam.stride = 36;
framesPerBlock = 5;

if ADAPTIVEMOTIONCOMPENSATION
    MOTIONCOMPENSATION = 1;
    ADAPT_WEIGHT = 8;
else
    ADAPT_WEIGHT = 0;
end

%% initialization
image_dims = [size(im_LR,1) size(im_LR,2)]*DParam.upscaleFactor;
nrFrames = size(im_LR,3);
framesOrder = 1:framesPerBlock;
framesOrder(floor(framesPerBlock/2)+1) = [];
framesOrder = [floor(framesPerBlock/2)+1 framesOrder];

nrSets = nrFrames-framesPerBlock+1;
tmpImgH5 = zeros([image_dims(2) image_dims(1) framesPerBlock+2 nrSets],'single');

%% loop through framesets (one frameset consists of 5 frames)
for i=1:nrSets
    display(['preprocessing frame ' num2str(i) ' of ' num2str(nrSets)])
    tic;
    %% loop through frames
    imgHi_bic_center = [];
    for layer=framesOrder

        %% load image
        fileIdx = i+layer-1;
        imgLo = im_LR(:,:,fileIdx);
        if isinteger(imgLo)
            imgLo = im2double(imgLo);
        end

        imgHi_bic = imresize(imgLo,DParam.upscaleFactor,'bicubic');

        %% motion compensation
        if MOTIONCOMPENSATION
            if layer == floor(framesPerBlock/2)+1 %center frame
                imgHi_bic_center = imgHi_bic;
            else

                [~,mc_err,~,~,imgHi_mc]=track_CLG_TV_adv(imgHi_bic_center,imgHi_bic,DParam,DParam);

                if ADAPT_WEIGHT > 0 %adaptive motion compensation
                    ratio_mc = exp(-mc_err./ADAPT_WEIGHT);
                    ratio_bic = ones(size(imgHi_bic)) - ratio_mc;

                    imgHi_bic = ratio_bic.*imgHi_bic_center + ratio_mc.*imgHi_mc;

                else
                    imgHi_bic = imgHi_mc;
                end
            end
        end

        %% convert into caffe-shape
        Xl_tmp = imgHi_bic;

        Xl_tmp = permute(Xl_tmp,[2 1 3 4]);
        tmpImgH5(:,:,layer,i) = Xl_tmp;

        %% calculate temporal image gradients -> not used for VSRnet
        if layer == framesPerBlock
            %first derriv.
            tmpImgH5(:,:,6,i) = tmpImgH5(:,:,2,i) - tmpImgH5(:,:,4,i);
            %second derriv.g
            tmpImgH5(:,:,7,i) = (tmpImgH5(:,:,1,i) + tmpImgH5(:,:,5,i) - 2*tmpImgH5(:,:,3,i));
        end

    end
    toc;
end

%% return values
input_data = tmpImgH5;
input_data = {input_data};
idx_gt = (1:nrSets) + floor(framesPerBlock/2);

end

author: thomas-koehler

repo: SupER


text:
function jnii=jnifticreate(varargin)
%
%    jnii=jnifticreate
%       or
%    jnii=jnifticreate('header1', value1, 'header2', value2, ...)
%    jnii=jnifticreate(img, 'header1', value1, ...)
%
%    Create a default JNIfTI structure with default header and image volume
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        img: set the jnii.NIFTIData section
%        'header_i': the header subfield name defined in the JNIfTI
%                    specification, see https://github.com/fangq/jnifti
%        value_i: set the value for the specified JNIfTI header field
%
%    output:
%        jnii: without any input, jnii gives the default jnii header
%              if img is given, jnii also includes the NIFTIData field
%
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

jnii=struct(encodevarname('_DataInfo_'),struct(),'NIFTIHeader',struct(), 'NIFTIData', []);

% jnii.NIFTIHeader.NIIHeaderSize=  0;
% jnii.NIFTIHeader.A75DataTypeName=   'uint8';
% jnii.NIFTIHeader.A75DBName=      '';
% jnii.NIFTIHeader.A75Extends=     0;
% jnii.NIFTIHeader.A75SessionError='';
% jnii.NIFTIHeader.A75Regular=     0;
jnii.NIFTIHeader.DimInfo.Freq=   0;
jnii.NIFTIHeader.DimInfo.Phase=  0;
jnii.NIFTIHeader.DimInfo.Slice=  0;
jnii.NIFTIHeader.Dim=            [];
jnii.NIFTIHeader.Param1=         0;
jnii.NIFTIHeader.Param2=         0;
jnii.NIFTIHeader.Param3=         0;
jnii.NIFTIHeader.Intent=         '';
jnii.NIFTIHeader.DataType=       'uint8';
jnii.NIFTIHeader.BitDepth=       8;
jnii.NIFTIHeader.FirstSliceID=   1;
jnii.NIFTIHeader.VoxelSize=      [1,1,1,1];
jnii.NIFTIHeader.Orientation=    struct('x','r','y','a','z','s');
% jnii.NIFTIHeader.NIIByteOffset=  0;
jnii.NIFTIHeader.ScaleSlope=     0;
jnii.NIFTIHeader.ScaleOffset=    0;
jnii.NIFTIHeader.LastSliceID=    1;
jnii.NIFTIHeader.SliceType=      '';
jnii.NIFTIHeader.Unit=           struct('L','mm','T', 's');
jnii.NIFTIHeader.MaxIntensity=   255;
jnii.NIFTIHeader.MinIntensity=   0;
jnii.NIFTIHeader.SliceTime=      1;
jnii.NIFTIHeader.TimeOffset=     0;
% jnii.NIFTIHeader.A75GlobalMax=  255;
% jnii.NIFTIHeader.A75GlobalMin=   0;
jnii.NIFTIHeader.Description=    '';
% jnii.NIFTIHeader.AuxFile=        '';
jnii.NIFTIHeader.QForm=          0;
jnii.NIFTIHeader.SForm=          1;
jnii.NIFTIHeader.Quatern.b=      0;
jnii.NIFTIHeader.Quatern.c=      0;
jnii.NIFTIHeader.Quatern.d=      0;
jnii.NIFTIHeader.QuaternOffset.x=0;
jnii.NIFTIHeader.QuaternOffset.y=0;
jnii.NIFTIHeader.QuaternOffset.z=0;
jnii.NIFTIHeader.Affine(1,:)=    [1 0 0 0];
jnii.NIFTIHeader.Affine(2,:)=    [0 1 0 0];
jnii.NIFTIHeader.Affine(3,:)=    [0 0 1 0];
jnii.NIFTIHeader.Name=           'default';
jnii.NIFTIHeader.NIIFormat=      'jnifti';
% jnii.NIFTIHeader.NIIExtender=    [0,0,0,0];

datainfo.JNIFTIVersion='0.5';
datainfo.Comment='Created by JNIFTY Toolbox (https://github.com/NeuroJSON/jnifty)';
datainfo.AnnotationFormat='https://github.com/NeuroJSON/jnifti/blob/master/JNIfTI_specification.md';
datainfo.SerialFormat='http://json.org';
datainfo.Parser=struct('Python',[], ...
                       'MATLAB',[], ...
                       'JavaScript', 'https://github.com/NeuroJSON/jsdata',...
                       'CPP', 'https://github.com/NeuroJSON/json',...
                       'C', 'https://github.com/NeuroJSON/ubj');
datainfo.Parser.Python={'https://pypi.org/project/jdata','https://pypi.org/project/bjdata'};
datainfo.Parser.MATLAB={'https://github.com/NeuroJSON/jnifty','https://github.com/NeuroJSON/jsonlab'};
datainfo.JNIFTIVersion='0.5';
jnii.(encodevarname('_DataInfo_'))=datainfo;

if(nargin==0)
    return;
end

img=[];
pid=1;
if(~ischar(varargin{1}))
    img=varargin{1};
    pid=2;
end

if(~isempty(varargin))
     for i=pid:2:length(varargin)
         jnii.NIFTIHeader.(varargin{i})=varargin{i+1};
     end
end

if(~isnumeric(img) && ~islogical(img))
    error('img input must be a numerical or logical array');
end

jnii.NIFTIHeader.Dim=size(img);
jnii.NIFTIHeader.DataType=class(img);
info=whos('img');
jnii.NIFTIHeader.BitDepth=info.bytes/numel(img)*8;
jnii.NIFTIHeader.MinIntensity=min(img(:));
jnii.NIFTIHeader.MaxIntensity=max(img(:));

jnii.NIFTIData=img;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [output, fs] = sp_johansson_impl_n3(inp_fs, input, time_skew)
% Reconstruction of a signal from non-uniform samples
% Based on Sindhi and Prabhu's implementation of Johansson's
% method

% index_mapping is used to reorded filters in case the first input
% signal is actually delayed comparing to the second and not vice versa
N = length(time_skew); % number of samplers
fs = inp_fs * N; % full sampling frequency

T_inp = 1/inp_fs;

% index_mapping is used to reorded filters in case the first input
% signal is actually delayed comparing to the second and not vice versa
[time_skew, index_mapping] = sort(time_skew);
x11 = cell2mat(input')';
x11 = x11(index_mapping, :);

TQ = 1; % Nyquist sampling period

% Decimation Periods - in our case all ADCs sample with the same rate
T = [3*TQ 3*TQ 3*TQ];

K = 0.5*lcm(2*T(1), 2*T(2))/TQ; % number of samples in recurrent period
K = 0.5*lcm(2*K, 2*T(3))/TQ;
capT = K*TQ; % the full sampling period - of all samplers
M = capT./T;
capM = lcm(M(1), M(2));
capM = lcm(capM, M(3));
excess = ceil((K-1)/capM);

maxf = K/(excess*capM+1);
K1 = K/maxf;

ML = min(cellfun('length', input)); % number of slices
w_c = 0.85;

LF = capM*2*K1+1;  %% min length of LF should be capM*2*K1
n = -(LF-1)/2:1:(LF-1)/2;

taus = time_skew / T_inp * capT; % ADC delays in seconds
tausI = sort([taus(1) taus(2) taus(3)]);
tauI = zeros(K,ML);
for p = 1:K
    tauI(p,:) = tausI(p)+(0:ML-1)*capT;
end;

r = tausI-TQ*(0:K-1);
r = r.';
w_o = w_c*pi*TQ;
hJ = zeros(K,LF);
C = zeros(1,LF);
Nt = (LF-1)/2;
for i = 1:K
    C = -2*sin(w_o*(n-r(1+(mod(i-1-n,K)))'))./(pi*(n-r(1+(mod(i-1-n,K)))'));
    C(isnan(C))=-2*w_o/pi;
    C = C.';
    S = zeros(LF,LF);
    for k = 1:LF
        S(k,:) = sin(w_o*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')))./(pi*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')));
    end;
    S(isnan(S))=w_o/pi;
    hJ(i,:) = -0.5*S\C;
end;

x1 = reshape(x11,1,size(x11,2)*K);
y1 = zeros(K,length(x1));
for j=1:K
    y1(j,:) = filter(hJ(j,:),1,x1);
    y1(j,:) = upsample(downsample(y1(j,:),K,j-1),K)/K;
    y1(j,:) = filter([zeros(1,j-1),1],1,y1(j,:));
end;
y = K*0.25*sum(y1,1);
delayJ = (size(hJ,2)-1)/2;
y = real(y(1+delayJ:end));

output = y;
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function test_pull1643

% MEM 8gb
% WALLTIME 00:40:00
% DEPENDENCY ft_connectivity_plm ft_connectivityanalysis

% create some dummy data
ntrl = 100;
nsmp = 1000;
fs   = 1000;
nchan = 3;
for k = 1:ntrl
  trial{1,k} = randn(nchan,nsmp);
  time{1,k}  = (0:999)/fs;
end
for k = 1:nchan
  label{k,1} = sprintf('chan%02d',k);
end
data.trial = trial;
data.time  = time;
data.label = label;

cfg = [];
cfg.method = 'plm';
plm = ft_connectivityanalysis(cfg,data);

% So far the test code only 'checks' whether it runs through, no checks are
% done on the validity of the output structure (in terms of how it looks
% w.r.t. fields etc., nor on the validity of the numeric results)


text:
function cMap=grayjet(varargin)

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

ns=7;
cMap=[flipud(gray(ns)); jet(ns)];

[cMap]=resampleColormap(cMap,n);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function imgOut = bilateralFilterSI(img, imgEdges, sigma_s, sigma_r)
%
%		 imgOut = bilateralFilterSI(img, imgEdges, sigma_s, sigma_r)
%
%        This function implements a bilateral filter without
%        approximations. Note this function is very slow!
%
%		 Input:
%			-img: is an image to be filtered.
%           -imgEdges: is an edge image, where its edges will be transfered
%           to img. Note set this to [] if you do not want to transfer
%           edges.
%           -sigma_s: is the spatial sigma value.
%           -sigma_r: is the range sigma value.
%
%		 Output:
%			-imgOut: is the filtered image.
%
%     Copyright (C) 2014-18  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

output_class = class(img);

if ~isa(img, 'double')
    img = double(img);
end

if ~isa(imgEdges, 'double')
    imgEdges = double(imgEdges);
end

if ~isa(sigma_s, 'double')
    sigma_s = double(sigma_s);
end

if ~isa(sigma_r, 'double')
    sigma_r = double(sigma_r);
end

imgOut = bilateralFilterS(img, imgEdges, sigma_s, sigma_r);

if strcmp(output_class, class(imgOut)) == 0
    imgOut = cast(imgOut, output_class);
end

end

author: banterle


text:
function nitfHdr = readDESegSubhdr(nitfHdr,fid,numDEs)
%
% <<<<<<<<<< CLASSIFICATION: UNCLASSIFIED >>>>>>>>>>
%
%READDESEGSUBHDR Read and parse data extension segments

% Written by Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN1)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% DETERMINE NITF VERSION
nitfVer = nitfHdr.filehdr.FHDR(6:8);

%% VERSION 2.1

if strcmp(nitfVer ,'2.1')
    nitfHdr.(numDEs).DE = fread(fid,2,'uint8=>char')';
    nitfHdr.(numDEs).DESID = strtrim(fread(fid,25,'uint8=>char')');
    nitfHdr.(numDEs).DESVER = fread(fid,2,'uint8=>char')';
    nitfHdr.(numDEs).DECLAS = fread(fid,1,'uint8=>char')';
    nitfHdr.(numDEs).DESCLSY = fread(fid,2,'uint8=>char')';
    nitfHdr.(numDEs).DESCODE = strtrim(fread(fid,11,'uint8=>char')');
    nitfHdr.(numDEs).DESCTLH = strtrim(fread(fid,2,'uint8=>char')');
    nitfHdr.(numDEs).DESREL =  strtrim(fread(fid,20,'uint8=>char')');
    nitfHdr.(numDEs).DESDCTP =  strtrim(fread(fid,2,'uint8=>char')');
    nitfHdr.(numDEs).DESDCDT =  strtrim(fread(fid,8,'uint8=>char')');
    nitfHdr.(numDEs).DESDCXM =  strtrim(fread(fid,4,'uint8=>char')');
    nitfHdr.(numDEs).DESDG =  fread(fid,1,'uint8=>char')';
    nitfHdr.(numDEs).DESDGDT =  strtrim(fread(fid,8,'uint8=>char')');
    nitfHdr.(numDEs).DESCLTX =  strtrim(fread(fid,43,'uint8=>char')');
    nitfHdr.(numDEs).DESCATP =  fread(fid,1,'uint8=>char')';
    nitfHdr.(numDEs).DESCAUT = strtrim(fread(fid,40,'uint8=>char')');
    nitfHdr.(numDEs).DESCRSN = fread(fid,1,'uint8=>char')';
    nitfHdr.(numDEs).DESSRDT = strtrim(fread(fid,8,'uint8=>char')');
    nitfHdr.(numDEs).DESCTLN = strtrim(fread(fid,15,'uint8=>char')');

    if strcmp(nitfHdr.(numDEs).DESID,'TRE_OVERFLOW')
        nitfHdr.(numDEs).DESOFLW = fread(fid,6,'uint8=>char')';
        nitfHdr.(numDEs).DESITEM = str2double(fread(fid,3,'uint8=>char')');
    end

    nitfHdr.(numDEs).DESSHL = str2double(fread(fid,4,'uint8=>char')');

    if strcmp(nitfHdr.(numDEs).DESID,'TRE_OVERFLOW')
        if nitfHdr.(numDEs).DESSHL > 0
            nitfHdr.(numDEs).DESSHF = fread(fid,nitfHdr.(numDEs).DESSHL,'uint8=>char')';
        end
        TRE = fread(fid,6,'uint8=>char')';
        fseek(fid, -6, 'cof');
        % Check to see if reader function exists for this TRE.
        if exist(strcat('read',TRE,'.m'),'file')
            nitfHdr.(numDEs) = eval(strcat('read',TRE,'(fid)'));
        else
            % Read undefined data
            nitfHdr.(numDEs) = readUndefinedTRE(fid);
        end
    elseif strcmp(nitfHdr.(numDEs).DESID,'STREAMING_FILE_HEADER')
        nitfHdr.(numDEs).SFH_L1 = str2double(fread(fid,7,'uint8=>char')');
        nitfHdr.(numDEs).SFH_DELIM1 = fread(fid,4,'uint8=>char')';
        nitfHdr.(numDEs).SFH_DR = str2double(fread(fid,nitfHdr.(numDEs).SFH_L1,'uint8=>char')');
        nitfHdr.(numDEs).SFH_DELIM2 = fread(fid,4,'uint8=>char')';
        nitfHdr.(numDEs).SFH_L2 = str2double(fread(fid,7,'uint8=>char')');
    elseif strcmp(nitfHdr.(numDEs).DESID,'SICD_XML')
        SICD_xml_string = fread(fid,nitfHdr.filehdr.LD(str2double(numDEs(10:end))),'uint8=>char')';
        nitfHdr.(numDEs).DESDATA = sicdxml2struct( xmlread( java.io.StringBufferInputStream( SICD_xml_string ) ) );
    else % Skip unrecognized DES types
        nitfHdr.(numDEs).DESUSH_OFFSET = ftell(fid); % Record offset, so we can parse user subheaders later
        fseek(fid, nitfHdr.(numDEs).DESSHL, 'cof');
        fseek(fid, nitfHdr.filehdr.LD(str2double(numDEs(10:end))), 'cof');
    end

%% VERSION 2.0

elseif strcmp(nitfVer, '2.0')
    nitfHdr.(numDEs).DE = fread(fid,2,'uint8=>char')';
    nitfHdr.(numDEs).DESTAG = strtrim(fread(fid,25,'uint8=>char')');
    nitfHdr.(numDEs).DESVER = fread(fid,2,'uint8=>char')';

    % nitfHdr.(numDEs).DESSG
    nitfHdr.(numDEs).DESSG.DESCLAS = fread(fid,1,'uint8=>char')';
    nitfHdr.(numDEs).DESSG.DESCODE = strtrim(fread(fid,40,'uint8=>char')');
    nitfHdr.(numDEs).DESSG.DESCTLH = strtrim(fread(fid,40,'uint8=>char')');
    nitfHdr.(numDEs).DESSG.DESREL =  strtrim(fread(fid,40,'uint8=>char')');
    nitfHdr.(numDEs).DESSG.DESCAUT = strtrim(fread(fid,20,'uint8=>char')');
    nitfHdr.(numDEs).DESSG.DESCTLN = strtrim(fread(fid,20,'uint8=>char')');
    nitfHdr.(numDEs).DESSG.DESDWNG = strtrim(fread(fid,6,'uint8=>char')');
    if strcmp(nitfHdr.(numDEs).DESSG.DESDWNG,'999998')
        nitfHdr.(numDEs).DESSG.DESDEVT = strtrim(fread(fid,40,'uint8=>char')');
    end

    if strcmpi(nitfHdr.(numDEs).DESTAG,'Registered Extensions') ||...
            strcmpi(nitfHdr.(numDEs).DESTAG,'Controlled Extensions')
        nitfHdr.(numDEs).DESOFLW = strtrim(fread(fid,6,'uint8=>char')');
        nitfHdr.(numDEs).DESITEM = fread(fid,3,'uint8=>char')';
    end

    nitfHdr.(numDEs).DESSHL = str2double(fread(fid,4,'uint8=>char')');
    if nitfHdr.(numDEs).DESSHL > 0
        nitfHdr.(numDEs).DESSHF = fread(fid,nitfHdr.(numDEs).DESSHL,'uint8=>char')';
    end

    if strcmpi(nitfHdr.(numDEs).DESTAG,'Registered Extensions') ||...
            strcmpi(nitfHdr.(numDEs).DESTAG,'Controlled Extensions')
        % READ TRE
        TRE = fread(fid,6,'uint8=>char')';
        fseek(fid, -6, 'cof');

        % Check to see if reader function exists for this TRE.
        if exist(strcat('read',TRE,'.m'),'file')
            nitfHdr.(numDEs) = eval(strcat('read',TRE,'(fid)'));
        else
            % Read undefined data
            nitfHdr.(numDEs) = readUndefinedTRE(fid);
        end
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [TmodelC, Cspawn, Stats] = mergeModelsBorg(CmodelIn, TmodelIn, rxnList, metList, Stats, score, mode, varargin)
% Checks Tmodel for duplicate reactions and other mistakes,
% that may have occured during reaction and metabolite matching. It
% resolves these problems and merges the models, and confirms that Cmodel
% is the same after being removed from the merged model. It also provides
% some statistics on the merging process and resulting combined model.
% Called by `driveModelBorgifier`, calls `addToTmodel`, `compareMatricies`, `reactionCompare`, `cleanTmodel`, organizeModelCool`, `TmodelStats`, `readCbTmodel`.
%
% USAGE:
%
%    [TmodelC, Cspawn, Stats, CMODEL] = mergeModelsBorg(CmodelIn, TmodelIn, rxnList, metList, Stats, score)
%
% INPUTS:
%    CmodelIn:     Comparison model
%    TmodelIn:     Template model
%    rxnList:      Array which designates matched and new reactions.
%    metList:      Array which desginates matched and new metabolites.
%    Stats:        Structure that comes from reactionCompare. Weighting
%                  information can be used and additional information addended.
% OPTIONAL INPUTS:
%    score:        The originalPress the any key to continue.
 scoring matrix, which may be used to correct
%                  problematic reaction upon recomparison.
%    mode:         {('p'),'a'} 'Revisit only [p]roblematic reactions or [a]ll
%                  reactions that problematic metabolites are involved in.
%    'Verbose':    Print statements on progress.
%
%
% OUTPUTS:
%    TmodelC:      Combined `C` and `Tmodel`.
%    Cspawn:       Cmodel after it has been removed from the `TmodelC`
%    Stats:        Structure of information regarding the merging and also stats
%                  on the combined model.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

if exist('score','var')
    global SCORE CMODEL TMODEL % Declare variables.
    SCORE = score ;
else
    global SCORE CMODEL TMODEL % Declare variables.
end

if ~exist('mode','var')
    modeProvided=0;
    mode='p';
else
    modeProvided=1;
end
CMODEL = CmodelIn ;
TMODEL = TmodelIn ;
% SCORE = score ;
TmodelC = [] ;
Cspawn = [] ;

verbose = false ;
if ~isempty(varargin)
    if sum(strcmp('Verbose', varargin))
        verbose = true ;
    end
end

% Allow pausing for some of the UI.
% pause on

% set flags for the two checks.
reviewMets = 1 ;
checkSimilarity = 1 ;

%% Check that metList is okay.
% Make sure all mets have been declared and that 2 mets from C aren't
% matched to the same met in T. Continue until everything is a-okay.

while reviewMets
    % Find mets that have not been declared.
    reviewMets = find(metList == 0) ;

    % Find mets in C that have been matched to the same metabolite in T.
    uniqMetIndex = unique(metList) ;
    count = histc(metList, uniqMetIndex) ;
    dups = uniqMetIndex(count > 1) ;
    for iDup = 1:length(dups)
        metIndex = find(metList == dups(iDup)) ;
        reviewMets(end + 1:length(reviewMets) + length(metIndex)) = metIndex ;
    end

    % Review all above mets with GUI.
    if ~isempty(reviewMets)
        fprintf('Problems within metList, resolve with GUI.\n')
        reviewMetsAgain = input('Press c to continue anyways or any other key to re-check metabolite matching.\n', 's') ;
        if ~strcmpi(reviewMetsAgain, 'c')
            RxnInfo.rxnIndex = 0 ; % Tells GUI just look at mets.
            RxnInfo.rxnList = rxnList ;
            RxnInfo.metList = metList ;
            RxnInfo.metIndex = reviewMets ;
            metList = metCompare(RxnInfo) ;
        else
            reviewMets = 0 ; % just go on without dealing with problem.
            checkSimilarity = 0 ; % matricies will be messed up, so do not do checking below
            fprintf('Skipped resolving, will not check fidelity of matricies.\n') ;
        end
    end
end

% Force good numbering for new metabolites
TmetNum = length(TMODEL.mets) ;
metList(metList > TmetNum) = TmetNum+(1:sum(metList > TmetNum)) ;

% merge and extract C model
CMODELoriginal = CMODEL ;
if verbose
    [TmodelC, CMODEL] = addToTmodel(CMODEL, TMODEL, rxnList, metList, 'NoClean', 'Verbose') ;
else
    [TmodelC, CMODEL] = addToTmodel(CMODEL, TMODEL, rxnList, metList, 'NoClean') ;
end

% Create Cmodel again from the combined model.
if verbose
    Cspawn = readCbTmodel(CMODEL.description, TmodelC, 'y', 'Verbose') ;
else
    Cspawn = readCbTmodel(CMODEL.description, TmodelC, 'y') ;
end

%% Reconsider reactions that don't have the same stoich between Cmodel and Cspawn
try
while checkSimilarity
    % Combine models and create spawn of Cmodel.
    % Combine the models. Problematic mets must have been resolved
    FluxCompare = compareMatricies(CMODELoriginal, Cspawn) ;

    % If there are differences pause and let the user know what is up.
    if sum(sum(FluxCompare.diffS)) || (sum(rxnList == -1) > 0)
        % Plots the differences between the S matricies.
        figure
        subplot(2, 2, 1)
        spy(FluxCompare.CmodelS)
        title([CMODEL.description ' before merging.'])
        subplot(2, 2, 2)
        spy(FluxCompare.CspawnS)
        title([CMODEL.description ' after merging.'])
        subplot(2, 1, 2)
        spy(FluxCompare.diffS)
        title('Difference.')

        % Pause and tell the users what is happening.
        fprintf('Difference between the matrices. Reactions from C\n')
        fprintf('that do not have the same stoich before and after\n')
        fprintf('merging need to be reviewed again.\n')
        fprintf('Press the any key to continue.\n')

        % Find the wrong reactions and metabolies,
        % mark them as undeclared, call GUI, and loop.
        diffmets = logical(sum(abs(FluxCompare.diffS), 2)) ;
        metList(FluxCompare.CmetsSorti(diffmets)) = 0 ;
        if ~modeProvided
            fprintf('Revisit only [p]roblematic reactions or [a]ll reactions that problematic metabolites are involved in?\n')
            mode = input('default = p ', 's') ;
        end
        if strcmp(mode, 'a')
            diffrxns = logical(sum(abs(FluxCompare.diffS), 1) + sum(FluxCompare.CmodelS(diffmets, :)) + sum(FluxCompare.CspawnS(diffmets, :))) ;
        else
            diffrxns = logical(sum(abs(FluxCompare.diffS), 1)) ;
        end
        rxnList(FluxCompare.CrxnsSorti(diffrxns)) = -1 ;

        [rxnList, metList, Stats] = reactionCompare(CMODEL, TMODEL, SCORE, rxnList, metList, Stats) ;

        % remerge and extract C model for second round of checking.
        [TmodelC, CMODEL] = addToTmodel(CMODEL, TMODEL, rxnList, metList, 'NoClean') ;
        Cspawn = readCbTmodel(CMODEL.description, TmodelC, 'y') ;
    else
        % Set the flag to not check for differences again.
        fprintf('Matricies are now equal before and after merging.\n')
        % Turn off flag.
        checkSimilarity = 0 ;
        break
    end
end
catch % return safely if re-matching is aborted by user
    return
end

% Add final rxn and metList to Stats.
Stats.rxnList = rxnList ;
Stats.metList = metList ;

%% Organize and Get stats on combined Tmodel.
if verbose
    TmodelC = cleanTmodel(TmodelC, 'Verbose') ;
else
    TmodelC = cleanTmodel(TmodelC) ;
end
try
    TmodelC = organizeModelCool(TmodelC) ;
catch
    disp('TmodelC too large for cool organization of stoichiometric matrix.')
end
if verbose
    Cspawn = readCbTmodel(CMODEL.description, TmodelC, 'y', 'Verbose') ;
else
    Cspawn = readCbTmodel(CMODEL.description, TmodelC, 'y') ;
end
Stats = TmodelStats(TmodelC, Stats) ;
end


text:
% ATAN   Inverse tangent, result in radians.
%    ATAN(X) is the arctangent of the elements of X.
%
%    See also ATAN2, TAN, ATAND, ATAN2D.
%
%    Reference page in Doc Center
%       doc atan
%
%    Other functions named atan
%
%       codistributed/atan    gpuArray/atan    sym/atan    ts/atan
%

author: jmaih


text:
function fx2 = p08_fx2 ( x )

%*****************************************************************************80
%
%% P08_FX2 evaluates the second derivative of the function for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  fx2 = - cos ( x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = subsasgn( x, S, y )
error( 'Subscripted assignment not allowed for cvx problem objects.' );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function varargout = OceanMap(varargin)

%OceanMap - GUI: Manually input values of an m x n matrix and save as *.MAT file.
%
% Syntax:
%       OCEANMAP, by itself, creates a new OCEANMAP or raises the existing
%       singleton*.
%
%       H = OCEANMAP returns the handle to a new OCEANMAP or the handle to
%       the existing singleton*.
%
%       OCEANMAP('CALLBACK',hObject,eventData,handles,...) calls the local
%       function named CALLBACK in OCEANMAP.M with the given input arguments.
%
%       OCEANMAP('Property','Value',...) creates a new OCEANMAP or raises the
%       existing singleton*.  Starting from the left, property value pairs are
%       applied to the GUI before OceanMap_OpeningFunction gets called.  An
%       unrecognized property name or invalid value makes property application
%       stop.  All inputs are passed to OceanMap_OpeningFcn via varargin.
%
%       *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%       instance to run (singleton)".
%
% Inputs:
%       none
%
% Outputs:
%       saved file (*.MAT) containing m x n array
%
% Example:
%       Not really applicable- load the GUI and see!
%
% Notes:
%       The number of rows and columns can be updated actively, without
%       losing data already contained in the matrix
%
%       You should call OceanMap from the command line, as clicking on the
%       .fig file doesnt work: The ActiveX FlexArray doesn't initialise.
%
%       Currently, Multidimensional Arrays are not supported... I'd be
%       pleased to hear if anyone modifies the software for extra
%       dimensions.
%
%       The output matrix DOES NOT CONTAIN row and column numbers
%
% Background Information:
%       OceanMap was originally created for the purpose of tracking a
%       submersible vehicle using a grid painted onto the floor of a water
%       tank. Technicians use this software to input the grid pattern
%       painted and convert it to matrix form. Other software performed
%       postprocessing to find the submarine's path by comparing grid
%       numbers (recorded by the sub's onboard camera) and comparing them
%       to the output matrix from OceanMap.
%
%       Clearly, this software can be used for far less specialised purposes.

% Other m-files required:
%       none
%
% Subfunctions:
%       All create-fcns and callbacks etc associated with GUI Objects (see Help GUIDE)
%       labelupdate(handles)
%                           - refreshes row/column numbering
%       savethedamnfile(filename, savethismatrix)
%                           - saves file according to filename input in GUI
%
% Files required:
%       OceanMap.m OceanMap.fig OceanMap_activex1
%
% See also:
%       GridNav GUI file (for postprocessing of results)
%       GUIDE, GUIDATA, GUIHANDLES
%
% Author: Thomas Clark, B.A.(Cantab) General Engineering
% email address: thc29@cam.ac.uk
% Created August 2005; Last revision: 25-August-2005
% Using Matlab 7.0.1.24704 (R14) Service Pack 1
%
% Thanks to Denis Gilbert (Maurice Lamontagne Institute, Dept. of Fisheries
% and Oceans Canada) for header file format (from www.mathworks.com).
%
% GUI Initialisation code Copyright 2002-2003 The MathWorks, Inc.

% --------------------------BEGIN CODE-------------------------------------------------------BEGIN CODE

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @OceanMap_OpeningFcn, ...
                   'gui_OutputFcn',  @OceanMap_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

% -------------------------OPENING FUNCTION----------------------------------------OPENING FUNCTION (inc. activex1)

% --- Executes just before OceanMap is made visible.-

function OceanMap_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 OceanMap (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% Initialising row and column numbers
labelupdate(handles);

% initialising global variables
global NUMROWS
NUMROWS = 200;
global NUMCOLS
NUMCOLS = 200;

% ------------------LABELLING UPDATE FUNCTION-------------------------------------------LABELLING UPDATE FUNCTION
function labelupdate(handles)

% Label the rows
set(handles.activex1,'Col',0);
for k=1:get(handles.activex1,'Rows')-1
    set(handles.activex1,'Row',k);
    set(handles.activex1,'Text',['Row ' num2str(k)]);
end;

% Label the Columns
k = 0;
set(handles.activex1,'Row',0);
set(handles.activex1,'Col',0);
set(handles.activex1,'Text','Column #:');
for k = 1:get(handles.activex1,'Cols')-1
    set (handles.activex1,'Col',k);
    set (handles.activex1,'Text',[num2str(k)]);
end;

% --- Outputs from this function are returned to the command line.
function varargout = OceanMap_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;

% ------------------FILENAME TEXT BOX 'EDIT 1'-------------------------------------------FILENAME EDIT TEXT BOX

function fileditbox_Callback(hObject, eventdata, handles)
% hObject    handle to fileditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of fileditbox as text
%        str2double(get(hObject,'String')) returns contents of fileditbox as a double

global FILENAME;
FILENAME = get(hObject,'String');

% --- Executes during object creation, after setting all properties.
function fileditbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to fileditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

global FILENAME;
FILENAME = get(hObject,'String');

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), ...
                                get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -------------------DATA EXPORT-----------------------------------------------------------------DATA EXPORT

% --- Executes on button press in exportbutton.
function exportbutton_Callback(hObject, eventdata, handles)
% hObject    handle to exportbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Recalling the desired filename in double format from the Edit1 textbox
global FILENAME;

% Initialising the OceanMap storage matrix:
gridmatrix = zeros([get(handles.activex1,'rows')-1],...
                                    [get(handles.activex1,'Cols')-1]);

% Obtaining and storing the OceanMap pattern (excluding row and col #s).
%   N.B. Null entries returned as 0, Invalid Entries as NaN
%        No imaginary numbers used. i,j, are scratch variables.

for i = 2:get(handles.activex1,'Cols') % For Loop 1
    for j = 2:get(handles.activex1,'rows') % For Loop 2
        set(handles.activex1,'Col',(i-1));
        set(handles.activex1,'Row',(j-1));
        gridmatrix((j-1),(i-1)) =  get(handles.activex1,'Value');
    end; % End For 2
end; %End For 1

savethefile(FILENAME,gridmatrix);

% Saving the variable gridmatrix in a .m file corresponding to the input
% filename.
function savethefile(filename,savethismatrix)

save(filename, 'savethismatrix');

% -------------------------OPTIONS PANEL------------------------------------------------------OPTIONS PANEL

function rowseditbox_Callback(hObject, eventdata, handles)
% hObject    handle to rowseditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of rowseditbox as text
%        str2double(get(hObject,'String')) returns contents of rowseditbox as a double
global NUMROWS
NUMROWS = str2num(get(hObject,'String'))+1;

% --- Executes during object creation, after setting all properties.
function rowseditbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to rowseditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'),...
                                get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function colseditbox_Callback(hObject, eventdata, handles)
% hObject    handle to colseditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of colseditbox as text
%        str2double(get(hObject,'String')) returns contents of colseditbox as a double
global NUMCOLS
NUMCOLS = str2num(get(hObject,'String'))+1;

% --- Executes during object creation, after setting all properties.
function colseditbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to colseditbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'),...
                                get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in updatebutton.
function updatebutton_Callback(hObject, eventdata, handles)
% hObject    handle to updatebutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Update the number of columns and rows in the table
global NUMCOLS
global NUMROWS
set(handles.activex1,'Rows', NUMROWS);
set(handles.activex1,'Cols', NUMCOLS);
labelupdate(handles);

% -------------------------END OF CODE-------------------------------------------------------END OF CODE

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over exportbutton.
function exportbutton_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to exportbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


text:
function ediv = q2div(xy,mv,flowsol)
%q2div  computes norm of divergence of Q2 flow solution
%   ediv = q2div(xy,mv,flowsol);
%   input
%          xy         vertex coordinate vector
%          mv         Q2 element mapping matrix
%          flowsol    Q2 solution vector
%
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      x=xy(:,1); y=xy(:,2); nvtx=length(x);
      nel=length(mv(:,1));
	  usol=flowsol(1:nvtx); vsol=flowsol(nvtx+1:2*nvtx);
      fprintf('computing divergence of discrete velocity solution ...  ')

%
% initialise global matrices
       ediv=zeros(nel,1);
%
% set up 3x3 Gauss points
      gpt=sqrt(0.6);
      s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
      s(2) =  gpt; t(2) = -gpt; wt(2)=25/81;
      s(3) =  gpt; t(3) =  gpt; wt(3)=25/81;
      s(4) = -gpt; t(4) =  gpt; wt(4)=25/81;
      s(5) =  0.0; t(5) = -gpt; wt(5)=40/81;
      s(6) =  gpt; t(6) =  0.0; wt(6)=40/81;
      s(7) =  0.0; t(7) =  gpt; wt(7)=40/81;
      s(8) = -gpt; t(8) =  0.0; wt(8)=40/81;
      s(9) =  0.0; t(9) =  0.0; wt(9)=64/81;
%
%
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(mv(:,ivtx));
        yl_v(:,ivtx) = y(mv(:,ivtx));
        end
        for idx = 1:9
		xsl(:,idx) = usol(mv(:,idx));
		ysl(:,idx) = vsol(mv(:,idx));
		end
% loop over 3x3 Gauss points
        for igpt = 1:9
        sigpt=s(igpt);
        tigpt=t(igpt);
        wght=wt(igpt);
%  evaluate derivatives etc
         [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
         [psi,dpsidx,dpsidy] = qderiv(sigpt,tigpt,xl_v,yl_v);
		 evec=zeros(nel,1);
            for j = 1:9
	evec(:) = evec(:) + wght * (xsl(:,j) .* dpsidx(:,j) + ysl(:,j) .* dpsidy(:,j));
		    end
	     ediv(:)=ediv(:) + evec(:).*evec(:);
% end of Gauss point loop
         end
%
% end of element loop
%
      err_div = sqrt(sum(ediv)); ediv = sqrt(ediv);
fprintf('done\n')
fprintf('estimated velocity divergence error:  %10.6e \n',err_div)
return


text:
function S = sub(A, i, j)
%SUB     Principal submatrix.
%        SUB(A,i,j) is A(i:j,i:j).
%        SUB(A,i)  is the leading principal submatrix of order i,
%        A(1:i,1:i), if i>0, and the trailing principal submatrix
%        of order ABS(i) if i<0.

if nargin == 2
   if i >= 0
      S = A(1:i, 1:i);
   else
      n = min(size(A));
      S = A(n+i+1:n, n+i+1:n);
   end
else
   S = A(i:j, i:j);
end


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function y_hat = useLR(reg_model,X)
%y_hat = useLR(reg_model,X)
%   X: m * n
%   reg_model: (n+1)*ny
%   y: m * ny

m = size(X,1);
X = [X ones(m,1)];
y_hat = X * reg_model;

end

author: zhusz


text:
function [ Aout ] = estimate_airlight( img, Amin, Amax, N, spacing, K, thres )
%Estimate airlight of an image, using a 3*2D Hough transform, where each
%point votes for a given location using a fixed set of angles.
%
%   This is an implementation of our paper:
%   Dana Berman, Tali Treibitz, Shai Avidan,
%   "Air-light Estimation using Haze-Lines", ICCP 2017
%
%   Input arguments:
%   img - input image (mandatory)
%   Amin (Amax) - minimal (maximal) value if the air-light. Optional
%               Can be used to reduce the search space and save time
%               Either a scalar (identical for all color channels) or a 3
%               value vector (different range for each color channel)
%   N - number of colors clusters to use (the image is converted to an
%       indexed image of at most N different cluster). Optional
%   spacing - air-light candidates' resolution, 1/M in the paper. Optional
%   K - angular resolution. Optional
%   thres - cone resolution, optional, default recommended

%% Verify input params, set defaults when necessary (same as published results)
if ~exist('thres','var') || isempty(thres), thres = 0.01 ; end;
if ~exist('spacing','var') || isempty(spacing), spacing = 0.02 ; end; %1/M in the paper
if ~exist('n_colors','var') || isempty(N), N = 1000 ; end; %number of colors clusters
if ~exist('K','var') || isempty(K), K = 40 ; end; %number of angles

% Define search range for the air-light. The search range is different for each
% color channel. These values were used in all of our experiments.
if ~exist('Amin','var') || isempty(Amin), Amin = [0,0.05,0.1]; end;
if ~exist('Amax','var') || isempty(Amax), Amax = 1; end;

% Air-light search range, accept a scalar if identical for all color channels
if isscalar(Amin), Amin = repmat(Amin,1,3); end
if isscalar(Amax), Amax = repmat(Amax,1,3); end

%% Convert input image to an indexed image
[img_ind, points] = rgb2ind(img, N);
[h,w,~] = size(img);
% Remove empty clusters
idx_in_use = unique(img_ind(:));
idx_to_remove = setdiff(0:(size(points,1)-1),idx_in_use);
points(idx_to_remove+1,:) = [];
img_ind_sequential = zeros(h,w);
for kk = 1:length(idx_in_use)
    img_ind_sequential(img_ind==idx_in_use(kk)) = kk;
end
% Now the min value of img_ind_sequential is 1 rather then 0, and the indices
% correspond to points

% Count the occurences if each index - this is the clusters' weight
[points_weight,~] = histcounts(img_ind_sequential(:),size(points,1));
points_weight = points_weight./(h*w);
if ~ismatrix(points), points = reshape(points,[],3); end % verify dim

%% Define arrays of candidate air-light values and angles
angle_list = reshape(linspace(0, pi, K),[],1);
% Use angle_list(1:end-1) since angle_list(end)==pi, which is the same line
% in 2D as since angle_list(1)==0
directions_all = [sin(angle_list(1:end-1)) , cos(angle_list(1:end-1)) ];

% Air-light candidates in each color channel
ArangeR = Amin(1):spacing:Amax(1);
ArangeG = Amin(2):spacing:Amax(2);
ArangeB = Amin(3):spacing:Amax(3);

%% Estimate air-light in each pair of color channels
% Estimate RG
Aall = generate_Avals(ArangeR, ArangeG);
[~, AvoteRG] = vote_2D(points(:,1:2), points_weight, directions_all, Aall, thres );
% Estimate GB
Aall = generate_Avals(ArangeG, ArangeB);
[~, AvoteGB] = vote_2D(points(:,2:3), points_weight, directions_all, Aall, thres );
% Estimate RB
Aall = generate_Avals(ArangeR, ArangeB);
[~, AvoteRB] = vote_2D(points(:,[1,3]), points_weight, directions_all, Aall, thres);

%% Find most probable airlight from marginal probabilities (2D arrays)
% Normalize (otherwise the numbers are quite large)
max_val = max( [max(AvoteRB(:)) , max(AvoteRG(:)) , max(AvoteGB(:)) ]);
AvoteRG2 = AvoteRG./max_val;
AvoteGB2 = AvoteGB./max_val;
AvoteRB2 = AvoteRB./max_val;
% Generate 3D volumes from 3 different 2D arrays
A11 = repmat( reshape(AvoteRG2, length(ArangeG),length(ArangeR))', 1,1,length(ArangeB));
tmp = reshape(AvoteRB2, length(ArangeB),length(ArangeR))';
A22 = repmat(reshape(tmp, length(ArangeR),1,length(ArangeB)) , 1,length(ArangeG),1);
tmp2 = reshape(AvoteGB2, length(ArangeB),length(ArangeG))';
A33 = repmat(reshape(tmp2, 1, length(ArangeG),length(ArangeB)) , length(ArangeR),1,1);
AvoteAll = A11.*A22.*A33;
[~, idx] = max(AvoteAll(:));
[idx_r,idx_g,idx_b] = ind2sub([length(ArangeR),length(ArangeG),length(ArangeB)],idx);
Aout = [ArangeR(idx_r), ArangeG(idx_g), ArangeB(idx_b)];

end % function estimate_airlight_2D

%% Sub functions

function Aall = generate_Avals(Avals1, Avals2)
%Generate a list of air-light candidates of 2-channels, using two lists of
%values in a single channel each
%Aall's length is length(Avals1)*length(Avals2)
Avals1 = reshape(Avals1,[],1);
Avals2 = reshape(Avals2,[],1);
A1 = kron(Avals1, ones(length(Avals2),1));
A2 = kron(ones(length(Avals1),1), Avals2);
Aall = [A1, A2];
end % function generate_Avals

function [Aout, Avote2] = vote_2D(points, points_weight, directions_all, Aall, thres)
n_directions = size(directions_all,1);
accumulator_votes_idx = false(size(Aall,1), size(points,1), n_directions);
for i_point = 1:size(points,1)
    for i_direction = 1:n_directions
		 % save time and ignore irelevant points from the get-go
        idx_to_use = find( (Aall(:, 1) > points(i_point, 1)) & (Aall(:, 2) > points(i_point, 2)));
        if isempty(idx_to_use), continue; end

        % calculate distance between all A options and the line defined by
        % i_point and i_direction. If the distance is smaller than a thres,
        % increase the cell in accumulator
        dist1 = sqrt(sum([Aall(idx_to_use, 1)-points(i_point, 1), Aall(idx_to_use, 2)-points(i_point, 2)].^2,2));
        %dist1 = dist1 - min(dist1);
        dist1 = dist1./sqrt(2) + 1;

        dist =  -points(i_point, 1)*directions_all(i_direction,2) + ...
            points(i_point, 2)*directions_all(i_direction,1) + ...
            Aall(idx_to_use, 1)*directions_all(i_direction,2) - ...
            Aall(idx_to_use, 2)*directions_all(i_direction,1);
        idx = abs(dist)<2*thres.*dist1;
        if ~any(idx), continue; end

        idx_full = idx_to_use(idx);
        accumulator_votes_idx(idx_full, i_point,i_direction) = true;
    end
end
% use only haze-lined that are supported by 2 points or more
accumulator_votes_idx2 = (sum(uint8(accumulator_votes_idx),2))>=2;
accumulator_votes_idx = bsxfun(@and, accumulator_votes_idx ,accumulator_votes_idx2);
accumulator_unique = zeros(size(Aall,1),1);
for iA = 1:size(Aall,1)
    idx_to_use = find(Aall(iA, 1) > points(:, 1) & (Aall(iA, 2) > points(:, 2)));
    points_dist = sqrt((Aall(iA,1) - points(idx_to_use,1)).^2+(Aall(iA,2) - points(idx_to_use,2)).^2);
    points_weight_dist = points_weight(idx_to_use).*(5.*exp(-reshape(points_dist,1,[]))+1);
    accumulator_unique(iA) = sum(points_weight_dist(any(accumulator_votes_idx(iA,idx_to_use,:),3)));
end
[~, Aestimate_idx] = max(accumulator_unique);
Aout = Aall(Aestimate_idx,:);
Avote2 = accumulator_unique;

end % function vote_2D

author: danaberman

repo: non-local-dehazing


text:
function plotRelativeAmps(view)
%
% plotAmps(view)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in the current ROI.
%
% gmb  5/25/98

% Compute means across scans, for all pixels in the
% currently selected ROI.
[meanAmps,meanPhs] = vectorMeans(view);

selectGraphWin

% Header
ROIname = view.ROIs(view.selectedROI).name;
headerStr = ['Mean of amplitudes, ROI ',ROIname];
set(gcf,'Name',headerStr);

% Plot the bar graph
clf
fontSize = 14;
mybar(meanAmps);
xlabel('Scan','FontSize',fontSize);
ylabel('Mean Amplitude','FontSize',fontSize);
ylim =get(gca,'YLim');
set(gca,'YLim',ylim*1.1);
set(gca,'FontSize',fontSize);

% Save the data in gca('UserData')
data.y = meanAmps;
set(gca,'UserData',data);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% sim_press_shaped.m
% Robin Simpson and Jamie Near, 2014.
%
% USAGE:
% out = sim_press_shaped(n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dx,dy,Gx,Gy,flipAngle,centreFreq)
%
% DESCRIPTION:
% This function simulates the PRESS experiment.  The excitation is
% simulated as an instantaneous rotation, and the refocusing pulses are
% simulated as shaped rotations.
%
% It employs coherence selection to only include desired coherence orders
% of the spin system being simulated. For the PRESS experiment, the desired
% coherence after the excitation pulse is -1, +1 after the first refocusing
% pulse and -1 after the final refocusing pulse.
%
% Finally, this code simulates the spectrum at a given point in space (x,y),
% given the values of the slice selection gradients (Gx, and Gy).  The pulse
% waveform is assumed to be the same for both refocusing pulses.  In order
% to fully simulate the MEGA-PRESS experiment, you have to run this
% simulation many times at various points in space (x,y), and then add
% together the resulting spectra.
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% tau1      = echo time 1 in [ms].
% tau2      = echo time 2 in [ms].
% RF        = RF pulse definition structure for refoc pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to first refocusing pulse) [cm]
% dy        = position offset in y-direction (corresponding to second refocusing pulse) [cm]
% Gx        = gradient strength for first selective refocusing pulse [G/cm]
% Gy        = gradient strength for second selective refocusing pulse [G/cm]
% flipAngle = flip angle of refocusing pulses [degrees] (Optional.  Default = 180 deg)
% centreFreq= centre frequency of the spectrum in [ppm] (Optional.  Default = 2.3)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using PRESS
%             sequence.

function out = sim_press_shaped(n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dx,dy,Gx,Gy,flipAngle,centreFreq)

if nargin<15
    centreFreq=2.3;
    if nargin<14
        flipAngle=180;
    end
end

if tau1<tp/1000
    error('ERROR:  Echo-time 1 cannot be less than duration of refocusing pulse! ABORTING!!');
end
if tau2<tp/1000
    error('ERROR:  Echo-time 2 cannot be less than duration of refocusing pulse! ABORTING!!');
end

%Set water to centre
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from tau1 and tau2;
delays=zeros(2);
delays(1)=tau1-tp;
delays(2)=tau2-tp;
if sum(delays<0)
    error(['ERROR! The following taus are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                                    %EXCITE
d=sim_COF(H,d,-1);                       %Keep only -1 coherences
d=sim_evolve(d,H,delays(1)/2000);                         %Evolve by delays(1)/2
d=sim_shapedRF(d,H,RF,tp,flipAngle,90,dx,Gx);             %1st shaped 180 degree refocusing pulse
d=sim_COF(H,d,1);                        %Keep only +1 coherences
d=sim_evolve(d,H,(delays(1)+delays(2))/2000);             %Evolve by (delays(1)+delays(2))/2
d=sim_shapedRF(d,H,RF,tp,flipAngle,90,dy,Gy);             %2nd shaped 180 degree refocusing pulse
d=sim_COF(H,d,-1);                       %Keep only -1 coherences
d=sim_evolve(d,H,delays(2)/2000);                         %Evolve by delays(2)/2
[out,dout]=sim_readout(d,H,n,sw,linewidth,90);            %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='press';
out.te=tau1+tau2;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

author: CIC-methods

repo: FID-A

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


text:
function rhs = globNedRHSVIFE3D(fun, mesh, fem, femI, dind, vind)

%% USAGE: generate global load vector on a tetrahedral mesh
%
% INPUTS:
% fun --- the load function from PDE (pde.f)
% mesh --- a struct data contains mesh information.
% fem --- global DoF for test function space
% dind --- derivative info for test function
%            d = [0,0,0]: function value
%            d = [1,0,0]: Dx value
%            d = [0,1,0]: Dy value
%            d = [0,0,1]: Dz value
%
% OUTPUTS:
% rhs --- global rhs vector

% Last Modified: 08/07/2020 by Xu Zhang

%% 1. RHS on non-interface elements
dof = size(fem.g2ldof,2);  nloc = dof;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);
AN = fem.area(ntID); gw = fem.gw;
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:);
X = zeros(nloc*ntN, 1);

feEvalBas = @EvalNed1Bas3D;
IfeEvalBas = @evalNed1IFEBas3D;

fN = feval(fun,gxN,gyN,gzN);
ind = 0;
I = reshape(femI.g2ldof(ntID,1:6),nloc*ntN,1);
for i = 1:dof
    ibas = feEvalBas(fem.bas, ntID, gxN, gyN, gzN, i, dind, vind);
    X(ind+1:ind+ntN) = AN.*sum((ibas.*fN).*gw',2).*fem.t_e_orit(ntID,i);
    ind = ind + ntN;
end
rhsN = sparse(I,1,X,length(femI.gdof),1);

%% 2. RHS on interface elements
AI1 = femI.area1; gw = femI.gw; gxI1 = femI.gx1; gyI1 = femI.gy1; gzI1 = femI.gz1;
ntI1 = size(femI.tType1,1); % not number of interface element, but quadrature element
nloc1 = 12; X1 = zeros(nloc1*ntI1, 1);
fI1 = feval(fun,gxI1,gyI1,gzI1);
ind = 0;
I1 = reshape(femI.tType1,nloc1*ntI1,1);
for i = 1:nloc1
    ibas = femI.basUType1(:,vind,i);
    X1(ind+1:ind+ntI1) = AI1.*sum((ibas.*fI1).*gw',2);
    ind = ind + ntI1;
end
rhsI1 = sparse(I1,1,X1,size(femI.gdof,1),1);

AI2 = femI.area2; gw = femI.gw; gxI2 = femI.gx2; gyI2 = femI.gy2; gzI2 = femI.gz2;
ntI2 = size(femI.tType2,1); % not number of interface element, but quadrature element
nloc2 = 14; X2 = zeros(nloc2*ntI2, 1);
fI2 = feval(fun,gxI2,gyI2,gzI2);
ind = 0;
I2 = reshape(femI.tType2,nloc2*ntI2,1);
for i = 1:nloc2
    ibas = femI.basUType2(:,vind,i);
    X2(ind+1:ind+ntI2) = AI2.*sum((ibas.*fI2).*gw',2);
    ind = ind + ntI2;
end
rhsI2 = sparse(I2,1,X2,size(femI.gdof,1),1);

rhs = rhsN + rhsI1 + rhsI2;

repo: ifem


text:
function DIC_3Djoined_results = joinPairs_old(DIC3DAllPairsResults)
%% function for joining multiple 3D-DIC pair results into one structure (plotting in step3)
%
% Inputs:
% * DIC3DAllPairsResults
%
% Outputs:
% * DIC_3Djoined_results
%

%%
nPairs=numel(DIC3DAllPairsResults);
nFrames=size(DIC3DAllPairsResults{1}.Points3D,1);

% create cells with all the variables for joining
Fi=cell(nPairs,1);
FCi=cell(nPairs,1);
Fisoi=cell(nPairs,1);
Vi=cell(nPairs,nFrames);
corrCombi=cell(nPairs,nFrames);
FaceCorrCombi=cell(nPairs,nFrames);
DispVi=cell(nPairs,nFrames);
DispMi=cell(nPairs,nFrames);
VFci=cell(nPairs,nFrames);
VFcTi=cell(nPairs,nFrames);
Ji=cell(nPairs,nFrames);
Lamda1i=cell(nPairs,nFrames);
Lamda2i=cell(nPairs,nFrames);
Emgni=cell(nPairs,nFrames);
emgni=cell(nPairs,nFrames);
Epc1i=cell(nPairs,nFrames);
Epc2i=cell(nPairs,nFrames);
Epc1veci=cell(nPairs,nFrames);
Epc1vecCuri=cell(nPairs,nFrames);
Epc2veci=cell(nPairs,nFrames);
Epc2vecCuri=cell(nPairs,nFrames);
epc1i=cell(nPairs,nFrames);
epc2i=cell(nPairs,nFrames);
epc1veci=cell(nPairs,nFrames);
epc2veci=cell(nPairs,nFrames);

for ip=1:nPairs
    DIC_3Dpair_struct=DIC3DAllPairsResults{ip};

    Fi{ip}=DIC_3Dpair_struct.Faces;
    FCi{ip}=DIC_3Dpair_struct.FaceColors;

    for it=1:nFrames
        Vi{ip,it}=DIC_3Dpair_struct.Points3D{it};
        corrCombi{ip,it}=DIC_3Dpair_struct.corrComb{it};
        FaceCorrCombi{ip,it}=DIC_3Dpair_struct.FaceCorrComb{it};
        DispVi{ip,it}=DIC_3Dpair_struct.Disp.DispVec{it};
        DispMi{ip,it}=DIC_3Dpair_struct.Disp.DispMgn{it};
        VFci{ip,it}=DIC_3Dpair_struct.FaceCentroids{it};
        Ji{ip,it}=DIC_3Dpair_struct.Deform.J{it};
        Lamda1i{ip,it}=DIC_3Dpair_struct.Deform.Lamda1{it};
        Lamda2i{ip,it}=DIC_3Dpair_struct.Deform.Lamda2{it};
        Emgni{ip,it}=DIC_3Dpair_struct.Deform.Emgn{it};
        emgni{ip,it}=DIC_3Dpair_struct.Deform.emgn{it};
        Epc1i{ip,it}=DIC_3Dpair_struct.Deform.Epc1{it};
        Epc2i{ip,it}=DIC_3Dpair_struct.Deform.Epc2{it};
        Epc1veci{ip,it}=DIC_3Dpair_struct.Deform.Epc1vec{it};
        Epc1vecCuri{ip,it}=DIC_3Dpair_struct.Deform.Epc1vecCur{it};
        Epc2veci{ip,it}=DIC_3Dpair_struct.Deform.Epc2vec{it};
        Epc2vecCuri{ip,it}=DIC_3Dpair_struct.Deform.Epc2vecCur{it};
        epc1i{ip,it}=DIC_3Dpair_struct.Deform.epc1{it};
        epc2i{ip,it}=DIC_3Dpair_struct.Deform.epc2{it};
        epc1veci{ip,it}=DIC_3Dpair_struct.Deform.epc1vec{it};
        epc2veci{ip,it}=DIC_3Dpair_struct.Deform.epc2vec{it};
        Fisoi{ip,it}=DIC_3Dpair_struct.FaceIsoInd{it};
    end

end

% join variables that are time invariant
[Fj,~,~]=joinElementSets(Fi,Vi(:,1),[]);
[~,~,FCj]=joinElementSets(Fi,Vi(:,1),FCi);

for it=1:nFrames
    % join variables that are time variant and scalar

    [~,Vj{1,it},DispMj{1,it}]=joinElementSets(Fi,Vi(:,it),DispMi(:,it));
    [~,~,corrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),corrCombi(:,it));
    [~,~,FaceCorrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),FaceCorrCombi(:,it));
    [~,~,Jj{1,it}]=joinElementSets(Fi,Vi(:,it),Ji(:,it));
    [~,~,Lamda1j{1,it}]=joinElementSets(Fi,Vi(:,it),Lamda1i(:,it));
    [~,~,Lamda2j{1,it}]=joinElementSets(Fi,Vi(:,it),Lamda2i(:,it));
    [~,~,Emgnj{1,it}]=joinElementSets(Fi,Vi(:,it),Emgni(:,it));
    [~,~,emgnj{1,it}]=joinElementSets(Fi,Vi(:,it),emgni(:,it));
    [~,~,Epc1j{1,it}]=joinElementSets(Fi,Vi(:,it),Epc1i(:,it));
    [~,~,Epc2j{1,it}]=joinElementSets(Fi,Vi(:,it),Epc2i(:,it));
    [~,~,epc1j{1,it}]=joinElementSets(Fi,Vi(:,it),epc1i(:,it));
    [~,~,epc2j{1,it}]=joinElementSets(Fi,Vi(:,it),epc2i(:,it));
    [~,~,Fisoj{1,it}]=joinElementSets(Fi,Vi(:,it),Fisoi(:,it));

    % join variables that are time variant and vectors
    DispVj{1,it}=DispVi{1,it};
    VFcj{1,it}=VFci{1,it};
    Epc1vecj{1,it}=Epc1veci{1,it};
    Epc1vecCurj{1,it}=Epc1vecCuri{1,it};
    Epc2vecj{1,it}=Epc2veci{1,it};
    Epc2vecCurj{1,it}=Epc2vecCuri{1,it};
    epc1vecj{1,it}=epc1veci{1,it};
    epc2vecj{1,it}=epc2veci{1,it};
    for ip=2:nPairs
        DispVj{1,it}=[DispVi{1,it}; DispVi{ip,it}];
        VFcj{1,it}=[VFcj{1,it}; VFci{ip,it}];
        Epc1vecj{1,it}=[Epc1vecj{1,it}; Epc1veci{ip,it}];
        Epc1vecCurj{1,it}=[Epc1vecCurj{1,it}; Epc1vecCuri{ip,it}];
        Epc2vecj{1,it}=[Epc2vecj{1,it}; Epc2veci{ip,it}];
        Epc2vecCurj{1,it}=[Epc2vecCurj{1,it}; Epc2vecCuri{ip,it}];
        epc1vecj{1,it}=[epc1vecj{1,it}; epc1veci{ip,it}];
        epc2vecj{1,it}=[epc2vecj{1,it}; epc2veci{ip,it}];
    end
end

DIC_3Djoined_results=struct;
DIC_3Djoined_results.Faces=Fj;
DIC_3Djoined_results.FaceColors=FCj;
DIC_3Djoined_results.FaceIsoInd=Fisoj;
DIC_3Djoined_results.FaceCentroids=VFcj;
DIC_3Djoined_results.Points3D=Vj;
DIC_3Djoined_results.corrComb=corrCombj;
DIC_3Djoined_results.FaceCorrComb=FaceCorrCombj;

DIC_3Djoined_results.Disp.DispVec=DispVj;
DIC_3Djoined_results.Disp.DispMgn=DispMj;

DIC_3Djoined_results.Deform.J=Jj;
DIC_3Djoined_results.Deform.Lamda1=Lamda1j;
DIC_3Djoined_results.Deform.Lamda2=Lamda2j;
DIC_3Djoined_results.Deform.Emgn=Emgnj;
DIC_3Djoined_results.Deform.emgn=emgnj;
DIC_3Djoined_results.Deform.Epc1=Epc1j;
DIC_3Djoined_results.Deform.Epc2=Epc2j;
DIC_3Djoined_results.Deform.Epc1vec=Epc1vecj;
DIC_3Djoined_results.Deform.Epc1vecCur=Epc1vecCurj;
DIC_3Djoined_results.Deform.Epc2vec=Epc2vecj;
DIC_3Djoined_results.Deform.Epc2vecCur=Epc2vecCurj;
DIC_3Djoined_results.Deform.epc1=epc1j;
DIC_3Djoined_results.Deform.epc2=epc2j;
DIC_3Djoined_results.Deform.epc1vec=epc1vecj;
DIC_3Djoined_results.Deform.epc2vec=epc2vecj;

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
%% adjacentdircount
% Below is a demonstration of the features of the |adjacentdircount| function

%% Syntax
% |[Lc]=adjacentdircount(L,n);|

%% Description
% The |adjacentdircount| function calculates the number of adjacent entries
% of an input logic array allong a specified direction.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
fontSize=10;

%% Example: Using |adjacentdircount| to compute "connectivity"
% The |adjacentdircount| function computes the number of adjacent elements
% in the input array that have the value 1. As such it computes a type of
% connectivity or "thickness" of structures in the logic arrays allong a
% certain dimension.

%%
% Lets consider the following array:

L=[0 1 1 0 0 1;...
   0 0 1 0 0 1;...
   1 0 0 1 1 1]; %A logic array or an array containing ones

%%
% Now view the output of the |adjacentdircount| function for n=1:
n=1; %Row direction (i.e. thickness orthogonal to rows)
[Lc]=adjacentdircount(L,n);
Lc

%%
% and similarly for for n=2:
n=2; %Column direction (i.e. thickness orthogonal to columns)
[Lc]=adjacentdircount(L,n);
Lc

%% Example: Use for 3D arrays
% Create a 3D example array, here L is copied upside down into a second
% "slice" and again for the third slice.
L(:,:,2)=flipud(L);
L(:,:,3)=flipud(L(:,:,1));

%%
%
n=3; %Slice direction (i.e. thickness orthogonal to slices)
[Lc]=adjacentdircount(L,n);
Lc

%% Example: Studying object line feature thickness

% Simulating a striped image in 3D
np=45;
[X,Y,Z]=meshgrid(linspace(-5*pi,5*pi,np));
R=sqrt(X.^2+Y.^2+Z.^2); %Radii
M=sin(R);

%%
% First consider the mid-slice for 2D analysis
midInd=round(size(M,3)/2);
m=M(:,:,midInd);

%%
% Convert to a logic array for thickness analysis
L=m>0;

%%
% Now use |adjacentdircount| to compute "connectivity" allong first and
% second dimensions

[Lc1]=adjacentdircount(L,1); %1st dim
[Lc2]=adjacentdircount(L,2); %1st dim

%%
% Plotting the results
hf=cFigure;
hold on;
subplot(1,3,1);
imagesc(L);
colormap gjet; colorbar;
axis equal; axis tight;
set(gca,'FontSize',fontSize);

subplot(1,3,2);
imagesc(Lc1);
colormap gjet; colorbar;
axis equal; axis tight;
set(gca,'FontSize',fontSize);

subplot(1,3,3);
imagesc(Lc2);
colormap gjet; colorbar;
axis equal; axis tight;
set(gca,'FontSize',fontSize);

drawnow;

%%
% Now lets expand to 3D analysis
L=M>0;
[Lc1]=adjacentdircount(L,1); %1st dim
[Lc2]=adjacentdircount(L,2); %1st dim
[Lc3]=adjacentdircount(L,3); %1st dim

%Construct logic array setting voxels to display using |ind2patch| function
logicPatch=L; %Initialise as L
logicPatch(:,:,midInd:end)=0; %Crop off half to visualize interior

[F,V,C]=ind2patch(logicPatch,L,'vb');
[F1,V1,C1]=ind2patch(logicPatch,Lc1,'vb');
[F2,V2,C2]=ind2patch(logicPatch,Lc2,'vb');
[F3,V3,C3]=ind2patch(logicPatch,Lc3,'vb');

cFigure;
subplot(2,2,1);
xlabel('J');ylabel('I'); zlabel('K'); hold on;
patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',C,'EdgeColor','k');
axis equal; view(3); axis tight; axis vis3d; grid off;
colormap gjet; colorbar;
set(gca,'FontSize',fontSize);

subplot(2,2,2);
xlabel('J');ylabel('I'); zlabel('K'); hold on;
patch('Faces',F1,'Vertices',V1,'FaceColor','flat','CData',C1,'EdgeColor','k');
axis equal; view(3); axis tight; axis vis3d; grid off;
colormap gjet; colorbar;
set(gca,'FontSize',fontSize);

subplot(2,2,3);
xlabel('J');ylabel('I'); zlabel('K'); hold on;
patch('Faces',F2,'Vertices',V2,'FaceColor','flat','CData',C2,'EdgeColor','k');
axis equal; view(3); axis tight; axis vis3d; grid off;
colormap gjet; colorbar;
set(gca,'FontSize',fontSize);

subplot(2,2,4);
xlabel('J');ylabel('I'); zlabel('K'); hold on;
patch('Faces',F3,'Vertices',V3,'FaceColor','flat','CData',C3,'EdgeColor','k');
axis equal; view(3); axis tight; axis vis3d; grid off;
colormap gjet; colorbar;
set(gca,'FontSize',fontSize);

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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


text:
function [gm samc] = mcurvature_vec(x,y,z)
% Description: The function calculates mean curvature and
% Surface Avg Mean Curvature (SAMC) of a surface formed by x, y & z.
% The input are the coordinate matrices x, y & z. The
% matrices can be formed using meshgrid or similar functions.
% This code is a vectorized form of original code (mcurvature) posted
% on Matlab File Exchange.

% The mean curvature is calculated according to the formula:

% If x:U->R^3 is a regular patch, then the mean curvature is given by
%
%         H = (eG-2fF+gE)/(2(EG-F^2)),
%
% where E, F, and G are coefficients of the first fundamental form and
% e, f, and g are coefficients of the second fundamental form

% Reference: Gray, A. "The Gaussian and Mean Curvatures." §16.5 in
% Modern Differential Geometry of Curves and Surfaces with Mathematica,
% 2nd ed. Boca Raton, FL: CRC Press, pp. 373-380, 1997 (p. 377).

% Inspired by:
% Title:  	Mean Curvature
% Author: 	Ahmed Elnaggar
% Summary: 	Calculate the Mean curvature of a given surface (x,y,z).

gm = zeros(size(z));
[xu,xv]     =   gradient(x);
[xuu,xuv]   =   gradient(xu);
[xvu,xvv]   =   gradient(xv);

[yu,yv]     =   gradient(y);
[yuu,yuv]   =   gradient(yu);
[yvu,yvv]   =   gradient(yv);

[zu,zv]     =   gradient(z);
[zuu,zuv]   =   gradient(zu);
[zvu,zvv]   =   gradient(zv);

Xu(:,:,1) = xu;
Xu(:,:,2) = yu;
Xu(:,:,3) = zu;

Xv(:,:,1) = xv;
Xv(:,:,2) = yv;
Xv(:,:,3) = zv;

Xuu(:,:,1) = xuu;
Xuu(:,:,2) = yuu;
Xuu(:,:,3) = zuu;

Xuv(:,:,1) = xuv;
Xuv(:,:,2) = yuv;
Xuv(:,:,3) = zuv;

Xvv(:,:,1) = xvv;
Xvv(:,:,2) = yvv;
Xvv(:,:,3) = zvv;

E = dot(Xu,Xu,3);
F           =   dot(Xu,Xv,3);
G           =   dot(Xv,Xv,3);
m           =   cross(Xu,Xv,3);
temp(:,:,1) = sqrt(sum(m.*m,3));
temp(:,:,2) = temp(:,:,1);
temp(:,:,3) = temp(:,:,1);
n           =   m./temp;
L           =   dot(Xuu,n,3);
M           =   dot(Xuv,n,3);
N           =   dot(Xvv,n,3);
gm          =   ((E.*N)+(G.*L)-(2.*F.*M))./(2.*(E.*G - F.^2));

dim = size(z);
samc = 1/ (dim(1) * dim(2)) * sum (gm(:).^2);


text:
function [out,out_se,out_dist] = loadM

    a = dir('model*mat');
    names = str2mat(a.name);
    eval(['load ' names(1,:)])

    out.energy = zeros(1,size(M.energy,2));
    out.conEnergy = zeros(1,size(M.conEnergy,2));
    out.eff_fitness = zeros(1,size(M.eff_fitness,2));
    out.eff = zeros(1,size(M.eff,2));
    out.se_contrasts = zeros(1,size(M.se_contrasts,2));
    out.resample_eff_loss = zeros(1,size(M.resample_eff_loss,2));
    out.smoothing_eff_loss = zeros(1,size(M.smoothing_eff_loss,2));
    out.vif = zeros(1,size(M.vif,2));
    out.conVif = zeros(1,size(M.conVif,2));
    out.cond = zeros(1,size(M.cond,2));
    out.conCond = zeros(1,size(M.conCond,2));
    out.conColin = zeros(1,size(M.conColin,2));
    out.hrf_eff_fitness = zeros(1,size(M.hrf_eff_fitness,2));
    out.hrf_eff = zeros(1,size(M.hrf_eff,2));
    out.hrf_eff_avg = zeros(1,size(M.hrf_eff_avg,2));
    out.cBal = zeros(1,size(M.cBal,2));

    out_se = out;
    out_dist = out;

    N = fieldnames(out);

    for j = N'
        eval([j{1} ' = [];']);
    end

    n = size(names,1);

warning off

for i = 1:n

    eval(['load ' names(i,:)])
    Mr = M;

    Mr = modelDiagnostics2(Mr,'noplot');

    for j = N'
        eval(['out.' j{1} ' = out.' j{1} ' + Mr.' j{1} ';'])
        eval(['out_dist.' j{1} ' = [out_dist.' j{1} '; (Mr.' j{1} ')];'])
    end

    if mod(i,100) == 0, fprintf(1,'.'), end

end

warning on

for j = N'
    eval(['out.' j{1} ' = out.' j{1} ' / n;'])
    eval(['out_se.' j{1} ' = std(out_dist.' j{1} ') ./ sqrt(n);'])
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore