text:
function x = mrdivide (b,F,is_inverse)
%MRDIVIDE x = b/A using the factorization F = factorize(A)
%
% Example
%   F = factorize(A) ;
%   x = b/F ;               % same as x=b/A
%
% See also factorize.

% Copyright 2009, Timothy A. Davis, University of Florida

if (nargin < 3)
    is_inverse = F.is_inverse ;
end

if (is_inverse)

    % x=b/inverse(A) is a double inverse, so it becomes simply x=b*A
    x = b*F.A ;

else

    bT = b' ;
    kind = F.kind ;
    switch kind

        case 1

            % minimum 2-norm solution of a sparse underdetermined problem
            % Q-less econonmy sparse QR factorization: (A*q)'*(A*q) = R'*R
            A = F.A ;
            R = F.R ;
            q = F.q ;
            x = A * (q * (R \ (R' \ (q' * bT)))) ;
            e = A * (q * (R \ (R' \ (q' * (bT - A' * x))))) ;
            x = (x + e)' ;

        case 2

            % minimum 2-norm solution of a dense underdetermined problem
            % dense economy QR factorization: A = Q*R
            % x = (Q * (R' \ b'))' ;
            Q = F.Q ;
            R = F.R ;
            if (issparse (bT))
                bT = full (bT) ;
            end
            opUT.UT = true ;
            opUT.TRANSA = true ;
            x = (Q * linsolve (R, bT, opUT))' ;

        case 3

            % least-squares solution of a sparse overdetermined problem
            % Q-less economy sparse QR factorization: (p*A)*(p*A)' = R'*R
            A = F.A ;
            R = F.R ;
            p = F.p ;
            x = p' * (R \ (R' \ (p * (A * bT)))) ;
            e = p' * (R \ (R' \ (p * (A * (bT - A' * x))))) ;
            x = (x + e)' ;

        case 4

            % least-squares solution of a dense overdetermined problem
            % dense economy QR factorization: A' = Q*R
            % x = (R \ (Q' * b'))'
            Q = F.Q ;
            R = F.R ;
            if (issparse (bT))
                bT = full (bT) ;
            end
            opU.UT = true ;
            x = linsolve (R, Q' * bT, opU)' ;

        case 5

            % sparse Cholesky factorization: q*A*q' = L*L'
            L = F.L ;
            q = F.q ;
            x = (q * (L' \ (L \ (q' * bT))))' ;

        case 6

            % dense Cholesky factorization: A = R'*R
            % x = (R \ (R' \ b'))'
            R = F.R ;
            if (issparse (bT))
                bT = full (bT) ;
            end
            opU.UT = true ;
            opUT.UT = true ;
            opUT.TRANSA = true ;
            x = linsolve (R, linsolve (R, bT, opUT), opU)' ;

        case 7

            % sparse LU factorization: p*A*q = L*U
            L = F.L ;
            U = F.U ;
            p = F.p ;
            q = F.q ;
            x = (p' * (L' \ (U' \ (q' * bT))))' ;

        case 8

            % dense LU factorization: p*A = L*U
            % x = (P' * (L' \ (U' \ b')))'
            L = F.L ;
            U = F.U ;
            p = F.p ;
            if (issparse (bT))
                bT = full (bT) ;
            end
            opUT.UT = true ;
            opUT.TRANSA = true ;
            opLT.LT = true ;
            opLT.TRANSA = true ;
            x = (p' * linsolve (L, linsolve (U, bT, opUT), opLT))' ;

    end
end

author: fieldtrip

repo: fieldtrip


text:
function varargout = process_mni_normalize( varargin )
% PROCESS_MNI_NORMALIZE: Compute deformation fields to the MNI ICBM152 space.
%
% USAGE:     OutputFiles = process_mni_normalize('Run', sProcess, sInputs)
%         [isOk, errMsg] = process_mni_normalize('Compute',            T1File, Method)
%                          process_mni_normalize('ComputeInteractive', T1File, Method=[ask], isUnload=1)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription()
    % Description the process
    sProcess.Comment     = 'MNI normalization';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 10;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ImportAnatomy#MNI_transformation';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    sProcess.isSeparator = 1;
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'NewSubject';
    % Method
    sProcess.options.method.Comment = {...
         'maff8:
Affine registration using SPM mutual information algorithm.
Estimates a simple 4x4 linear transformation to the MNI space.
Included in Brainstorm.', ...
         'segment:
Non-linear normalization and tissue classification with SPM12.
Estimates forward and inverse deformation fields to MNI space (IXI549).
Requires the installation of SPM12.'; ...
         'maff8', 'segment'};
    sProcess.options.method.Type    = 'radio_label';
    sProcess.options.method.Value   = 'maff8';
    % Use T2 when available
    sProcess.options.uset2.Comment = 'Use T2 when available ("segment" only)';
    sProcess.options.uset2.Type    = 'checkbox';
    sProcess.options.uset2.Value   = 0;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs)
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end
    % Get method
    Method = sProcess.options.method.Value;
    % Use T2
    UseT2 = sProcess.options.uset2.Value;

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy) || isempty(sSubject.iAnatomy)
        bst_report('Error', sProcess, [], ['No MRI available for subject "' SubjectName '".']);
        return
    end
    % Get default MRI as T1
    T1File = sSubject.Anatomy(sSubject.iAnatomy).FileName;
    % Look for T2 MRI
    T2File = [];
    if UseT2 && (length(sSubject.Anatomy) > 1)
        iT2 = find(~cellfun(@(c)isempty(strfind(c,'t2')), lower({sSubject.Anatomy.Comment})));
        % Warning when multiple T2 images
        if (length(iT2) > 1)
            bst_report('Warning', sProcess, [], ['Subject "' sSubject.Name '" has multiple anatomy volumes with tag "T2": not using T2 to avoid confusion.' 10 ...
                'To use a T2 MRI for better volume segmentation, rename or delete the extra T2-labelled files from the subject anatomy.']);
            iT2 = [];
        end
        % Get T2 filename
        if ~isempty(iT2)
            T2File = sSubject.Anatomy(iT2).FileName;
        end
    end

    % ===== COMPUTE MNI TRANSFORMATION =====
    % Call normalize function
    [sMriT1, errMsg] = Compute(T1File, Method, T2File);
    % Error handling
    if ~isempty(errMsg)
        bst_report('Error', sProcess, [], errMsg);
    end

    OutputFiles = {'import'};
end

%% ===== COMPUTE =====
function [sMriT1, errMsg] = Compute(T1File, Method, T2File)
    % Parse inputs
    if (nargin < 3) || isempty(T2File)
        T2File = [];
    end
    % Compute normalization
    [sMriT1, errMsg] = bst_normalize_mni(T1File, Method, T2File);
end

%% ===== COMPUTE/INTERACTIVE =====
function sMriT1 = ComputeInteractive(T1File, Method, isUnload, T2File)
    % Parse inputs
    if (nargin < 4) || isempty(T2File)
        T2File = [];
    end
    if (nargin < 3) || isempty(isUnload)
        isUnload = 1;
    end
    if (nargin < 2) || isempty(Method)
        Method = [];
    end
    % Unloading everything
    if isUnload
        bst_memory('UnloadAll', 'Forced');
    end
    % Ask method
    if isempty(Method)
        sProcess = GetDescription();
        Method = java_dialog('question', ['' sprintf('%s

', sProcess.options.method.Comment{1,:})], ...
            'MNI normalization method', [], {sProcess.options.method.Comment{2,:}, 'Cancel'}, sProcess.options.method.Comment{2,1});
        % Cancel
        if isempty(Method) || strcmpi(Method, 'Cancel')
            sMriT1 = [];
            return;
        end
    end
    % For Segment method: look for T2 MRI in the subject
    if strcmpi(Method, 'segment') && isempty(T2File)
        % Get subject info
        sSubject = bst_get('MriFile', T1File);
        % Find any possible T2
        iT2 = [];
        if (length(sSubject.Anatomy) > 1)
            iT2 = find(~cellfun(@(c)isempty(strfind(c,'t2')), lower({sSubject.Anatomy.Comment})));
            % Warning when multiple T2 images
            if (length(iT2) > 1)
                res = java_dialog('question', ...
                    ['Subject "' sSubject.Name '" has multiple anatomy volumes with tag "T2".' 10 10 ...
                     'If you want to use a T2 MRI for better volume segmentation, ' 10 ...
                     'rename or delete the extra T2-labelled files from the subject anatomy.' 10 10 ...
                     'If you click OK, no T2 image will be used in the segmentation.' 10 10], 'MNI normalization', [], {'OK', 'Cancel'}, 'OK');
                if isempty(res) || isequal(res, 'Cancel')
                    return;
                end
                iT2 = [];
            % Confirmation for using T2 image
            elseif (length(iT2) == 1)
                isT2 = java_dialog('confirm', ...
                    ['Subject "' sSubject.Name '" seems to include one T2 MRI, named "' sSubject.Anatomy(iT2).Comment '".' 10 10 ...
                     'Use this T2 image to refine the volume segmentation?' 10 10], 'MNI normalization');
                if ~isT2
                    iT2 = [];
                end
            end
        end
        % Get T2 filename
        if ~isempty(iT2)
            T2File = sSubject.Anatomy(iT2).FileName;
        end
    end
    % Open progress bar
    bst_progress('start', 'MNI normalization', 'Initialization...');
    % Call non-interactive function
    [sMriT1, errMsg] = Compute(T1File, Method, T2File);
    % Error handling
    if ~isempty(errMsg)
        sMriT1 = [];
        bst_error(errMsg, 'MNI normalization', 0);
    end
    % Close progress bar
    bst_progress('stop');
end

author: brainstorm-tools


text:
% finputcheck() - check Matlab function {'key','value'} input argument pairs
%
% Usage: >> result = finputcheck( varargin, fieldlist );
%        >> [result varargin] = finputcheck( varargin, fieldlist, ...
%                                              callingfunc, mode, verbose );
% Input:
%   varargin  - Cell array 'varargin' argument from a function call using 'key',
%               'value' argument pairs. See Matlab function 'varargin'
%   fieldlist - A 4-column cell array, one row per 'key'. The first
%               column contains the key string, the second its type(s),
%               the third the accepted value range, and the fourth the
%               default value.  Allowed types are 'boolean', 'integer',
%               'real', 'string', 'cell' or 'struct'.  For example,
%                       {'key1' 'string' { 'string1' 'string2' } 'defaultval_key1'}
%                       {'key2' {'real' 'integer'} { minint maxint } 'defaultval_key2'}
%  callingfunc - Calling function name for error messages. {default: none}.
%  mode        - ['ignore'|'error'] ignore keywords that are either not specified
%                in the fieldlist cell array or generate an error.
%                {default: 'error'}.
%  verbose     - ['verbose', 'quiet'] print information. Default: 'verbose'.
%
% Outputs:
%   result     - If no error, structure with 'key' as fields and 'value' as
%                content. If error this output contain the string error.
%   varargin   - residual varagin containing unrecognized input arguments.
%                Requires mode 'ignore' above.
%
% Note: In case of error, a string is returned containing the error message
%       instead of a structure.
%
% Example (insert the following at the beginning of your function):
%	result = finputcheck(varargin, ...
%               { 'title'         'string'   []       ''; ...
%                 'percent'       'real'     [0 1]    1 ; ...
%                 'elecamp'       'integer'  [1:10]   [] });
%   if isstr(result)
%       error(result);
%   end
%
% Note:
%   The 'title' argument should be a string. {no default value}
%   The 'percent' argument should be a real number between 0 and 1. {default: 1}
%   The 'elecamp' argument should be an integer between 1 and 10 (inclusive).
%
%   Now 'g.title' will contain the title arg (if any, else the default ''), etc.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 10 July 2002

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 10 July 2002, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [g, varargnew] = finputcheck( vararg, fieldlist, callfunc, mode, verbose )

	if nargin < 2
		help finputcheck;
		return;
	end;
	if nargin < 3
		callfunc = '';
	else
		callfunc = [callfunc ' ' ];
	end;
    if nargin < 4
        mode = 'do not ignore';
    end;
    if nargin < 5
        verbose = 'verbose';
    end;
	NAME = 1;
	TYPE = 2;
	VALS = 3;
	DEF  = 4;
	SIZE = 5;

	varargnew = {};
	% create structure
	% ----------------
	if ~isempty(vararg)
		for index=1:length(vararg)
			if iscell(vararg{index})
				vararg{index} = {vararg{index}};
			end;
		end;
		try
			g = struct(vararg{:});
		catch
            vararg = removedup(vararg, verbose);
            try
                g = struct(vararg{:});
            catch
                g = [ callfunc 'error: bad ''key'', ''val'' sequence' ]; return;
            end;
		end;
	else
		g = [];
	end;

	for index = 1:size(fieldlist,NAME)
		% check if present
		% ----------------
		if ~isfield(g, fieldlist{index, NAME})
			g = setfield( g, fieldlist{index, NAME}, fieldlist{index, DEF});
		end;
		tmpval = getfield( g, {1}, fieldlist{index, NAME});

		% check type
		% ----------
        if ~iscell( fieldlist{index, TYPE} )
            res = fieldtest( fieldlist{index, NAME},  fieldlist{index, TYPE}, ...
                           fieldlist{index, VALS}, tmpval, callfunc );
            if isstr(res), g = res; return; end;
        else
            testres = 0;
            tmplist = fieldlist;
            for it = 1:length( fieldlist{index, TYPE} )
                if ~iscell(fieldlist{index, VALS})
                     res{it} = fieldtest(  fieldlist{index, NAME},  fieldlist{index, TYPE}{it}, ...
                                           fieldlist{index, VALS}, tmpval, callfunc );
                else res{it} = fieldtest(  fieldlist{index, NAME},  fieldlist{index, TYPE}{it}, ...
                                           fieldlist{index, VALS}{it}, tmpval, callfunc );
                end;
                if ~isstr(res{it}), testres = 1; end;
            end;
            if testres == 0,
                g = res{1};
                for tmpi = 2:length(res)
                    g = [ g 10 'or ' res{tmpi} ];
                end;
                return;
            end;
        end;
	end;

    % check if fields are defined
	% ---------------------------
	allfields = fieldnames(g);
	for index=1:length(allfields)
		if isempty(strmatch(allfields{index}, fieldlist(:, 1)', 'exact'))
			if ~strcmpi(mode, 'ignore')
				g = [ callfunc 'error: undefined argument ''' allfields{index} '''']; return;
			end;
			varargnew{end+1} = allfields{index};
			varargnew{end+1} = getfield(g, {1}, allfields{index});
		end;
	end;

function g = fieldtest( fieldname, fieldtype, fieldval, tmpval, callfunc );
	NAME = 1;
	TYPE = 2;
	VALS = 3;
	DEF  = 4;
	SIZE = 5;
    g = [];

    switch fieldtype
     case { 'integer' 'real' 'boolean' 'float' },
      if ~isnumeric(tmpval) && ~islogical(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be numeric' ]; return;
      end;
      if strcmpi(fieldtype, 'boolean')
          if tmpval ~=0 && tmpval ~= 1
              g = [ callfunc 'error: argument ''' fieldname ''' must be 0 or 1' ]; return;
          end;
      else
          if strcmpi(fieldtype, 'integer')
              if ~isempty(fieldval)
                  if (any(isnan(tmpval(:))) && ~any(isnan(fieldval))) ...
                          && (~ismember(tmpval, fieldval))
                      g = [ callfunc 'error: wrong value for argument ''' fieldname '''' ]; return;
                  end;
              end;
          else % real or float
              if ~isempty(fieldval) && ~isempty(tmpval)
                  if any(tmpval < fieldval(1)) || any(tmpval > fieldval(2))
                      g = [ callfunc 'error: value out of range for argument ''' fieldname '''' ]; return;
                  end;
              end;
          end;
      end;

     case 'string'
      if ~isstr(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a string' ]; return;
      end;
      if ~isempty(fieldval)
          if isempty(strmatch(lower(tmpval), lower(fieldval), 'exact'))
              g = [ callfunc 'error: wrong value for argument ''' fieldname '''' ]; return;
          end;
      end;

     case 'cell'
      if ~iscell(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a cell array' ]; return;
      end;

     case 'struct'
      if ~isstruct(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a structure' ]; return;
      end;

     case '';
     otherwise, error([ 'finputcheck error: unrecognized type ''' fieldname '''' ]);
    end;

% remove duplicates in the list of parameters
% -------------------------------------------
function cella = removedup(cella, verbose)
% make sure if all the values passed to unique() are strings, if not, exist
%try
    [tmp indices] = unique(cella(1:2:end));
    if length(tmp) ~= length(cella)/2
        myfprintf(verbose,'Note: duplicate ''key'', ''val'' parameter(s), keeping the last one(s)\n');
    end;
    cella = cella(sort(union(indices*2-1, indices*2)));
%catch
    % some elements of cella were not string
%    error('some ''key'' values are not string.');
%end;

function myfprintf(verbose, varargin)

if strcmpi(verbose, 'verbose')
    fprintf(varargin{:});
end;

author: PatternRecognition

repo: OpenBMI


%==============================================================================
% This code is part of the Finite Element Method app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIRFEM
%==============================================================================
%
% classdef TetraMesh1 < handle
%
% Finite Element Mesh based on tetrahedral subdivision of rectangular mesh.
%
% Each Cell is divided into 24 tetrahedra:
%
%
%  To construct an instance of this class type:
%
%  >> Mesh = TriMesh1(omega,m)
%
% Input:
% 	omega - description of spatial domain
%   m     - number of cells
%
% Properties:
%   xn     - node list
%   tri    - triangle list
%   dim    - space dimension
%   omega  - description of spatial domain
%   m      - number of cells
%   type   - type of partition
%   vol    - volume of triangles
%   nnodes - number of nodes
%   ntri   - number of triangles
%   dx1    - partial derivative operator
%   dx2    - partial derivative operator
%   dx3    - partial derivative operator
%   GRAD   - gradient operator
%   P1     - projection operator for node 1
%   P2     - projection operator for node 2
%   P3     - projection operator for node 3
%   P4     - projection operator for node 4
%   PC     - projection operator for Barycentrum
%   P      - prolongation operator
%   Pt     - restriction operator
%
%  Methods:
%   mfPu   - matrix free prolongation/restriction
%   mfPi   - matrix free edge projector
%   getP   - builds prolongation operator
%   tri2cc - averaging
%
%
% see also
% =========================================================================
classdef TetraMesh1 < handle

    properties
        % ===================================
        % node list
        % ===================================
        xn
        % ===================================
        % triangle list
        % ===================================
        tri
        % ===================================
        % space dimension
        % ===================================
        dim  = 3;
        % ===================================
        % description of computational domain
        % ===================================
        omega
        % ===================================
        % number of cells
        % ===================================
        m
        % ===================================
        % type of partition
        % ===================================
        type = 1;
        % ===================================
        % function handle to myself
        % ===================================
        me = @TetraMesh1;
        % ===================================
        % volume of triangles
        % ===================================
        vol
        % ===================================
        % number of nodes
        % ===================================
        nnodes
        % ===================================
        % number of triangles in mesh
        % ===================================
        ntri
    end

    properties (Access = public, Dependent) % These will be created when first callend and stores persistently
        % ===================================
        % dx1  - partial derivative operator
        % ===================================
        dx1
        % ===================================
        % dx2  - Partial derivative operator
        % ===================================
        dx2
        % ===================================
        % dx3  - Partial derivative operator
        % ===================================
        dx3
        % ===================================
        % GRAD - Gradient operator
        %
        %         | dx1 |
        %         |     |
        %  GRAD = | dx2 |
        %         |     |
        %         | dx3 |
        %
        % ===================================
        GRAD
        % ===================================
        % B - Vector gradient operator
        %
        %         | GRAD    0   0   |
        %         |                 |
        %   B =   |   0   GRAD  0   |
        %         |                 |
        %         |   0     0  GRAD |
        %
        % ===================================
        B
        % ===================================
        % P1 - Projection operator on Node 1
        % ===================================
        P1
        % ===================================
        % P2 - Projection operator on Node 2
        % ===================================
        P2
        % ===================================
        % P3 - Projection operator on Node 3
        % ===================================
        P3
        % ===================================
        % P4 - Projection operator on Node 4
        % ===================================
        P4
        % ===================================
        % PC - Projection operator on Barycenter
        % ===================================
        PC
        % ===================================
        % P  - Prolongation operator
        % ===================================
        P
        % ===================================
        % Pt - Restriction operator
        % ===================================
        Pt
        % ===================================
        % Boundary indices
        % ===================================
        boundaryIdx
        % ===================================
        % Boundary projector
        % ===================================
        boundaryProj
        mfdx1
        mfdx2
        mfdx3
        mfGRAD
    end

    properties (Access = private)
        % These are where the dependent data is actually stored
        dx1_
        dx2_
        dx3_
        GRAD_
        B_
        P1_
        P2_
        P3_
        P4_
        PC_
        P_
        Pt_
        mfdx1_
        mfdx2_
        mfdx3_
        mfGRAD_
        boundaryIdx_
        boundaryProj_
    end

    methods
        function this = TetraMesh1(omega,m)
            if nargin==0,
                help(mfilename);
                this.runMinimalExample;
                return;
            end
            h          = (omega(2:2:end)-omega(1:2:end))./m;
            this.omega = omega;
            this.m     = m;

            % get nodes (combine nodal, stg-1,stg-2,stg-3,cc grids)
            xc = @(i) linspace(omega(2*i-1)+h(i)/2,omega(2*i)-h(i)/2,m(i))'; % cell centers
            xn = @(i) linspace(omega(2*i-1),omega(2*i),m(i)+1)'; % nodes
            % add nodal grid
            nn = reshape(1:prod(m+1),m+1);
            nn = reshape(nn(1:end-1,1:end-1,1:end-1),1,[]); % get indices of bottom left vertices
            this.xn = reshape(getNodalGrid(omega,m),[],3);
            % add stg-1 grid
            ns  = prod(ones(length(m),1)*m+eye(length(m)),2); %staggered
            ns1 = reshape(size(this.xn,1)+(1:ns(1)),[m(1)+1 m(2) m(3)]);
            ns1 = reshape(ns1(1:end-1,:,:),1,[]); % get indices of bottom left vertices
            [x1,x2,x3] = ndgrid(xn(1),xc(2),xc(3));
            this.xn = [this.xn; x1(:), x2(:), x3(:)];
            % add stg-2 grid
            ns2 = reshape(size(this.xn,1)+(1:ns(2)),[m(1) m(2)+1 m(3)]);
            ns2 = reshape(ns2(:,1:end-1,:),1,[]); % get indices of bottom left vertices
            [x1,x2,x3] = ndgrid(xc(1),xn(2),xc(3));
            this.xn = [this.xn; x1(:), x2(:), x3(:)];
            % add stg-3 grid
            ns3 = reshape(size(this.xn,1)+(1:ns(3)),[m(1) m(2) m(3)+1]);
            ns3 = reshape(ns3(:,:,1:end-1),1,[]); % get indices of bottom left vertices
            [x1,x2,x3] = ndgrid(xc(1),xc(2),xn(3));
            this.xn = [this.xn; x1(:), x2(:), x3(:)];
            % add cc grid
            nc = reshape(size(this.xn,1)+(1:prod(m)),m);
            nc = reshape(nc,1,[]); % get indices of bottom left vertices
            this.xn = [this.xn; reshape(getCellCenteredGrid(omega,m),[],3)];

            % specify triangles
            iyn = m(1)+1; izn = prod(m(1:2)+1); izc = prod(m(1:2));
            this.tri   = [
                nn;           nn+1;         ns2;      nc;
                nn+1;         nn+izn+1;     ns2;      nc;
                nn+izn+1;     nn+izn;       ns2;      nc;
                nn+izn;       nn;           ns2;      nc; % 4
                nn+1;         nn+iyn+1;     ns1+1;    nc;
                nn+iyn+1;     nn+iyn+1+izn; ns1+1;    nc;
                nn+iyn+1+izn; nn+1+izn;     ns1+1;    nc;
                nn+1+izn;     nn+1;         ns1+1;    nc; % 8
                nn+1;         nn;           ns3;      nc;
                nn;           nn+iyn;       ns3;      nc;
                nn+iyn;       nn+1+iyn;     ns3;      nc;
                nn+1+iyn;     nn+1;         ns3;      nc; % 12

                nn+iyn;       nn;           ns1;      nc;
                nn;           nn+izn;       ns1;      nc;
                nn+izn;       nn+iyn+izn;   ns1;      nc;
                nn+iyn+izn;   nn+iyn;       ns1;      nc; % 16

                nn+iyn+1;     nn+iyn;       ns2+m(1); nc;
                nn+iyn;       nn+iyn+izn;   ns2+m(1); nc;
                nn+iyn+izn;   nn+iyn+1+izn; ns2+m(1); nc;
                nn+iyn+1+izn; nn+iyn+1;     ns2+m(1); nc; % 20
                nn+izn;       nn+izn+1;     ns3+izc;  nc;
                nn+izn+1;     nn+izn+1+iyn; ns3+izc;  nc;
                nn+iyn+1+izn; nn+iyn+izn;   ns3+izc;  nc;
                nn+iyn+izn;   nn+izn;       ns3+izc;  nc; % 24
                ];
            this.tri   = reshape(this.tri,4,[])';
            this.nnodes = size(this.xn,1);
            this.ntri  = size(this.tri,1);
            this.vol = prod((omega(2:2:end)-omega(1:2:end))./m)/24*ones(this.ntri,1);
        end

        function runMinimalExample(~)
            omega = [0 4 2 6 0 3]; m = [3 4 6];
            Mesh  = feval(mfilename,omega,m);
        end

        function x = mfPi(this,x,i)
            % =============================================================
            % function x = mfPi(this,x,i)
            %
            % matrix free edge projector
            % =============================================================
            switch i
                case 1
                    P = this.P1;
                case 2
                    P = this.P2;
                case 3
                    P = this.P3;
                case 4
                    P = this.P4;
                case 'C'
                    P = this.PC;
            end
            if size(x,1) == this.ntri,
                % ajoint
                x = P'*x;
            else
                x = P * x;
            end
        end

        function x = tri2cc(this,x)
            % =============================================================
            % function x = tri2cc(x)
            %
            % averaging or adjoint
            % =============================================================
            if numel(x)==this.ntri,
                x = mean(reshape(x,24,[]),1);
            else
                x = reshape(x,1,[]);
                x = (1/24)*repmat(x,[24 1]);
                x = x(:);
            end

        end

        function Pu = mfPuNodal(this,yn,m,flag)
            % =============================================================
            % function Pu = mfPuNodal(~,yn,m,flag)
            %
            % matrix free prolongation/restriction for nodal quantities
            % =============================================================
            yn = reshape(yn,[],this.dim);
            v = @(x) x(:);
            switch flag
                case 'Pu' % coarse --> fine
                    mf = 2*m;

                    % indices of coarse grid nodes (split by nodal, stg-i, c)
                    ns   = prod(ones(length(m),1)*m+eye(length(m)),2); %staggered
                    inc  = reshape(1:prod(m+1),m+1);
                    is1c = reshape(prod(m+1)+(1:ns(1)),[m(1)+1 m(2) m(3)]);
                    is2c = reshape(prod(m+1)+ns(1)+ (1:ns(2)),[m(1) m(2)+1 m(3)]);
                    is3c = reshape(prod(m+1)+ns(1)+ns(2)+(1:ns(3)),[m(1) m(2) m(3)+1]);
                    icc  = reshape(prod(m+1)+ns(1)+ns(2)+ns(3)+(1:prod(m)),m);
                    % indices of fine grid nodes
                    ns   = prod(ones(length(mf),1)*mf+eye(length(mf)),2); %staggered
                    inf  = reshape(1:prod(mf+1),mf+1);
                    is1f = reshape(prod(mf+1)+(1:ns(1)),[mf(1)+1 mf(2) mf(3)]);
                    is2f = reshape(prod(mf+1)+ns(1)+ (1:ns(2)),[mf(1) mf(2)+1 mf(3)]);
                    is3f = reshape(prod(mf+1)+ns(1)+ns(2)+(1:ns(3)),[mf(1) mf(2) mf(3)+1]);
                    icf  = reshape(prod(mf+1)+ns(1)+ns(2)+ns(3)+(1:prod(mf)),mf);

                    % allocate space
                    Pu = zeros(icf(end),this.dim);

                    % include existing nodes  (nodal-->nodal)
                    Pu(v(inf(1:2:end,1:2:end,1:2:end)),:)  = yn(inc(:),:);
                    % include existing nodes  (stg-1-->nodal)
                    Pu(v(inf(1:2:end,2:2:end,2:2:end)),:) = yn(is1c(:),:);
                    % include existing nodes  (stg-2-->nodal)
                    Pu(v(inf(2:2:end,1:2:end,2:2:end)),:) = yn(is2c(:),:);
                    % include existing nodes  (stg-3-->nodal)
                    Pu(v(inf(2:2:end,2:2:end,1:2:end)),:) = yn(is3c(:),:);
                    % include existing nodes  (cc-->nodal)
                    Pu(v(inf(2:2:end,2:2:end,2:2:end)),:) = yn(icc(:),:);

                    % average to get edge-stg-1
                    Pu(v(inf(2:2:end,1:2:end,1:2:end)),:) = ...
                        .5*(yn(v(inc(1:end-1,:,:)),:) + yn(v(inc(2:end,:,:)),:)) ;
                    % average to get edge-stg-2
                    Pu(v(inf(1:2:end,2:2:end,1:2:end)),:) = ...
                        .5*(yn(v(inc(:,1:end-1,:)),:) + yn(v(inc(:,2:end,:)),:)) ;
                    % average to get edge-stg-3
                    Pu(v(inf(1:2:end,1:2:end,2:2:end)),:) = ...
                        .5*(yn(v(inc(:,:,1:end-1)),:) + yn(v(inc(:,:,2:end)),:)) ;
                    % get face-stg-1
                    Pu(is1f(:),:) = .25*(   Pu(v(inf(1:end,1:end-1,1:end-1)),:) ...
                        + Pu(v(inf(1:end,2:end  ,1:end-1)),:) ...
                        + Pu(v(inf(1:end,2:end  ,2:end  )),:) ...
                        + Pu(v(inf(1:end,1:end-1,2:end  )),:));

                    % get face-stg-2
                    Pu(is2f(:),:) = .25*(   Pu(v(inf(1:end-1,1:end  ,1:end-1)),:) ...
                        + Pu(v(inf(2:end  ,1:end  ,1:end-1)),:) ...
                        + Pu(v(inf(2:end  ,1:end  ,2:end  )),:) ...
                        + Pu(v(inf(1:end-1,1:end  ,2:end  )),:));
                    % get face-stg-3
                    Pu(is3f(:),:) = .25*(   Pu(v(inf(1:end-1,1:end-1,1:end)),:) ...
                        + Pu(v(inf(2:end  ,1:end-1,1:end)),:) ...
                        + Pu(v(inf(2:end  ,2:end  ,1:end)),:) ...
                        + Pu(v(inf(1:end-1,2:end  ,1:end)),:));

                    % get new cell-centers
                    Pu(icf(:),:)  = .125*(   Pu(v(inf(1:end-1,1:end-1,1:end-1)),:) ...
                        + Pu(v(inf(2:end  ,1:end-1,1:end-1)),:) ...
                        + Pu(v(inf(2:end  ,2:end  ,1:end-1)),:) ...
                        + Pu(v(inf(1:end-1,2:end  ,1:end-1)),:) ...
                        + Pu(v(inf(1:end-1,1:end-1,2:end)),:) ...
                        + Pu(v(inf(2:end  ,1:end-1,2:end)),:) ...
                        + Pu(v(inf(2:end  ,2:end  ,2:end)),:) ...
                        + Pu(v(inf(1:end-1,2:end  ,2:end)),:));

                case 'PTu' % fine --> coarse
                    % include parent nodes
                    mf = m;
                    m  = m/2;
                    % indices of coarse grid nodes (split by nodal, stg-i, c)
                    ns   = prod(ones(length(m),1)*m+eye(length(m)),2); %staggered
                    inc  = reshape(1:prod(m+1),m+1);
                    is1c = reshape(prod(m+1)+(1:ns(1)),[m(1)+1 m(2) m(3)]);
                    is2c = reshape(prod(m+1)+ns(1)+ (1:ns(2)),[m(1) m(2)+1 m(3)]);
                    is3c = reshape(prod(m+1)+ns(1)+ns(2)+(1:ns(3)),[m(1) m(2) m(3)+1]);
                    icc  = reshape(prod(m+1)+ns(1)+ns(2)+ns(3)+(1:prod(m)),m);
                    % indices of fine grid nodes
                    ns   = prod(ones(length(mf),1)*mf+eye(length(mf)),2); %staggered
                    inf  = reshape(1:prod(mf+1),mf+1);
                    is1f = reshape(prod(mf+1)+(1:ns(1)),[mf(1)+1 mf(2) mf(3)]);
                    is2f = reshape(prod(mf+1)+ns(1)+ (1:ns(2)),[mf(1) mf(2)+1 mf(3)]);
                    is3f = reshape(prod(mf+1)+ns(1)+ns(2)+(1:ns(3)),[mf(1) mf(2) mf(3)+1]);
                    icf  = reshape(prod(mf+1)+ns(1)+ns(2)+ns(3)+(1:prod(mf)),mf);

                    % allocate space
                    Pu = zeros(icc(end),this.dim);

                    % push weights from cell-centers
                    yn(v(inf(1:end-1,1:end-1,1:end-1)),:)  = yn(v(inf(1:end-1,1:end-1,1:end-1)),:) + .125* yn(icf(:),:);
                    yn(v(inf(2:end  ,1:end-1,1:end-1)),:)  = yn(v(inf(2:end  ,1:end-1,1:end-1)),:) + .125* yn(icf(:),:);
                    yn(v(inf(2:end  ,2:end  ,1:end-1)),:)  = yn(v(inf(2:end  ,2:end  ,1:end-1)),:) + .125* yn(icf(:),:);
                    yn(v(inf(1:end-1,2:end  ,1:end-1)),:)  = yn(v(inf(1:end-1,2:end  ,1:end-1)),:) + .125* yn(icf(:),:);
                    yn(v(inf(1:end-1,1:end-1,2:end)),:)    = yn(v(inf(1:end-1,1:end-1,2:end)),:)   + .125* yn(icf(:),:);
                    yn(v(inf(2:end  ,1:end-1,2:end)),:)    = yn(v(inf(2:end  ,1:end-1,2:end)),:)   + .125* yn(icf(:),:);
                    yn(v(inf(2:end  ,2:end  ,2:end)),:)    = yn(v(inf(2:end  ,2:end  ,2:end)),:)   + .125* yn(icf(:),:);
                    yn(v(inf(1:end-1,2:end  ,2:end)),:)    = yn(v(inf(1:end-1,2:end  ,2:end)),:)   + .125* yn(icf(:),:);
                    % push weights from face-stg-1
                    yn(v(inf(1:end,1:end-1,1:end-1)),:)   = yn(v(inf(1:end,1:end-1,1:end-1)),:) + .25 * yn(is1f(:),:);
                    yn(v(inf(1:end,2:end  ,1:end-1)),:)   = yn(v(inf(1:end,2:end  ,1:end-1)),:) + .25 * yn(is1f(:),:);
                    yn(v(inf(1:end,2:end  ,2:end  )),:)   = yn(v(inf(1:end,2:end  ,2:end  )),:) + .25 * yn(is1f(:),:);
                    yn(v(inf(1:end,1:end-1,2:end  )),:)   = yn(v(inf(1:end,1:end-1,2:end  )),:) + .25 * yn(is1f(:),:);
                    % push weights from face-stg-2
                    yn(v(inf(1:end-1,1:end  ,1:end-1)),:) = yn(v(inf(1:end-1,1:end  ,1:end-1)),:) + .25 * yn(is2f(:),:);
                    yn(v(inf(2:end  ,1:end  ,1:end-1)),:) = yn(v(inf(2:end  ,1:end  ,1:end-1)),:) + .25 * yn(is2f(:),:);
                    yn(v(inf(2:end  ,1:end  ,2:end  )),:) = yn(v(inf(2:end  ,1:end  ,2:end  )),:) + .25 * yn(is2f(:),:);
                    yn(v(inf(1:end-1,1:end  ,2:end  )),:) = yn(v(inf(1:end-1,1:end  ,2:end  )),:) + .25 * yn(is2f(:),:);
                    % push weights from face-stg-3
                    yn(v(inf(1:end-1,1:end-1,1:end)),:)   = yn(v(inf(1:end-1,1:end-1,1:end)),:)  + .25 *  yn(is3f(:),:);
                    yn(v(inf(2:end  ,1:end-1,1:end)),:)   = yn(v(inf(2:end  ,1:end-1,1:end)),:)  + .25 *  yn(is3f(:),:);
                    yn(v(inf(2:end  ,2:end  ,1:end)),:)   = yn(v(inf(2:end  ,2:end  ,1:end)),:)  + .25 *  yn(is3f(:),:);
                    yn(v(inf(1:end-1,2:end  ,1:end)),:)   = yn(v(inf(1:end-1,2:end  ,1:end)),:)  + .25 *  yn(is3f(:),:);

                    % include existing nodes  (nodal-->nodal)
                    Pu(inc(:),:)  = yn(v(inf(1:2:end,1:2:end,1:2:end)),:) ;
                    % include existing nodes  (stg-1-->nodal)
                    Pu(is1c(:),:) = yn(v(inf(1:2:end,2:2:end,2:2:end)),:);
                    % include existing nodes  (stg-2-->nodal)
                    Pu(is2c(:),:) = yn(v(inf(2:2:end,1:2:end,2:2:end)),:);
                    % include existing nodes  (stg-3-->nodal)
                    Pu(is3c(:),:) = yn(v(inf(2:2:end,2:2:end,1:2:end)),:);
                    % include existing nodes  (cc-->nodal)
                    Pu(icc(:),:)  = yn(v(inf(2:2:end,2:2:end,2:2:end)),:);

                    % average to get edge-stg-1
                    Pu(v(inc(1:end-1,:,:)),:) = Pu(v(inc(1:end-1,:,:)),:) + .5 * yn(v(inf(2:2:end,1:2:end,1:2:end)),:);
                    Pu(v(inc(2:end,:,:)),:)   = Pu(v(inc(2:end,:,:)),:)   + .5 * yn(v(inf(2:2:end,1:2:end,1:2:end)),:);

                    % average to get edge-stg-2
                    Pu(v(inc(:,1:end-1,:)),:) = Pu(v(inc(:,1:end-1,:)),:) + .5 * yn(v(inf(1:2:end,2:2:end,1:2:end)),:);
                    Pu(v(inc(:,2:end  ,:)),:) = Pu(v(inc(:,2:end  ,:)),:) + .5 * yn(v(inf(1:2:end,2:2:end,1:2:end)),:);

                    % average to get edge-stg-3
                    Pu(v(inc(:,:,1:end-1)),:) = Pu(v(inc(:,:,1:end-1)),:) + .5 * yn(v(inf(1:2:end,1:2:end,2:2:end)),:);
                    Pu(v(inc(:,:,2:end  )),:) = Pu(v(inc(:,:,2:end  )),:) + .5 * yn(v(inf(1:2:end,1:2:end,2:2:end)),:);

            end
            Pu = Pu(:);

        end
        function Pu = mfInterNodal(this,yn,m)
            % =============================================================
            % function Pu = mfInterNodal(~,yn,m,)
            %
            % interpolates to nodal grid of a finer resolution.
            % =============================================================
            yn = reshape(yn,[],this.dim);
            v = @(x) x(:);
            mf = 2*m;

            % indices of coarse grid nodes (split by nodal, stg-i, c)
            ns   = prod(ones(length(m),1)*m+eye(length(m)),2); %staggered
            inc  = reshape(1:prod(m+1),m+1);
            is1c = reshape(prod(m+1)+(1:ns(1)),[m(1)+1 m(2) m(3)]);
            is2c = reshape(prod(m+1)+ns(1)+ (1:ns(2)),[m(1) m(2)+1 m(3)]);
            is3c = reshape(prod(m+1)+ns(1)+ns(2)+(1:ns(3)),[m(1) m(2) m(3)+1]);
            icc  = reshape(prod(m+1)+ns(1)+ns(2)+ns(3)+(1:prod(m)),m);
            % indices of fine grid nodes
            inf  = reshape(1:prod(mf+1),mf+1);

            % allocate space
            Pu = zeros(inf(end),this.dim);

            % include existing nodes  (nodal-->nodal)
            Pu(v(inf(1:2:end,1:2:end,1:2:end)),:)  = yn(inc(:),:);
            % include existing nodes  (stg-1-->nodal)
            Pu(v(inf(1:2:end,2:2:end,2:2:end)),:) = yn(is1c(:),:);
            % include existing nodes  (stg-2-->nodal)
            Pu(v(inf(2:2:end,1:2:end,2:2:end)),:) = yn(is2c(:),:);
            % include existing nodes  (stg-3-->nodal)
            Pu(v(inf(2:2:end,2:2:end,1:2:end)),:) = yn(is3c(:),:);
            % include existing nodes  (cc-->nodal)
            Pu(v(inf(2:2:end,2:2:end,2:2:end)),:) = yn(icc(:),:);

            % average to get edge-stg-1
            Pu(v(inf(2:2:end,1:2:end,1:2:end)),:) = ...
                .5*(yn(v(inc(1:end-1,:,:)),:) + yn(v(inc(2:end,:,:)),:)) ;
            % average to get edge-stg-2
            Pu(v(inf(1:2:end,2:2:end,1:2:end)),:) = ...
                .5*(yn(v(inc(:,1:end-1,:)),:) + yn(v(inc(:,2:end,:)),:)) ;
            % average to get edge-stg-3
            Pu(v(inf(1:2:end,1:2:end,2:2:end)),:) = ...
                .5*(yn(v(inc(:,:,1:end-1)),:) + yn(v(inc(:,:,2:end)),:)) ;

            Pu = Pu(:);

        end

        function P = getPuNodal(~,m)
            % =============================================================
            % function P = getPuNodal(~,m)
            %
            % returns prolongation operator for input of cell-width m
            % =============================================================

            v = @(x) x(:);

            mf = 2*m;
            % indices of coarse grid nodes (split by nodal, stg-i, c)
            ns   = prod(ones(length(m),1)*m+eye(length(m)),2); %staggered
            inc  = reshape(1:prod(m+1),m+1);
            is1c = reshape(prod(m+1)+(1:ns(1)),[m(1)+1 m(2) m(3)]);
            is2c = reshape(prod(m+1)+ns(1)+ (1:ns(2)),[m(1) m(2)+1 m(3)]);
            is3c = reshape(prod(m+1)+ns(1)+ns(2)+(1:ns(3)),[m(1) m(2) m(3)+1]);
            icc  = reshape(prod(m+1)+ns(1)+ns(2)+ns(3)+(1:prod(m)),m);
            % indices of fine grid nodes
            ns   = prod(ones(length(mf),1)*mf+eye(length(mf)),2); %staggered
            inf  = reshape(1:prod(mf+1),mf+1);
            is1f = reshape(prod(mf+1)+(1:ns(1)),[mf(1)+1 mf(2) mf(3)]);
            is2f = reshape(prod(mf+1)+ns(1)+ (1:ns(2)),[mf(1) mf(2)+1 mf(3)]);
            is3f = reshape(prod(mf+1)+ns(1)+ns(2)+(1:ns(3)),[mf(1) mf(2) mf(3)+1]);
            icf  = reshape(prod(mf+1)+ns(1)+ns(2)+ns(3)+(1:prod(mf)),mf);

            % allocate space
            I = []; J = []; W = [];

            % include existing nodes  (nodal-->nodal)
            ii = inf(1:2:end,1:2:end,1:2:end);
            jj = inc;
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % include existing nodes  (stg-1-->nodal)
            ii = inf(1:2:end,2:2:end,2:2:end);
            jj = is1c;
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % include existing nodes  (stg-2-->nodal)
            ii = inf(2:2:end,1:2:end,2:2:end);
            jj = is2c;
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % include existing nodes  (stg-3-->nodal)
            ii = inf(2:2:end,2:2:end,1:2:end);
            jj = is3c;
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % include existing nodes  (cc-->nodal)
            ii = inf(2:2:end,2:2:end,2:2:end);
            jj = icc;
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % average to get edge-stg-1
            ii   = inf(2:2:end,1:2:end,1:2:end);
            jj   = [reshape(inc(1:end-1,:,:),[],1); reshape(inc(2:end,:,:),[],1)];
            ww   = .5*ones(size(jj));
            I    = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % average to get edge-stg-2
            ii   = inf(1:2:end,2:2:end,1:2:end);
            jj   = [reshape(inc(:,1:end-1,:),[],1); reshape(inc(:,2:end,:),[],1)];
            ww   = .5*ones(size(jj));
            I    = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % average to get edge-stg-3
            ii   = inf(1:2:end,1:2:end,2:2:end);
            jj   = [reshape(inc(:,:,1:end-1),[],1); reshape(inc(:,:,2:end),[],1)];
            ww   = .5*ones(size(jj));
            I    = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % coarse --> fine(nodal only)
            P1  = sparse(I,J,W,prod(mf+1),icc(end));
            % allocate space
            I = []; J = []; W = [];

            % get face-stg-1
            ii  = is1f;
            jj  = [
                v(inf(1:end,1:end-1,1:end-1));
                v(inf(1:end,2:end  ,1:end-1));
                v(inf(1:end,2:end  ,2:end  ));
                v(inf(1:end,1:end-1,2:end  ));
                ];
            ww   = .25*ones(size(jj));
            I    = [I; ii(:);ii(:);ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % get face-stg-2
            ii  = is2f;
            jj  = [
                v(inf(1:end-1,1:end  ,1:end-1));
                v(inf(2:end  ,1:end  ,1:end-1));
                v(inf(2:end  ,1:end  ,2:end  ));
                v(inf(1:end-1,1:end  ,2:end  ));
                ];
            ww   = .25*ones(size(jj));
            I    = [I; ii(:);ii(:);ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            % get face-stg-3
            ii  = is3f;
            jj  = [
                v(inf(1:end-1,1:end-1,1:end));
                v(inf(2:end  ,1:end-1,1:end));
                v(inf(2:end  ,2:end  ,1:end));
                v(inf(1:end-1,2:end  ,1:end));
                ];
            ww   = .25*ones(size(jj));
            I    = [I; ii(:);ii(:);ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % get new cell-centers
            ii  = icf;
            jj  = [
                v(inf(1:end-1,1:end-1,1:end-1));
                v(inf(2:end  ,1:end-1,1:end-1));
                v(inf(2:end  ,2:end  ,1:end-1));
                v(inf(1:end-1,2:end  ,1:end-1));
                v(inf(1:end-1,1:end-1,2:end));
                v(inf(2:end  ,1:end-1,2:end));
                v(inf(2:end  ,2:end  ,2:end));
                v(inf(1:end-1,2:end  ,2:end));
                ];
            ww   = .125*ones(size(jj));
            I    = [I; repmat(ii(:),[8,1])]; J = [J; jj(:)]; W = [W; ww(:)];

            P2   = sparse(I-prod(mf+1),J,W,ns(1)+ns(2)+ns(3)+prod(mf),prod(mf+1));

            P    = [P1; P2*P1];
        end

        % ========== get methods ========================================
        function dx1 = get.dx1(this)
            if isempty(this.dx1_),
                [this.dx1_, this.dx2_, this.dx3_] = getGradientMatrixFEM(this,0);
            end
            dx1 = this.dx1_;
        end

        function dx2 = get.dx2(this)
            if isempty(this.dx2_),
                [this.dx1_, this.dx2_, this.dx3_] = getGradientMatrixFEM(this,0);
            end
            dx2 = this.dx2_;
        end

        function dx3 = get.dx3(this)
            if isempty(this.dx3_),
                [this.dx1_, this.dx2_, this.dx3_] = getGradientMatrixFEM(this,0);
            end
            dx3 = this.dx3_;
        end

        function GRAD = get.GRAD(this)
            if isempty(this.GRAD_),
                this.GRAD_ = [this.dx1;this.dx2;this.dx3];
            end
            GRAD = this.GRAD_;
        end
        function B = get.B(this)
            if isempty(this.B_),
                this.B_ = blkdiag(this.GRAD,this.GRAD,this.GRAD);
            end
            B = this.B_;
        end

        function P1 = get.P1(this)
            if isempty(this.P1_),
                A = speye(this.nnodes);
                this.P1_ = A(this.tri(:,1),:);
            end
            P1 = this.P1_;
        end
        function P2 = get.P2(this)
            if isempty(this.P2_),
                A = speye(this.nnodes);
                this.P2_ = A(this.tri(:,2),:);
            end
            P2 = this.P2_;
        end

        function P3 = get.P3(this)
            if isempty(this.P3_),
                A = speye(this.nnodes);
                this.P3_ = A(this.tri(:,3),:);
            end
            P3 = this.P3_;
        end

        function P4 = get.P4(this)
            if isempty(this.P4_),
                A = speye(this.nnodes);
                this.P4_ = A(this.tri(:,4),:);
            end
            P4 = this.P4_;
        end
        function PC = get.PC(this)
            if isempty(this.PC_),
                this.PC_ = (this.P1+this.P2+this.P3+this.P4)/4;
            end
            PC = this.PC_;
        end
        function P = get.P(this)
            if isempty(this.P_),
                this.P_ = this.getPuNodal(this.m);
            end
            P = this.P_;
        end
        function Pt = get.Pt(this)
            if isempty(this.Pt_),
                this.Pt_ = this.getPuNodal(this.m/2);
            end
            Pt = this.Pt_;
        end

        function mfdx1 = get.mfdx1(this)
            if isempty(this.mfdx1_),
                [this.mfdx1_, this.mfdx2_, this.mfdx3_] = getGradientMatrixFEM(this,1);
            end
            mfdx1 = this.mfdx1_;
        end

        function mfdx2 = get.mfdx2(this)
            if isempty(this.mfdx2_),
                [this.mfdx1_, this.mfdx2_, this.mfdx3_] = getGradientMatrixFEM(this,1);
            end
            mfdx2 = this.mfdx2_;
        end

        function mfdx3 = get.mfdx3(this)
            if isempty(this.mfdx3_),
                [this.mfdx1_, this.mfdx2_, this.mfdx3_] = getGradientMatrixFEM(this,1);
            end
            mfdx3 = this.mfdx3_;
        end

        function mfGRAD = get.mfGRAD(this)
            if isempty(this.mfGRAD_),
                this.mfGRAD_ = getGradientMatrixFEM(this,1);
            end
            mfGRAD = this.mfGRAD_;
        end

        function idx = get.boundaryIdx(this)
            if isempty(this.boundaryIdx_),
                mkvc = @(v) v(:);
                ns   = prod(ones(length(this.m),1)*this.m+eye(length(this.m)),2); %staggered

                % nodal points
                id = reshape(1:prod(this.m+1),this.m+1);
                idx = [ ...
                            mkvc(id([1,end],:,:));
                            mkvc(id(:,[1,end],:));
                            mkvc(id(:,:,[1,end]));
                            ];
                % stg-1 points
                id = prod(this.m+1) + reshape(1:prod(this.m+[1,0,0]),this.m+[1,0,0]);
                idx = [idx; mkvc(id([1,end],:,:)) ];
                % stg-2 points
                id = prod(this.m+1) +ns(1) + reshape(1:prod(this.m+[0,1,0]),this.m+[0,1,0]);
                idx = [idx; mkvc(id(:,[1,end],:)) ];
                % stg-3 points
                id = prod(this.m+1) +ns(1) + ns(2) + reshape(1:prod(this.m+[0,0,1]),this.m+[0,0,1]);
                idx = [idx; mkvc(id(:,:,[1,end])) ];

                this.boundaryIdx_ = unique( idx);
            end
            idx = this.boundaryIdx_;
        end

        function idx = get.boundaryProj(this)
            if isempty(this.boundaryProj_),
                idx = this.boundaryIdx;

                P = speye(size(this.xn,1));
                P = P(idx,:);
                P = kron(speye(3),P);

                this.boundaryProj_ = P;
            end
            idx = this.boundaryProj_;
        end

    end
end


text:
function moleFractionStats(modelT)
% Plots of mole fraction statistics.
%
% Plot a histogram of the number of metabolite species relevant between pH 5
% and 9, and a stacked bar chart of the mole fractions of reactants with
% significant (<0.99) distributions over more than one metabolite species.
%
% USAGE:
%
%    moleFractionStats(modelT)
%
% INPUT:
%    modelT:    structure with fields:
%
%                 * .S
%                 * .officialName
%                 * .mf
%
% .. Authors:
%       - Ronan M.T. Fleming
%       - Hulda SH Jul 7 2011 Enabled inclusion of up to 10 microspecies in bar graph

[nMets,nRxn]=size(modelT.S);

for m=1:nMets
    %nSpecies(m)=length(find(modelT.mf(m))>=0.001); % Imposed a 0.1% threshold on the abundance of species - Hulda
    nSpecies(m)=length(find(cell2mat(modelT.mf(m))>=0.001)); % Imposed a 0.1% threshold on the abundance of species - Hulda
    if (sum(cell2mat(modelT.mf(m)))-1)>1e-8
        %modelT.met(m)
        %modelT.met(m).mf
        error('Sum of mole fractions should be one')
    end
    if nnz(find(cell2mat(modelT.mf(m))<0))>1
        %modelT.met(m)
        %modelT.met(m).mf
        error('Mole fractions should be positive')
    end
end
[n, xout]=hist(nSpecies, [1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5]); % Uses as many bins as are needed, not just 3 - Hulda
n = n(1:find(n,1,'last'));
figure;
bar(0.5:1:(length(n) - 0.5),n);
title('Number of metabolite species for each reactant.')
ylabel('# Reactants','FontSize',14);
set(gca,'XTick',0.5:1:(length(n) - 0.5));
xTickLabels = {'One', 'Two', 'Three', 'Four', 'Five', 'Six', 'Seven', 'Eight', 'Nine', 'Ten'};
set(gca,'XTickLabel',xTickLabels(1:length(n)),'FontSize',14);

maxNSpecies=max(nSpecies);
%Y=zeros(nSpecies,maxNSpecies);
%Y=zeros(1,maxNSpecies);
Y=zeros(nMets,maxNSpecies);

ns=1;
p=1;
%sorted 1,2,3
while ns<=maxNSpecies
    for m=1:nMets
        if ns==nSpecies(m)
            %Y(p,1:nSpecies(m))=sort(modelT.mf{m}','descend');
            %Y(p,1:nSpecies(m))=cell2mat(modelT.mf(m));
            temp = cell2mat(modelT.mf(m));
            for i=1:nSpecies(m)
                Y(p,i)=temp(i);
            end
            p=p+1;
        end
    %     Y(m,1:nSpecies(m))=modelT.met(m).mf';
    end
    ns=ns+1;
end
 figure
 barh(Y,'stack')

Y=[];
%sorted by connectivity
significantMF=false(nMets,1);
Shat=sign(abs(modelT.S));
A=Shat*Shat';
metConnect=diag(A,0)';
size(metConnect);
[metConnectSorted,ix]=sort(metConnect);
p=1;
for m=1:nMets
    %only take cytoplasmic metabolites
    metAbbr=modelT.mets{ix(m)};
    if max(cell2mat(modelT.mf(ix(m))))<0.95 && strcmp('[c]',metAbbr(end-2:end))
%         Y(p,1:nSpecies(ix(m)))=sort(modelT.met(ix(m)).mf','descend');
        temp = cell2mat(modelT.mf(ix(m)));
        for i=1:length(temp)
            Y(p,i) = temp(i);
        end
        %Y(p,1:nSpecies(ix(m)))=modelT.mf(ix(m))';
        metInd(p)=ix(m);
%         metNames{p,1}=[modelT.met(ix(m)).officialName metAbbr(end-2:end)];
        metNames(p)=modelT.metNames(ix(m));
        p=p+1;
    end

end

figure;
left=0.1;
bottom=0.1;
width=2;
height=0.8;
%subplot(1,2,1,'Position',[left bottom width height])
%subplot(1,2,1,'Position',[left bottom width height])
%figure;
barh(Y,'stack');
set(gca,'XTick',[0:0.2:1]);
set(gca,'YTick',1:length(metInd));
set(gca,'YTickLabel',metNames);
ylim([0.5 length(metInd)-0.5]) ;
xlabel('Mole Fraction','FontSize',14);
set(gca,'FontSize',9)
title('Mole fraction statistics')
saveas(gcf,'moleFraction95','fig');

left=0.8;
%subplot(1,2,2,'Position',[left bottom width height])
figure;
barh(metConnect(metInd));
set(gca,'YTick',1:length(metInd));
set(gca,'YTickLabel',metNames);
ylim([0.5 length(metInd)-0.5]) ;
xlabel('Metabolite (ordered as mole fraction stat figure)','FontSize',14);
xlabel('Connectivity','FontSize',14);
set(gca,'FontSize',9);
title('Mole fraction statistics: connectivity')
saveas(gcf,'moleFraction95connectivity','fig');

author: opencobra

repo: cobratoolbox


text:
function [pp, mse] = fitSpline(tData, xData, tKnot)
% pp = fitSpline(tData, xData, tGrid)
%
% Fits a piece-wise cubic spline to {tData, xData} using the
% knot-points in tKnot.
%
% INPUTS:
%     tData = [1,nData] = time stamps for data (uniform spacing, monotonic!)
%     xData = [nFunc,nData] = function value at each time
%     tKnot = [1,nKnot] = time for each knot point (monotonic!)
%
%   OUTPUTS:
%     pp = struct to be evaluated by Matlab's ppval command
%
%   NOTES:   (strongly suggested)
%
%     tData([1,end]) == tKnot([1,end])
%

nFunc = size(xData,1);
nKnot = length(tKnot);

xGuess = interp1(tData',xData',tKnot')';

%Numerically estimate the local slope
vData = zeros(size(xData));
for i=1:nFunc
vData(i,:) = diffCenter(xData(i,:),mean(diff(tData)));
end

vGuess = interp1(tData',vData',tKnot')';
zGuess = reshape([xGuess;vGuess],2*nFunc*nKnot,1);
options = optimoptions('fminunc',...
    'MaxFunEvals',200*length(zGuess),...
    'Display','off',...
    'Algorithm','quasi-newton');

objFun = @(z)checkFit(z,tData,xData,tKnot,nFunc,nKnot);
zSoln = fminunc(objFun,zGuess,options);

% Check the solution and return the answer:
[mse, pp] = checkFit(zSoln,tData,xData,tKnot,nFunc,nKnot);

end

function [mse, pp] = checkFit(z,tData,xData,tKnot,nFunc,nKnot)

zKnot = reshape(z,2*nFunc,nKnot);
xKnot = zKnot(1:nFunc,:);
vKnot = zKnot((nFunc+1):end,:);

pp = pwch(tKnot,xKnot,vKnot);
xFit = ppval(pp,tData);
mse = sum(mean((xFit-xData).^2));

end

function dx = diffCenter(x,dt)
% dx = diffCenter(x,dt)
%
%   Computes the second-order finite difference approximation of x with
%   respect to t. A one-sided second order difference is used at the end
%   points, so size(dx) == size(x).
%
%   INPUTS:
%       x = [1 x n] vector of function values
%       dt = sampling period of x   (default = 1)
%
%   OUTPUTS:
%       d = dx/dt = first derivative of x wrt t
%
% See also: cumInt, diff

if nargin == 1
    dt = 1;
end

if length(dt) ~= 1
    error('Time-step (dt) must be a scalar');
end

n = length(x);
if n < 2
    error('length(x) must be at least 2');
elseif n == 2
    dx = diff(x)/dt;
else %Then enough points for second-order finite difference
    dx = zeros(size(x));
    Dx = diff(x);   %Matlab first-order finite difference
    dx(1) = (-3*x(1) + 4*x(2) - x(3))/(2*dt);
    dx(2:(end-1)) = (Dx(1:(end-1)) + Dx(2:end))/(2*dt);  %Mid-points are easy
    dx(end) = (x(end-2) - 4*x(end-1) + 3*x(end))/(2*dt);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function y = mean(w)
   %MEAN Average or mean value of waveform's data.
   %   Y = mean(waveform)
   %   returns a scalar containing the mean value of the waveform data
   %
   %   Input Arguments
   %       WAVEFORM: waveform object   N-DIMENSIONAL
   %
   %   Output
   %       Y: array of same size as WAVEFORM, with each element corresponding
   %          to the mean value of the matching waveform
   %
   %   NOTE: Values of NaN are ignored.
   %
   %   See also NANMEAN, WAVEFORM/MIN, WAVEFORM/MAX, WAVEFORM/MEDIAN, SORT.

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   y = nan(size(w));
   % if the statistics toolbox is installed, use the builtin nanmean function
   % to ignore NaN values during the variance calculation.
   if ~isempty(ver('stats'))
      for I = 1 : numel(w);
         y(I) = nanmean( w(I).data );
      end
   else
      % the statistics toolbox is not installed, so any nan values will have
      % to be dealt with (ignored) manually.
      for I = 1 : numel(w);
         d = w(I).data;
         d = d(~isnan(d)); %ignore NaN values
         if ~isempty(d)
            y(I) = mean( d );
         end
      end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [beta_gibbs,sigma_gibbs,favar,It,Bu]=favar_dogibbs(It,Bu,B,EPS,n,T,lags,data_endo,data_exo,const,favar,ar,arvar,lambda1,lambda3,lambda4,m,p,k,priorexo,Y,X,cband,Tstar)

%% the methodolgy closely follows Bernanke, Boivin, Eliasz (2005) and lends from the FAVAR model of Koop & Korobilis

% function [beta_gibbs sigma_gibbs]=ndgibbs(It,Bu,beta0,omega0,X,Y,y,Bhat,n,T,q)
% performs the Gibbs algorithm 1.5.2 for the normal-diffuse prior, and returns draws from posterior distribution
% inputs:  - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - vector 'beta0': vector of prior values for beta (defined in 1.3.4)
%          - matrix 'omega0': prior covariance matrix for the VAR coefficients (defined in 1.3.8)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 1.1.8)
%          - matrix 'Y': matrix of regressands for the VAR model (defined in 1.1.8)
%          - vector 'y': vectorised regressands for the VAR model (defined in 1.1.12)
%          - matrix 'Bhat': OLS VAR coefficients, in non vectorised form (defined in 1.1.9)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
% outputs: - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)

%% preliminary tasks
% initialise variables
nfactorvar=favar.nfactorvar;
numpc=favar.numpc;
favarX=favar.X(:,favar.plotX_index);
onestep=favar.onestep;
XY=favar.XY;

    % initial conditions XZ0~N(XZ0mean,XZ0var)
    favar.XZ0mean=zeros(n*lags,1);
    favar.XZ0var=favar.L0*eye(n*lags); %BBE set-up

L=favar.L;
Sigma=bear.nspd(favar.Sigma);
if onestep==1
indexnM=favar.indexnM;
end
XZ0mean=favar.XZ0mean;
XZ0var=favar.XZ0var;
favar_X=favar.X;
% load priors
L0=favar.L0*eye(n);
a0=favar.a0;
b0=favar.b0;
sigmahat=(1/T)*(EPS'*EPS);

% preallocation
beta_gibbs=zeros(size(B(:),1),It-Bu);
sigma_gibbs=zeros(size(sigmahat(:),1),It-Bu);
% X_gibbs=zeros(size(X(:),1),It-Bu);
% Y_gibbs=zeros(size(Y(:),1),It-Bu);
% FY_gibbs=zeros(size(data_endo(:),1),It-Bu);
L_gibbs=zeros(size(L(:),1),It-Bu);
R2_gibbs=zeros(size(favarX,2),It-Bu);

if onestep==0 %static factors in this case
    FY=data_endo;
    pbstring='two-step'; %string for the progress bar
elseif onestep==1
    pbstring='one-step'; %string for the progress bar
end

% state-space representation
if onestep==1
B_ss=[B';eye(n*(lags-1)) zeros(n*(lags-1),n)];
sigma_ss=[sigmahat zeros(n,n*(lags-1));zeros(n*(lags-1),n*lags)];
end

% create a progress bar
hbar = bear.parfor_progressbar(It,['Progress of the Gibbs sampler (',pbstring,').']);

%% start iterations
for ii=1:It
    if onestep==1
        % Sample latent factors using Carter and Kohn (1994)
        FY=bear.favar_kfgibbsnv(XY,XZ0mean,XZ0var,L,Sigma,B_ss,sigma_ss,indexnM);
        % demean generated factors
        FY=bear.favar_demean(FY);
        % Sample autoregressive coefficients B,in the twostep procedure FY is static, and we want to use updated B
        [B,~,~,X,~,Y]=bear.olsvar(FY,data_exo,const,lags);
        [arvar]=bear.arloop(FY,const,p,n);
    end

   % set 'prior' values (here, the dummy observations)
   [Y,X,Tstar]=bear.doprior(Y,X,n,m,p,Tstar,ar,arvar,lambda1,lambda3,lambda4,priorexo);
   % obtain posterior distribution parameters
   [Bcap,betacap,Scap,alphacap,phicap,alphatop]=bear.dopost(X,Y,Tstar,k,n);

% draw B from a matrix-variate student distribution with location Bcap, scale Scap and phicap and degrees of freedom alphatop
stationary=0;
while stationary==0
B=bear.matrixtdraw(Bcap,Scap,phicap,alphatop,k,n);
   [stationary]=bear.checkstable(B(:),n,lags,size(B,1)); %switches stationary to 0, if the draw is not stationary
end
if onestep==1
B_ss(1:n,:)=B';
end

% then draw sigma from an inverse Wishart distribution with scale matrix Scap and degrees of freedom alphacap (step 3)
sigma=bear.iwdraw(Scap,alphacap);
if onestep==1
sigma_ss(1:n,1:n)=sigma;
end

%% Sample Sigma and L
[Sigma,L]=bear.favar_SigmaL(Sigma,L,nfactorvar,numpc,onestep,n,favar_X,FY,a0,b0,T,lags,L0);

%% record the values if the number of burn-in iterations is exceeded
if ii>Bu
% values of vector beta
beta_gibbs(:,ii-Bu)=B(:);
% values of sigma (in vectorized form)
sigma_gibbs(:,ii-Bu)=sigma(:);

% save the factors and loadings
X_gibbs(:,ii-Bu)=X(:);
Y_gibbs(:,ii-Bu)=Y(:);
FY_gibbs(:,ii-Bu)=FY(:);
L_gibbs(:,ii-Bu)=L(:);

% compute R2 (Coefficient of Determination) for plotX variables (can be done after burn-in)
R2=bear.favar_R2(favarX,FY);
R2_gibbs(:,ii-Bu)=R2(:);

% compute posterior estimates, this is different here to the other prior
[beta_median,B_median,beta_std,beta_lbound,beta_ubound,sigma_median]=bear.doestimates(betacap,phicap,Scap,alphacap,alphatop,n,k,cband);
beta_median_gibbs(:,:,ii-Bu)=beta_median;
B_median_gibbs(:,:,ii-Bu)=B_median;
beta_std_gibbs(:,:,ii-Bu)=beta_std;
beta_lbound_gibbs(:,:,ii-Bu)=beta_lbound;
beta_ubound_gibbs(:,:,ii-Bu)=beta_ubound;
sigma_median_gibbs(:,:,ii-Bu)=sigma_median;
% if current iteration is still a burn iteration, do not record the result
else
end

% update progress by one iteration
hbar.iterate(1);

% go for next iteration
end

% in case we have thinning of the draws,
thin=abs(round(favar.thin)); % should be a positive integer
if thin~=1
    beta_gibbs=beta_gibbs(:,thin:thin:end);
    sigma_gibbs=sigma_gibbs(:,thin:thin:end);
    X_gibbs=X_gibbs(:,thin:thin:end);
    Y_gibbs=Y_gibbs(:,thin:thin:end);
    FY_gibbs=FY_gibbs(:,thin:thin:end);
    L_gibbs=L_gibbs(:,thin:thin:end);
    R2_gibbs=R2_gibbs(:,thin:thin:end);
    It=(1/thin)*It;
    Bu=(1/thin)*Bu;
end

% save in favar structure
favar.X_gibbs=X_gibbs;
favar.Y_gibbs=Y_gibbs;
favar.FY_gibbs=FY_gibbs;
favar.L_gibbs=L_gibbs;
favar.R2_gibbs=R2_gibbs;

favar.beta_median_gibbs=beta_median_gibbs;
favar.B_median_gibbs=B_median_gibbs;
favar.beta_std_gibbs=beta_std_gibbs;
favar.beta_lbound_gibbs=beta_lbound_gibbs;
favar.beta_ubound_gibbs=beta_ubound_gibbs;
favar.sigma_median_gibbs=sigma_median_gibbs;

% close progress bar
close(hbar);

author: european-central-bank

repo: BEAR-toolbox


text:
function [varargout]=tri2quad(varargin)

% function [Fq,Vq,Cq]=tri2quad(Ft,Vt,convertMethod)

%%
switch nargin
    case 2
        Ft=varargin{1};
        Vt=varargin{2};
        convertMethod=1;
    case 3
        Ft=varargin{1};
        Vt=varargin{2};
        convertMethod=varargin{3};
end

%Cope with 2D input
if size(Vt,2)==2
    Vt(:,3)=0;
end
%% Mid edge sets

edgeMat=[Ft(:,[1 2]); Ft(:,[2 3]);  Ft(:,[3 1])];

E_sort=sort(edgeMat,2); %Sorted edges matrix
[~,ind1,~]=unique(E_sort,'rows');
edgeMat=edgeMat(ind1,:);

numPoints = size(Vt,1);
numEdges = size(edgeMat,1);

% Get indices of the four edges associated with each face
A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
A = max(A,A'); %Copy symmetric

%Indices for A matrix
indA_12=Ft(:,1)+(Ft(:,2)-1)*numPoints;
indA_23=Ft(:,2)+(Ft(:,3)-1)*numPoints;
indA_31=Ft(:,3)+(Ft(:,1)-1)*numPoints;

%Get indices for vertex array
indV_12=full(A(indA_12));
indV_23=full(A(indA_23));
indV_31=full(A(indA_31));

%% Mid face

% Ft_sort=sort(Ft,2); %Sorted edges matrix
% [~,ind1,ind2]=unique(Ft_sort,'rows');
% Ft=Ft(ind1,:);

indV_midFace=(1:1:size(Ft,1))';
indOffset=numPoints+size(edgeMat,1);

indV_midFace123=(indV_midFace((1-1)*size(Ft,1)+(1:size(Ft,1))))+indOffset;

%% Create quad faces array

 Fq=[Ft(:,1) indV_12 indV_midFace123 indV_31;... %Corner quad 1
    indV_12  Ft(:,2) indV_23 indV_midFace123;... %Corner quad 2
     indV_midFace123 indV_23 Ft(:,3) indV_31;... %Corner quad 3
    ];

%% Create vertex arrays

%new mid-edge points
Vn=0.5*(Vt(edgeMat(:,1),:)+Vt(edgeMat(:,2),:));

switch convertMethod
    case 1
        %new mid-element points
        Vm=zeros(size(Ft,1),3);
        for q=1:1:size(Vt,2)
            X=Vt(:,q);
            if size(Ft,1)==1
                Vm(:,q)=mean(X(Ft)',2);
            else
                Vm(:,q)=mean(X(Ft),2);
            end
        end
    case 2
        %new mid-face points
        TR=triangulation(Ft,Vt);
        Vm = incenter(TR,(1:size(Ft,1))');
end

Vq=[Vt; Vn; Vm]; %Join point sets

CVq=[0*ones(size(Vt,1),1); 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1);];

%%

varargout{1}=Fq;
varargout{2}=Vq;
varargout{3}=CVq;

%%
% _*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:
% VL_DEMO_KMEANS_VS_BUILTIN

clear elaps energy ;

dimension  = 128 ;
numData    = 1000 ;
numCenters = 10 ;
numTrials  = 10 ;

for trial = 1:numTrials
  X = rand(dimension, numData) ;

  tic ;
  [idx, C, e] = kmeans(X', numCenters) ;
  e = sum(e) ;
  elaps.builtin(trial) = toc ;
  energy.builtin(trial) = e ;

  tic ;
  [idx_, C_, e] = vl_kmeans(X, numCenters, 'verbose') ;
  elaps.vl(trial) = toc ;
  energy.vl(trial) = e ;

  tic ;
  [idx_, C_, e] = vl_kmeans(X, numCenters, 'initialization', 'plusplus') ;
  elaps.vlpp(trial) = toc ;
  energy.vlpp(trial) = e ;
end

figure(1) ; clf ;
subplot(1,2,1) ; title('Energy') ;
mu = [mean(elaps.builtin) mean(elaps.vl) mean(elaps.vlpp)] ;
st = [std(elaps.builtin) std(elaps.vl) std(elaps.vlpp)] ;
bar(mu) ; hold on ;
errorbar(mu, st, 'linestyle', 'none', 'color', 'r', 'linewidth', 4) ;

subplot(1,2,2) ;
mu = [mean(energy.builtin) mean(energy.vl) mean(energy.vlpp)] ;
st = [std(energy.builtin) std(energy.vl) std(energy.vlpp)] ;
bar(mu) ; hold on ;
errorbar(mu, st, 'linestyle', 'none', 'color', 'r', 'linewidth', 4) ;

author: Cloud-CV

repo: object-proposals


text:
function [y,deriv] = identity_trans(w)
% This is an MV2DF . See MV2DF_API_DEFINITION.readme.
%
%  w --> w
%

if nargin==0
    test_this();
    return;
end

if isempty(w)
    map = @(w) w;
    y = linTrans(w,map,map);
    return;
end

if isa(w,'function_handle')
    f = identity_trans([]);
    y = compose_mv(f,w,[]);
    return;
end

f = identity_trans([]);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

function test_this()

f = identity_trans([]);
test_MV2DF(f,randn(5,1));

author: nesl


text:
function varargout = gabrielGraph(pts)
%GABRIELGRAPH  Gabriel Graph of a set of points.
%
%   EDGES = gabrielGraph(PTS)
%   Computes the Gabriel graph of the input set of points PTS. The Gabriel
%   graph is based on the euclidean Delaunay triangulation, and keeps only
%   edges whose circumcircle does not contain any other input point than
%   the edge extremities.
%
%   [NODES, EDGES] = gabrielGraph(PTS)
%   Also returns the initial set of points;
%
%   Example
%     pts = rand(100, 2);
%     edges = gabrielGraph(pts);
%     figure; drawPoint(pts);
%     hold on; axis([0 1 0 1]); axis equal;
%     drawGraph(pts, edges);
%
%   See also
%     graphs, drawGraph, delaunayGraph
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2012-01-22, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% compute Delaunay triangulation
if verLessThan('matlab', '8.1')
    % Code for versions before R2013a
    dt = DelaunayTri(pts); %#ok<DDELTRI>
else
    % Code for versions R2013a and later
    dt = delaunayTriangulation(pts);
end

% extract edges (N-by-2 array)
eds = dt.edges();

% radius of the circle circumscribed to each edge
rads = edgeLength([pts(eds(:,1), :) pts(eds(:,2), :)]) / 2;

% extract middle point of each edge
midPts = midPoint(pts(eds(:,1), :), pts(eds(:,2), :));

% distance between midpoints and all points
% closest points should be edge vertices
dists = minDistancePoints(midPts, pts);

% geometric tolerance (adapted to point set extent)
tol = max(max(pts) - min(pts)) * eps;

% keep only edges whose circumcircle does not contain any other point
keep = dists >= rads - tol;
edges = eds(keep, :);

% format output depending on number of output arguments
if nargout < 2
    varargout = {edges};
else
    varargout = {pts, edges};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Box volume maximization
% Boyd, Kim, Vandenberghe, and Hassibi, "A Tutorial on Geometric Programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% Maximizes volume of a box-shaped structure which has constraints
% on its total wall area, its total floor area, and which has lower
% and upper bounds on the aspect ratios. This leads to a GP:
%
%   maximize   h*w*d
%       s.t.   2(h*w + h*d) <= Awall, w*d <= Afloor
%              alpha <= h/w <= beta
%              gamma <= d/w <= delta
%
% where variables are the box height h, width w, and depth d.

% problem constants
alpha = 0.5; beta = 2; gamma = 0.5; delta = 2;

% varying parameters for an optimal trade-off curve
N = 10;
Afloor = logspace(1,3,N);
Awall  = [100 1000 10000];
opt_volumes = zeros(length(Awall),N);

disp('Computing optimal box volume for:')

% setup various GP problems with varying parameters
cvx_setpath
cvx_setspath
for k = 1:length(Awall)
  Awall_k = Awall(k);
  fprintf( 'Awall = %d:\n', Awall(k) );
  for n = 1:N
    % resolve the problem with varying parameters
    Afloor_n = Afloor(n);
    fprintf( '    Afloor = %7.2f: ', Afloor(n) );
    cvx_begin gp quiet
      variables h w d
      % objective function is the box volume
      maximize( h*w*d )
      subject to
        2*(h*w + h*d) <= Awall_k; %#ok
        w*d <= Afloor_n; %#ok
        alpha <= h/w <= beta; %#ok
        gamma <= d/w <= delta; %#ok
    cvx_end
    fprintf( 'max_volume = %3.2f\n', cvx_optval );
    opt_volumes(k,n) = cvx_optval;
  end
end
cvx_clearspath
cvx_clearpath

% plot the tradeoff curve
figure, clf
loglog(Afloor,opt_volumes(1,:), Afloor,opt_volumes(2,:), Afloor,opt_volumes(3,:));
xlabel('Afloor'); ylabel('V');

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

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

end

author: ecmadao

repo: Coding-Guide


text:
function [output binarizedContour] = CORF(I, sigma, t)
% CORF Contour detection based on a computational model of a simple cell.
%
% VERSION 22/04/2012
% CREATED BY: George Azzopardi and Nicolai Petkov, University of Groningen,
%             Johann Bernoulli Institute for Mathematics and Computer Science, Intelligent Systems
%
%If you use this script please cite the following paper:
%   George Azzopardi and Nicolai Petkov, "A CORF computational model of a
%     simple cell that relies on LGN input outperforms the Gabor function
%     model", 2012, DOI: 10.1007/s00422-012-0486-6
%
%   CORF achieves orientation selectivity by combining the output - at certain
%   positions with respect to the center of the CORF operator - of center-on
%   and center-off difference of Gaussians (DoG) functions by a weighted geometric mean.
%
%   CORF takes as input:
%      I -> intensity image
%      sigma -> the standard deviation of the outer Gaussian function of the
%               DoG operator
%      t -> high threshold used for hysteresis thresholding
%
%   CORF returns:
%      output -> maximum superposition of CORF responses that correspond to 12 orientations
%      binarizedContour -> A contour map that is obtaiend by first thinning
%                          output and then performs hysteresis thresholding
%                          with a high threshold t and a low threshold that
%                          is a fraction 0.5 of the given t.
%
%   Example: [o bc] = CORF(imread('rino.pgm'),2.5,0.3);
%
%   The image rino.pgm is taken from the RuG data set of 40 images of
%   natural scenes, which can be downloaded from:
%   http://www.cs.rug.nl/~imaging/databases/contour_database/contour_database.html

% configure CORF operator
operator = configureCORF(sigma,0.5);

% Set number of orientations
noriens = 12;
orienslist = 0:180/noriens:359;

% Preprocessing
inputImage = double(I);
if ndims(inputImage) == 3
    inputImage = double(rgb2gray(inputImage));
end
if max(inputImage(:)) > 1
    inputImage = inputImage ./ 255;
end

% Apply difference of Gaussians function.
padwidth = ceil(max(operator.params.rho));
DoGresponse(:,:,1) = applyDoG(inputImage, 0, padwidth, operator.params.sigma, operator.params.sigmaRatio, 0);
DoGresponse(:,:,2) = -DoGresponse(:,:,1);
DoGresponse(:,:,1) = DoGresponse(:,:,1) .* (DoGresponse(:,:,1) > 0);
DoGresponse(:,:,2) = DoGresponse(:,:,2) .* (DoGresponse(:,:,2) > 0);

% Apply the CORF oeprator at different orientations
data = [];
output = zeros(size(DoGresponse,1),size(DoGresponse,2),noriens);
for rot = 1:length(orienslist)
    rotatedOperator = operator;
    rotatedOperator.tuple(4,:) = rotatedOperator.tuple(4,:) + orienslist(rot)*pi/180;
    [output(:,:,rot) data] = getCORFresponse(DoGresponse,rotatedOperator,data);
end
output = output(padwidth+1:end-padwidth,padwidth+1:end-padwidth,:);

% Merge output of all orientations by using maximum superposition
[output, oriensMatrix] = calc_viewimage(output,1:size(output,3), orienslist*pi/180);

% Perform thinning
thinningOutput = thinning(output, oriensMatrix, 2);

%Perform hysteresis thresholding
t = t * max(thinningOutput(:));
binarizedContour = hysthresh(thinningOutput, t, 0.5*t);

function [output data] = getCORFresponse(DoGresponse,operator,data)

sz = [size(DoGresponse,1),size(DoGresponse,2)];

if isempty(data)
    data.params = [];
    data.tupleOutput = cell(0);
    data.location = [];
    data.paramsindex = 1;
    data.current.tupleOutput = zeros(sz(1),sz(2),size(operator.tuple,2));
end

% The weight vector is requried for the computation of the weighted
% geometric mean at the bottom of the function.
weightVector = zeros(1,size(operator.tuple,2));
sgm = max(operator.tuple(3,:)) / 3;

for i = 1:size(operator.tuple,2)
    polarity = operator.tuple(1,i);
    rho = operator.tuple(3,i);
    phi = operator.tuple(4,i);

    weightVector(i) = (exp(-rho^2/(2*sgm*sgm)));

    mem = find(ismember(data.params,[polarity rho],'rows'),1);
    if ~isempty(mem)
        % This tuple output is obtained by appropriate shifting of another
        % tuple output that was computed for the same values of polarity,
        % sigma and rho.
        [col row] = pol2cart(phi,rho);
        shiftrow = -(data.location{mem}(1)-fix(row));
        shiftcol = data.location{mem}(2)-fix(col);
        data.current.tupleOutput(:,:,i) = circshift(data.tupleOutput{mem},[shiftrow,shiftcol]);
    else
        DoG = DoGresponse(:,:,polarity+1);
        data.params(data.paramsindex,:) = [polarity rho];
        [col row] = pol2cart(phi,rho);

        r = (operator.params.d0 + operator.params.alpha*rho)/2;
        if r > 0
            smoothfilter = fspecial('gaussian',round([2*r+1,2*r+1]),r/3);
            data.current.tupleOutput(:,:,i) = conv2(DoG,smoothfilter,'same');
            data.current.tupleOutput(:,:,i) = circshift(data.current.tupleOutput(:,:,i),[fix(row),fix(-col)]);
        else
            data.current.tupleOutput(:,:,i) = DoG;
        end
        data.current.tupleOutput(:,:,i) = data.current.tupleOutput(:,:,i) .^ weightVector(i);

        % We use the following variables to reuse the computations obtained
        % for the same values of parameters: polarity, sigma and rho.
        data.tupleOutput{data.paramsindex} = data.current.tupleOutput(:,:,i);
        data.location{data.paramsindex} = [round(row) round(col)];
        data.paramsindex = data.paramsindex + 1;
    end
end

% compute the weighted geometric mean
output = prod(data.current.tupleOutput,3).^(1/sum(weightVector));

function operator = configureCORF(sigma, sigmaRatio)
% configureCORF: configure a CORF operator for the given parameters.

% The parameters alpha and d0 are set as reported in the above mentioned paper
operator.params.alpha = 0.9;
operator.params.d0 = 2;

operator.params.sigma = sigma;
operator.params.sigmaRatio = sigmaRatio;

% The following are the rho values that were used in the experiments
% reported in the above mentioned paper.
if sigma >= 1 && sigma < 2.5
    operator.params.rho = [14.38 6.9796 3.0310 1.4135];
elseif sigma >= 2.5 && sigma < 4
    operator.params.rho = [3.0515 6.1992 12.6488 24.62];
elseif sigma >= 4 && sigma <= 5
    operator.params.rho = [3.3021 4.7877 9.2467 18.08 34.43];
else
    error('The value of the parameter sigma is out of bounds');
end

maxRadius = ceil(max(operator.params.rho))+1;
stimulus = zeros((2*maxRadius));
stimulus(:,1:maxRadius) = 1;
center = [maxRadius maxRadius];

%Obtain the output of DoG function for the synthetic edge stimulus
DoGresponse(:,:,1) = applyDoG(stimulus, 0, maxRadius, sigma, sigmaRatio, 1);
DoGresponse(:,:,2) = -DoGresponse(:,:,1);
DoGresponse(:,:,1) = DoGresponse(:,:,1) .* (DoGresponse(:,:,1) > 0);
DoGresponse(:,:,2) = DoGresponse(:,:,2) .* (DoGresponse(:,:,2) > 0);

operator.tuple = [];
for r = 1:length(operator.params.rho)
    [polarity rho phi] = getTuple(DoGresponse,operator.params.rho(r),center);
    operator.tuple = [operator.tuple [polarity; repmat(sigma,1,length(polarity)); rho; phi]];
end

function [polarity rho phi] = getTuple(DoGresponse,radius,fp)
% Determine the values of the polar coordinates (rho,phi)

if radius <= 0
    error('Parameter radius must be greater than 0');
end

phi = []; polarity = [];
x = 1:360;

for pol = 1:size(DoGresponse,3)
    DoG = DoGresponse(:,:,pol);
    y = DoG(sub2ind(size(DoG),round(fp(1) + radius*cos(pi/2+x*pi/180)),round(fp(2) + radius*sin(pi/2+x*pi/180))));

    % Threshold low values
    y(y < 0.1*max(DoGresponse(:))) = 0;
    y = round(y*1000)/1000;

    BW     = bwlabel(imregionalmax(y));
    npeaks = max(BW(:));

    for i = 1:npeaks
        phi(end+1) = mean(x(BW == i)) * pi/180;
        polarity(end+1) = pol - 1;
    end
end
rho = repmat(radius,1,length(phi));

function output = applyDoG(inputImage,polarity,padwidth,sigma,sigmaRatio,crop)

% pad input image
paddedInputImage = padarray(inputImage,[padwidth padwidth],'both','symmetric');

% create DoG operator
sz = size(inputImage) + padwidth + padwidth;
g1 = fspecial('gaussian',sz,sigma);
g2 = fspecial('gaussian',sz,sigma*sigmaRatio);
if polarity == 1
    DoG = g2 - g1;
elseif polarity == 0
    DoG = g1 - g2;
else
    error('Polarity must be either 0 (on) or 1 (off)');
end

% compute DoG
output = fftshift(ifft2(fft2(DoG,sz(1),sz(2)) .* fft2(paddedInputImage)));

if crop == 1
    output = output(padwidth+1:end-padwidth,padwidth+1:end-padwidth);
end

function [result, oriensMatrix] = calc_viewimage(matrices, dispcomb, theta)
% VERSION 14/05/04
% CREATED BY: M.B. Wieling and N. Petkov, Groningen University,
%             Department of Computer Science, Intelligent Systems
%
% CALC_VIEWIMAGE: calculates the maximum-superposition of all the matrices stored
% in MATRICES (according to the L-infinity norm). It uses only the matrices for
% which the index is entered in DISPCOMB, e.g. if DISPCOMB contains the values
% 1,2,4: only the first, second and fourth matrix contained in MATRICES are used
% for the superposition. This method also calculates the orientationmatrix (ORIENSMATRIX)
% which stores the maximum orientation response of each point in the resulting matrix
% (RESULT) - for use in CALC_THINNING. A progressbar of the calculations is also shown.
%   CALC_VIEWIMAGE(MATRICES, DISPCOMB, THETA)
%   calculates the single viewing image according to the following parameters
%     MATRICES - the matrices which hold all the convolutions for each orientation
%     DISPCOMB  - the indexes of each matrix (in MATRICES) which should be used for the
%                 superposition
%     THETA - a list of all the orientations - to create ORIENSMATRIX

% initialize values
oriensMatrix = 0;
tmpMaxConv = -Inf;
result = -Inf;
cnt1 = 1;

if (size(dispcomb,2) == 1)
  result = matrices(:,:,dispcomb(1));
else

  % calculate the superposition (L-infinity norm)
  while (cnt1 <= size(dispcomb,2))
    % calculate the maximum orientation-response in each point (based on the absolute values)
    oriensMatrixtmp1 = (abs(matrices(:,:,dispcomb(cnt1))) > tmpMaxConv) .* theta(dispcomb(cnt1));
    oriensMatrixtmp2 = (abs(matrices(:,:,dispcomb(cnt1))) <= tmpMaxConv) .* oriensMatrix;
    oriensMatrix = oriensMatrixtmp1 + oriensMatrixtmp2;
    tmpMaxConv = max(abs(matrices(:,:,dispcomb(cnt1))), tmpMaxConv);

    % calculate the superposition
    result = max(result,abs(matrices(:,:,dispcomb(cnt1))));
    cnt1 = cnt1 + 1;
  end
end

function result = thinning(matrix, oriensMatrix, method)
% VERSION 27/02/05
% CREATED BY: M.B. Wieling and N. Petkov, Groningen University,
%             Department of Computer Science, Intelligent Systems
%
% THINNING: reduces the edges to a width of 1 pixel. This is done
%           by looking at the two pixels next to the current pixel.
%           the neighbourpixels are determined by the orientation of
%           the current point (this is stored in the exact same location
%           as in the matrix ORIENSMATRIX). A progressbar of the
%           calculations is also shown.
%   THINNING(MATRIX, ORIENSMATRIX, METHOD) thins the edge
%      MATRIX - the matrix which should be thinned
%      ORIENSMATRIX - the matrix which holds for every pixel of MATRIX the
%                     orientation (in the same position)
%      METHOD - the method of thinning: METHOD == 1: simple method, just
%               take the value of the nearest pixels to compare with.
%               e.g. if the orientation = 20 degrees, the points S & N are chosen
%               to compare to, if the orientation = 25 degrees the points SW & NE
%               are chosen (see below, P = current point)
%                             NW  N  NE
%                             W   P   E
%                             SW  S  SE
%               METHOD == 2: the values to compare with are calculated using
%               interpolation based on the surrounding pixels (NW, N, NE, E, SE, S, SW, W)
%
% corrected an error: <pi := <=pi

% create the coordinate-system
h = size(matrix,1); % height of image
w = size(matrix,2); % width of image
mx = max(h,w);
[xcoords, ycoords] = meshgrid(1:mx);
xcoords = xcoords(1:h, 1:w);
ycoords = ycoords(1:h, 1:w);

% set every value between 0 and pi
oriensMatrix = mod(oriensMatrix, pi);

% add a border of zeros to 'matrix' and 'oriensMatrix' to ease calculations.
matrixB(h+2,w+2) = 0;
matrixB(2:h+1,2:w+1) = matrix(:,:);
oriensMatrixB(h+2,w+2) = 0;
oriensMatrixB(2:h+1,2:w+1) = oriensMatrix(:,:);

if (method == 1) % simple thinning
  result = 0;
  for I=2:h+1 % the rows
    for J=2:w+1 % the columns
      orien = oriensMatrixB(I,J);

      % calculate dx and dy - this is correct as can be seen
      % in a picture
      dx = (orien < (3/8)*pi) - (orien > (5/8)*pi);
      dy = ((orien > (1/8)*pi) & (orien <= (1/2)*pi)) - ((orien > (1/2)*pi) & (orien < (7/8)*pi));

      % normally the same pixels would be checked if (dy and dx > 0) and (dy and dx < 0). because
      % different pixels should be checked with a gabor-orientation of 45 and 135 this difference should
      % be checked
      if (dy < 0) & (dx < 0)
       result(I,J) = ((matrixB(I,J) >= matrixB(I+dy, J+dx)) & (matrixB(I,J) >= matrixB(I-dy, J-dx))) * matrixB(I,J);
      else
       result(I,J) = ((matrixB(I,J) >= matrixB(I-dy, J+dx)) & (matrixB(I,J) >= matrixB(I+dy, J-dx))) * matrixB(I,J);
      end
    end
  end
else % linear thinning
  %hb = waitbar(0,'Applying linear thinning, please wait ... (Step 6/7)'); % display a progressbar
  result = 0;
  for I=2:h+1 % the rows
    for J=2:w+1 % the columns
      orien = oriensMatrixB(I,J);
      % get the values of the surrounding pixels
      north = matrixB(I-1, J);
      northeast = matrixB(I-1, J+1);
      east = matrixB(I, J+1);
      southeast = matrixB(I+1, J+1);
      south = matrixB(I+1, J);
      southwest = matrixB(I+1, J-1);
      west = matrixB(I, J-1);
      northwest = matrixB(I-1, J-1);

      % calculate the value of the points in one line (using interpolation)
      if (orien <= (1/4)*pi)
          fraction = orien/((1/4)*pi);
          pnt1 = (1-fraction) * east + (fraction) * northeast;
          pnt2 = (1-fraction) * west + (fraction) * southwest;
      elseif (orien <= (1/2)*pi)
          fraction = (orien-(1/4)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * northeast + (fraction) * north;
          pnt2 = (1-fraction) * southwest + (fraction) * south;
      elseif (orien <= (3/4)*pi)
          fraction = (orien-(1/2)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * north + (fraction) * northwest;
          pnt2 = (1-fraction) * south + (fraction) * southeast;
      elseif (orien <= pi)
          fraction = (orien-(3/4)*pi)/((1/4)*pi);
          pnt1 = (1-fraction) * northwest + (fraction) * west;
          pnt2 = (1-fraction) * southeast + (fraction) * east;
      else
        orien
      end
      result(I,J) = ( (matrixB(I,J) >= pnt1) & (matrixB(I,J) >= pnt2) ) * matrixB(I,J);
    end
    %waitbar(I/h); % update the progressbar
  end
  %close(hb);
end

% removing the borders
result = result(2:h+1, 2:w+1);

function bw = hysthresh(im, T1, T2)
% HYSTHRESH - Hysteresis thresholding
%
% Usage: bw = hysthresh(im, T1, T2)
%
% Arguments:
%             im  - image to be thresholded (assumed to be non-negative)
%             T1  - upper threshold value
%             T2  - lower threshold value
%
% Returns:
%             bw  - the thresholded image (containing values 0 or 1)
%
% Function performs hysteresis thresholding of an image.
% All pixels with values above threshold T1 are marked as edges
% All pixels that are adjacent to points that have been marked as edges
% and with values above threshold T2 are also marked as edges. Eight
% connectivity is used.
%
% It is assumed that the input image is non-negative
%
% Author: Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk @ csse uwa edu au   http://www.csse.uwa.edu.au/~pk
%
% December 1996  - Original version
% March    2001  - Speed improvements made (~4x)

%
% A stack (implemented as an array) is used to keep track of all the
% indices of pixels that need to be checked.
% Note: For speed the number of conditional tests have been minimised
% This results in the top and bottom edges of the image being considered to
% be connected.  This may cause some stray edges to be propagated further than
% they should be from the top or bottom.
%

if (T2 > T1 | T2 < 0 | T1 < 0)  % Check thesholds are sensible
  error('T1 must be >= T2 and both must be >= 0 ');
end

[rows, cols] = size(im);    % Precompute some values for speed and convenience.
rc = rows*cols;
rcmr = rc - rows;
rp1 = rows+1;

bw = im(:);                 % Make image into a column vector
pix = find(bw > T1);        % Find indices of all pixels with value > T1
npix = size(pix,1);         % Find the number of pixels with value > T1

stack = zeros(rows*cols,1); % Create a stack array (that should never
                            % overflow!)

stack(1:npix) = pix;        % Put all the edge points on the stack
stp = npix;                 % set stack pointer
for k = 1:npix
    bw(pix(k)) = -1;        % mark points as edges
end

% Precompute an array, O, of index offset values that correspond to the eight
% surrounding pixels of any point. Note that the image was transformed into
% a column vector, so if we reshape the image back to a square the indices
% surrounding a pixel with index, n, will be:
%              n-rows-1   n-1   n+rows-1
%
%               n-rows     n     n+rows
%
%              n-rows+1   n+1   n+rows+1

O = [-1, 1, -rows-1, -rows, -rows+1, rows-1, rows, rows+1];

while stp ~= 0            % While the stack is not empty
    v = stack(stp);         % Pop next index off the stack
    stp = stp - 1;

    if v > rp1 & v < rcmr   % Prevent us from generating illegal indices
			    % Now look at surrounding pixels to see if they
                            % should be pushed onto the stack to be
                            % processed as well.
       index = O+v;	    % Calculate indices of points around this pixel.
       for l = 1:8
	   ind = index(l);
	   if bw(ind) > T2   % if value > T2,
	       stp = stp+1;  % push index onto the stack.
	       stack(stp) = ind;
	       bw(ind) = -1; % mark this as an edge point
	   end
       end
    end
end

bw = (bw == -1);            % Finally zero out anything that was not an edge
bw = reshape(bw,rows,cols); % and reshape the image

author: Sable

repo: mcbench-benchmarks


text:
function [C,nc] = graph_coloring(A,nc)
  % GRAPH_COLORING Color a graph given an adjacency matrix. This heuristic is
  % designed to efficiently find a conservative coloring (i.e., using too many
  % colors). It does *not* attempt to find the minimal coloring. It is best used
  % when the provided nc is a few more than the optimal coloring. For example,
  % planar graphs can be colored with 4 colors, (slow) heuristics for finding the
  % optimal coloring will often find 6 colors, but this function will very
  % quickly find a 7-coloring (nc=7). On such a graph, if you passed nc=4, this
  % function will likely waste time attemping to find a 4-coloring, 5-coloring,
  % 6-coloring, and then very quickly find a 7-coloring. It is designed to give
  % up searching for parsimonious colorings and increase the number of colors
  % (throwing a warning).
  %
  % [C,nc] = graph_coloring(A,nc)
  %
  % Inputs:
  %   A  #V by #V adjacency matrix
  %   nc  desired number of colors (for planar/triangle meshes, 7 is a fast
  %     choice; for tetrahedral meshes, 13 appears to often be fast)
  % Outputs:
  %   C  #V list of color ids into (1:nc)
  %   nc  effective number of colors >= input nc

  n = size(A,1);
  [FI,FJ] = find(triu(A,1));
  % always marking the higher valence vertex as bad definitely makes things
  % worse 50x.
  %
  % if nc is much greater than optimal, always marking the lower valence vertex
  % doesn't much effect. But for nc close to optimal, this seems to make a very
  % big difference (300x)
  val = sum(A,2);
  swap = val(FI)<val(FJ);
  EI = FI.*swap + FJ.*~swap;
  EJ = FJ.*swap + FI.*~swap;

  onc = nc;
  % worst case just increase number of colors
  max_outer = 10;
  for outer = 1:max_outer
    % Might as well retry a few times for this nc
    for retry = 1:10
      C = ceil(nc*rand(n,1));
      for iter = 1:2*ceil(sqrt(size(A,1)))
        bad = unique(EI(C(EI)==C(EJ)));
        if numel(bad) == 0
          break;
        end
        C(bad) = nan;
        C(bad) = ceil(nc*rand(numel(bad),1));
      end
      if numel(bad) == 0
        break;
      end
    end
    if numel(bad) == 0
      break;
    end
    nc = nc+1;
    if nc>max_outer
      error('Failed to converge.\n');
    end
  end
  if nc ~= onc
    warning('Had to increasing number of colors')
  end

  %V2C = sparse(1:n,C,1,n,nc);
  % find all edges with the same color
  %A*V2C;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function IOU = IOU_with_GT(estimate, ground_truth)
    % strip off NaNs used to maintain all groundtruth in same matrix
    ground_truth = ground_truth(~isnan(ground_truth));
    % axis-aligned rectangle format
    if numel(ground_truth)==4
        if sum(estimate)==0
            IOU = 0;
        else
            IOU = bboxOverlapRatio(estimate, ground_truth);
        end
    else
        % rotated bounding box format (VOT)
        warning off;
        e_x = [estimate(1) estimate(1)+estimate(3) estimate(1)+estimate(3) estimate(1)];
        e_y = [estimate(2) estimate(2) estimate(2)+estimate(4) estimate(2)+estimate(4)];
        gt_x = ground_truth(1:2:7);
        gt_y = ground_truth(2:2:8);
        [intersection_x, intersection_y] = polybool('intersection',e_x,e_y,gt_x,gt_y);
        [union_x, union_y] = polybool('union',e_x,e_y,gt_x,gt_y);

        interseaction_area = polyarea(intersection_x, intersection_y);
        union_area = polyarea(union_x, union_y);

        IOU = interseaction_area / union_area;
        if isnan(IOU)
            IOU = 0;
        end
    end
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet