text:
classdef meshgen
    %   MESHGEN: Mesh generation class
    %   Handles input parameters to create a meshgen class object that can be
    %   used to build a msh class.
    %   Copyright (C) 2018  Keith Roberts & William Pringle
    %
    %   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/>.
    %
    %   Available options:
    %         ef            % edgefx class
    %         bou           % geodata class
    %         h0            % minimum edge length (optional if bou exists)
    %         bbox          % bounding box [xmin,ymin; xmax,ymax] (manual specification, no bou)
    %         proj          % structure containing the m_map projection info
    %         plot_on       % flag to plot (def: 1) or not (0)
    %         nscreen       % how many it to plot and write temp files (def: 5)
    %         itmax         % maximum number of iterations.
    %         pfix          % fixed node positions (nfix x 2 )
    %         egfix         % edge constraints
    %         outer         % meshing boundary (manual specification, no bou)
    %         inner         % island boundaries (manual specification, no bou)
    %         mainland      % the shoreline boundary (manual specification, no bou)
    %         fixboxes      % a flag that indicates which boxes will use fixed constraints
    %         memory_gb     % memory in GB allowed to use for initial rejector
    %         cleanup       % logical flag or string to trigger cleaning of topology (default is on).
    %         direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
    %         dj_cutoff     % the cutoff area fraction for disjoint portions to delete
    %         qual_tol      % tolerance for the accepted negligible change in quality
    %         enforceWeirs  % whether or not to enforce weirs in meshgen
    %         enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
    % delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen
    % improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality
    %         big_mesh      % set to 1 to remove the bou data from memory
    %
    properties
        fd            % handle to distance function
        fh            % handle to edge function
        h0            % minimum edge length
        edgefx        % edgefx class
        bbox          % bounding box [xmin,ymin; xmax,ymax]
        pfix          % fixed node positions (nfix x 2 )
        egfix         % edge constraints
        fixboxes      % a flag that indicates which boxes will use fixed constraints
        plot_on       % flag to plot (def: 1) or not (0)
        nscreen       % how many it to plot and write temp files (def: 5)
        bou           % geodata class
        ef            % edgefx class
        itmax         % maximum number of iterations.
        outer         % meshing boundary (manual specification, no bou)
        inner         % island boundaries (manual specification, no bou)
        mainland      % the shoreline boundary (manual specification, no bou)
        boubox        % the bbox as a polygon 2-tuple
        inpoly_flip   % used to flip the inpoly test to determine the signed distance.
        memory_gb     % memory in GB allowed to use for initial rejector
        cleanup       % logical flag or string to trigger cleaning of topology (default is on).
        direc_smooth  % logical flag to trigger direct smoothing of mesh in the cleanup
        dj_cutoff     % the cutoff area fraction for disjoint portions to delete
        grd = msh();  % create empty mesh class to return p and t in.
        big_mesh      % release bou data from memory
        qual          % mean, lower 3rd sigma, and the minimum element quality.
        qual_tol      % tolerance for the accepted negligible change in quality
        proj          % structure containing the m_map projection info
        anno          % Approx. Nearest Neighbor search object.
        annData       % datat contained with KD-tree in anno
        Fb            % bathymetry data interpolant
        enforceWeirs  % whether or not to enforce weirs in meshgen
        enforceMin    % whether or not to enfore minimum edgelength for all edgefxs
        delaunay_elim_on_exit % whether or not to run delaunay_elim on exit of meshgen
        improve_with_reduced_quality % whether or not to allow mesh improvements with decreases in mesh quality
    end

    methods

        % class constructor/default grd generation options
        function obj = meshgen(varargin)
            % Check for m_map dir
            M_MAP_EXISTS=0;
            if exist('m_proj','file')==2
              M_MAP_EXISTS=1 ;
            end
            if M_MAP_EXISTS~=1
              error('Where''s m_map? Chief, you need to read the user guide')
            end

            % Check for utilties dir
            UTIL_DIR_EXISTS=0 ;
            if exist('inpoly.m','file')
              UTIL_DIR_EXISTS=1;
            end
            if UTIL_DIR_EXISTS~=1
              error('Where''s the utilities directory? Chief, you need to read the user guide')
            end

            p = inputParser;
            % unpack options and set default ones, catch errors.

            defval = 0; % placeholder value if arg is not passed.
            % add name/value pairs
            addOptional(p,'h0',defval);
            addOptional(p,'bbox',defval);
            addOptional(p,'fh',defval);
            addOptional(p,'pfix',defval);
            addOptional(p,'egfix',defval);
            addOptional(p,'fixboxes',defval);
            addOptional(p,'inner',defval);
            addOptional(p,'outer',defval);
            addOptional(p,'mainland',defval);
            addOptional(p,'bou',defval);
            addOptional(p,'ef',defval);
            addOptional(p,'plot_on',defval);
            addOptional(p,'nscreen',defval);
            addOptional(p,'itmax',defval);
            addOptional(p,'memory_gb',1);
            addOptional(p,'cleanup',1);
            addOptional(p,'direc_smooth',1);
            addOptional(p,'dj_cutoff',0.25);
            addOptional(p,'big_mesh',defval);
            addOptional(p,'proj',defval);
            addOptional(p,'qual_tol',defval);
            addOptional(p,'enforceWeirs',1);
            addOptional(p,'enforceMin',1);
            addOptional(p,'delaunay_elim_on_exit',1);
            addOptional(p,'improve_with_reduced_quality',0);

            % parse the inputs
            parse(p,varargin{:});

            %if isempty(varargin); return; end
            % store the inputs as a struct
            inp = p.Results;

            % kjr...order these argument so they are processed in a predictable
            % manner. Process the general opts first, then the OceanMesh
            % classes...then basic non-critical options.
            inp = orderfields(inp,{'h0','bbox','enforceWeirs','enforceMin',...
                                   'delaunay_elim_on_exit','improve_with_reduced_quality',...
                                   'fh',...
                                   'inner','outer','mainland',...
                                   'bou','ef',... %<--OceanMesh classes come after
                                   'egfix','pfix','fixboxes',...
                                   'plot_on','nscreen','itmax',...
                                   'memory_gb','qual_tol','cleanup',...
                                   'direc_smooth','dj_cutoff',...
                                   'big_mesh','proj'});
            % get the fieldnames of the edge functions
            fields = fieldnames(inp);
            % loop through and determine which args were passed.
            % also, assign reasonable default values if some options were
            % not assigned.
            for i = 1 : numel(fields)
                type = fields{i};
                switch type
                    % parse aux options first
                    case('h0')
                        % Provide in meters
                        obj.h0 = inp.(fields{i});
                    case('fh')
                        if isa(inp.(fields{i}),'function_handle')
                            obj.fh = inp.(fields{i});
                        end
                        % can't check for errors here yet.
                    case('bbox')
                        obj.bbox= inp.(fields{i});
                        if iscell(obj.bbox)
                            % checking bbox extents
                            ob_min = obj.bbox{1}(:,1);
                            ob_max = obj.bbox{1}(:,2);
                            for ii = 2:length(obj.bbox)
                                if any(obj.bbox{ii}(:,1) < ob_min) || ...
                                        any(obj.bbox{ii}(:,2) > ob_max)
                                    error(['Outer bbox must contain all ' ...
                                        'inner bboxes: inner box #' ...
                                        num2str(ii) ' violates this'])
                                end
                            end
                        end

                        % if user didn't pass anything explicitly for
                        % bounding box make it empty so it can be populated
                        % from ef as a cell-array
                        if obj.bbox(1)==0
                            obj.bbox = [];
                        end
                    case('pfix')
                        obj.pfix= inp.(fields{i});
                        if obj.pfix(1)~=0
                            obj.pfix(:,:) = inp.(fields{i});
                        else
                            obj.pfix = [];
                        end
                        if obj.enforceWeirs
                            for j = 1 : length(obj.bou)
                                if  ~isempty(obj.bou{j}.weirPfix)
                                    obj.pfix = [obj.pfix ; obj.bou{j}.weirPfix];
                                end
                            end
                        end
                    case('egfix')
                        obj.egfix= inp.(fields{i});
                        if ~isempty(obj.egfix) && obj.egfix(1)~=0
                            obj.egfix = inp.(fields{i});
                        else
                            obj.egfix = [];
                        end
                        if obj.enforceWeirs
                            for j = 1 : length(obj.bou)
                                if ~isempty(obj.bou{j}.weirEgfix) && ~isempty(obj.egfix)
                                    obj.egfix = [obj.egfix ; obj.bou{j}.weirEgfix+max(obj.egfix(:))];
                                else
                                    obj.egfix =  obj.bou{j}.weirEgfix;
                                end
                            end
                        end
                        obj.egfix = renumberEdges(obj.egfix);
                    case('fixboxes')
                        obj.fixboxes= inp.(fields{i});
                    case('bou')
                        % got it from user arg
                        if obj.outer~=0, continue; end

                        obj.outer = {} ;
                        obj.inner = {} ;
                        obj.mainland = {} ;

                        obj.bou = inp.(fields{i});

                        % handle when not a cell
                        if ~iscell(obj.bou)
                            boutemp = obj.bou;
                            obj.bou = cell(1);
                            obj.bou{1} = boutemp;
                        end

                        % then the geodata class was provide, unpack
                        for ee = 1:length(obj.bou)
                            try
                                arg = obj.bou{ee} ;
                            catch
                                arg = obj.bou;
                            end
                            if isa(arg,'geodata')
                                obj.outer{ee} = obj.bou{ee}.outer;
                                obj.inner{ee} = obj.bou{ee}.inner;

                                % save bathy interpolant to meshgen
                                if ~isempty(obj.bou{ee}.Fb)
                                  obj.Fb{ee} = obj.bou{ee}.Fb ;
                                end

                                if ~isempty(obj.inner{ee}) && ...
                                   obj.inner{ee}(1)~= 0
                                    obj.outer{ee} = [obj.outer{ee};
                                        obj.inner{ee}];
                                end
                                obj.mainland{ee} = obj.bou{ee}.mainland;
                                obj.boubox{ee} = obj.bou{ee}.boubox;
                                obj.inpoly_flip{ee} = obj.bou{ee}.inpoly_flip;
                                if obj.big_mesh
                                    % release gdat's
                                    obj.bou{ee}.mainland= [];
                                    obj.bou{ee}.outer= [];
                                    if ~isempty(obj.bou{ee}.inner)
                                        obj.bou{ee}.inner= [];
                                    end
                                end
                            end
                        end

                    case('ef')
                        tmp = inp.(fields{i});
                        if isa(tmp, 'function_handle')
                            error('Please specify your edge function handle through the name/value pair fh');
                        end
                        obj.ef = tmp;

                        % handle when not a cell
                        if ~iscell(obj.ef)
                            eftemp = obj.ef;
                            obj.ef = cell(1);
                            obj.ef{1} = eftemp;
                        end

                        % Gather boxes from ef class.
                        for ee = 1 : length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.bbox{ee} = obj.ef{ee}.bbox;
                            end
                        end

                         % checking bbox extents
                        if iscell(obj.bbox)
                            ob_min = obj.bbox{1}(:,1);
                            ob_max = obj.bbox{1}(:,2);
                            for ii = 2:length(obj.bbox)
                                if any(obj.bbox{ii}(:,1) < ob_min) || ...
                                        any(obj.bbox{ii}(:,2) > ob_max)
                                    error(['Outer bbox must contain all ' ...
                                        'inner bboxes: inner box #' ...
                                        num2str(ii) ' violates this'])
                                end
                            end
                        end

                        % kjr 2018 June: get h0 from edge functions
                        for ee = 1:length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.h0(ee) = obj.ef{ee}.h0;
                            end
                        end

                        % kjr 2018 smooth the outer automatically
                        if length(obj.ef) > 1
                            % kjr 2020, ensure the min. sizing func is
                            % used
                            if obj.enforceMin
                                obj.ef = enforce_min_ef(obj.ef);
                            end
                            obj.ef = smooth_outer(obj.ef,obj.Fb);
                        end

                        % Save the ef interpolants into the edgefx
                        for ee = 1:length(obj.ef)
                            if isa(obj.ef{ee},'edgefx')
                                obj.fh{ee} = @(p)obj.ef{ee}.F(p);
                            end
                        end

                    case('plot_on')
                        obj.plot_on= inp.(fields{i});
                    case('big_mesh')
                        obj.big_mesh = inp.(fields{i});
                    case('nscreen')
                        obj.nscreen= inp.(fields{i});
                        if obj.nscreen ~=0
                            obj.nscreen = inp.(fields{i});
                            obj.plot_on = 1;
                        else
                            obj.nscreen = 5; % default
                        end
                    case('itmax')
                        obj.itmax= inp.(fields{i});
                        if obj.itmax ~=0
                            obj.itmax = inp.(fields{i});
                        else
                            obj.itmax = 100;
                            warning('No itmax specified, itmax set to 100');
                        end
                    case('qual_tol')
                        obj.qual_tol = inp.(fields{i});
                        if obj.qual_tol ~=0
                            obj.qual_tol = inp.(fields{i});
                        else
                            obj.qual_tol = 0.01;
                        end
                    case('inner')
                        if ~isa(obj.bou,'geodata')
                            obj.inner = inp.(fields{i});
                        end
                    case('outer')
                        if ~isa(obj.bou,'geodata')
                            obj.outer = inp.(fields{i});
                            if obj.inner(1)~=0
                                obj.outer = [obj.outer; obj.inner];
                            end
                        end
                    case('mainland')
                        if ~isa(obj.bou,'geodata')
                            obj.mainland = inp.(fields{i});
                        end
                    case('memory_gb')
                        if ~isa(obj.bou,'memory_gb')
                            obj.memory_gb = inp.(fields{i});
                        end
                    case('cleanup')
                        obj.cleanup = inp.(fields{i});
                        if isempty(obj.cleanup) || obj.cleanup == 0
                            obj.cleanup = 'none';
                        elseif obj.cleanup == 1
                            obj.cleanup = 'default';
                        end
                    case('dj_cutoff')
                        obj.dj_cutoff = inp.(fields{i});
                    case('direc_smooth')
                        obj.direc_smooth = inp.(fields{i});
                    case('proj')
                        obj.proj = inp.(fields{i});
                        % default CPP
                        if obj.proj == 0; obj.proj = 'equi'; end
                        if ~isempty(obj.bbox)
                            % kjr Oct 2018 use outer coarsest box for
                            % multiscale meshing
                            lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110;
                            lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110;
                            lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110;
                            lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110;
                        else
                            lon_mi = -180; lon_ma = 180;
                            lat_mi = -90; lat_ma = 90;
                        end
                        % Set up projected space
                        dmy = msh() ;
                        dmy.p(:,1) = [lon_mi; lon_ma];
                        dmy.p(:,2) = [lat_mi; lat_ma];
                        del = setProj(dmy,1,obj.proj) ;
                    case('enforceWeirs')
                        obj.enforceWeirs = inp.(fields{i});
                    case('enforceMin')
                        obj.enforceMin = inp.(fields{i});
                    case('delaunay_elim_on_exit')
                        obj.delaunay_elim_on_exit = inp.(fields{i});
                    case('improve_with_reduced_quality')
                        obj.improve_with_reduced_quality = inp.(fields{i});
                end
            end

            if isempty(varargin); return; end

            % error checking
            if isempty(obj.boubox) && ~isempty(obj.bbox)
                % Make the bounding box 5 x 2 matrix in clockwise order if
                % it isn't present. This case must be when the user is
                % manually specifying the PSLG.
                obj.boubox{1} = [obj.bbox(1,1) obj.bbox(2,1);
                                 obj.bbox(1,1) obj.bbox(2,2); ...
                                 obj.bbox(1,2) obj.bbox(2,2);
                                 obj.bbox(1,2) obj.bbox(2,1); ...
                                 obj.bbox(1,1) obj.bbox(2,1); NaN NaN];
            end
            if any(obj.h0==0), error('h0 was not correctly specified!'), end
            if isempty(obj.outer), error('no outer boundary specified!'), end
            if isempty(obj.bbox), error('no bounding box specified!'), end
            obj.fd = @dpoly;  % <-default distance fx accepts p and pv (outer polygon).
            % kjr build ANN object into meshgen
            obj = createANN(obj) ;

            global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST
            obj.grd.proj    = MAP_PROJECTION ;
            obj.grd.coord   = MAP_COORDS ;
            obj.grd.mapvar  = MAP_VAR_LIST ;

        end

        % Creates Approximate nearest neighbor objects on start-up
        function  obj = createANN(obj)

            box_vec = 1:length(obj.bbox);

            for box_num = box_vec
                if ~iscell(obj.outer)
                    dataset = obj.outer;
                    dataset(isnan(obj.outer(:,1)),:) = [];
                else
                    dataset = obj.outer{box_num};
                    dataset(isnan(obj.outer{box_num}(:,1)),:) = [];
                end
                if all(abs(obj.bbox{box_num}(1,:)) == 180)
                    % This line removes the line that can appear in the
                    % center for a global mesh
                    dataset(abs(dataset(:,1)) > 180-1e-6,:) = [];
                    dataset(abs(dataset(:,1)) < 1e-6,:) = [];
                end
                [dataset(:,1),dataset(:,2)] = m_ll2xy(dataset(:,1),dataset(:,2));
                dataset(isnan(dataset(:,1)),:) = [];
                dmy = ann(dataset');
                obj.anno{box_num} = dmy;
                obj.annData{box_num}=dataset;
            end
        end

        function  obj = build(obj)
            %DISTMESH2D 2-D Mesh Generator using Distance Functions.
            % Checking existence of major inputs
            %%
            warning('off','all')
            %%
            tic
            it = 1 ;
            Re = 6378.137e3;
            geps = 1e-3*min(obj.h0)/Re;
            deps = sqrt(eps);
            ttol=0.1; Fscale = 1.2; deltat = 0.1;
            delIT = 0 ; delImp = 2;
            imp = 10; % number of iterations to do mesh improvements (delete/add)

            % unpack initial points.
            p = obj.grd.p;
            if isempty(p)
                disp('Forming initial point distribution...');
                % loop over number of boxes
                for box_num = 1:length(obj.h0)
                    disp(['    for box #' num2str(box_num)]);
                    % checking if cell or not and applying local values
                    h0_l = obj.h0(box_num);
                    max_r0 = 1/h0_l^2;
                    if ~iscell(obj.bbox)
                        bbox_l = obj.bbox'; % <--we must tranpose this!
                    else
                        bbox_l = obj.bbox{box_num}'; % <--tranpose!
                    end
                    if ~iscell(obj.fh)
                        fh_l = obj.fh;
                    else
                        fh_l = obj.fh{box_num};
                    end
                    % Lets estimate the num_points the distribution will be
                    num_points = ceil(2/sqrt(3)*prod(abs(diff(bbox_l)))...
                                      /(h0_l/111e3)^2);
                    noblks = ceil(num_points*2*8/obj.memory_gb*1e-9);
                    len = abs(bbox_l(1,1)-bbox_l(2,1));
                    blklen = len/noblks;
                    st = bbox_l(1,1) ; ed = st + blklen; ns = 1;
                    %% 1. Create initial distribution in bounding box
                    %% (equilateral triangles)
                    for blk = 1 : noblks
                        if blk == noblks
                            ed = bbox_l(2,1);
                        end
                        ys = bbox_l(1,2);
                        ny = floor(1e3*m_lldist(repmat(0.5*(st+ed),2,1),...
                                   [ys;bbox_l(2,2)])/h0_l);
                        dy = diff(bbox_l(:,2))/ny;
                        ns = 1;
                        % start at lower left and make grid going up to
                        % north latitude
                        for ii = 1:ny
                            if st*ed < 0
                                nx = floor(1e3*m_lldist([st;0],...
                                     [ys;ys])/(2/sqrt(3)*h0_l)) + ...
                                     floor(1e3*m_lldist([0;ed],...
                                     [ys;ys])/(2/sqrt(3)*h0_l));
                            else
                                nx = floor(1e3*m_lldist([st;ed],...
                                     [ys;ys])/(2/sqrt(3)*h0_l));
                            end
                            ne = ns+nx-1;
                            if mod(ii,2) == 0
                                % no offset
                                x(ns:ne) = linspace(st,ed,nx);
                            else
                                % offset
                                dx = (ed-st)/nx;
                                x(ns:ne) = linspace(st+0.5*dx,ed,nx);
                            end
                            y(ns:ne) = ys;
                            ns = ne+1; ys = ys + dy;
                        end
                        st = ed;
                        ed = st + blklen;
                        p1 = [x(:) y(:)]; clear x y

                        %% 2. Remove points outside the region, apply the rejection method
                        p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:);     % Keep only d<0 points
                        r0 = 1./feval(fh_l,p1).^2;                          % Probability to keep point
                        p1 = p1(rand(size(p1,1),1) < r0/max_r0,:);          % Rejection method
                        p  = [p; p1];                                       % Adding p1 to p
                    end
                    if box_num == 1
                        % add points along the outermost polygon to fill
                        % outer extent more quickly.
                        outer_temp = obj.outer{1};
                        Inan = find(isnan(outer_temp(:,1)),1,'first');
                        p1 = outer_temp(1:Inan-1,:);
                        p1 = p1(feval(obj.fd,p1,obj,box_num) < geps,:);     % Keep only d<0 points
                        r0 = 1./feval(fh_l, p1).^2;                         % Probability to keep point
                        p1 = p1(rand(size(p1,1),1) < r0/max_r0,:);          % Rejection method
                        p = [p; p1];                                        % Adding p1 to p
                    end
                end
            else
                disp('User-supplied initial points!');
                obj.grd.b = [];
                h0_l = obj.h0(end); % finest h0 (in case of a restart of meshgen.build).
            end

            % remove pfix/egfix outside of desired subdomain
            nfix = size(obj.pfix,1);    % Number of fixed points
            negfix = size(obj.egfix,1); % Number of edge constraints
            if negfix > 0
                if length(obj.fixboxes)==1 && obj.fixboxes(1)==0
                  obj.fixboxes(1)=1 ;
                end
                pfixkeep = setdiff([1:nfix]',unique(obj.egfix(:)));
                % remove bars if midpoint is outside domain
                egfix_mid = (obj.pfix(obj.egfix(:,1),:) + ...
                             obj.pfix(obj.egfix(:,2),:))/2;
                for jj = 1 : length(obj.fixboxes)
                    if obj.fixboxes(jj)
                        iboubox = obj.boubox{jj};
                        inbar(:,jj) = inpoly(egfix_mid,iboubox(1:end-1,:));
                    end
                end
                inbar = sum(inbar,2) ;
                obj.egfix(~inbar,:) = [];
                tmppfix = obj.pfix([unique(obj.egfix(:)); pfixkeep],:);
                obj.pfix = tmppfix;
                obj.egfix = renumberEdges(obj.egfix);
                negfix = size(obj.egfix,1); % Number of edge constraints.
            end
            if nfix > 0
                if length(obj.fixboxes)==1 && obj.fixboxes(1)==0
                  obj.fixboxes(1)=1 ;
                end
                % remove pfix if outside domain
                for jj = 1 : length(obj.fixboxes)
                    if obj.fixboxes(jj)
                        inbox(:,jj) = inpoly(obj.pfix,obj.boubox{jj}(1:end-1,:));
                    end
                end
                inbox = sum(inbox,2) ;
                inbox(unique(obj.egfix(:))) = 1;
                obj.pfix(~inbox,:) = [];
                nfix = size(obj.pfix,1);    % Number of fixed points
            end
            if nfix >= 0, disp(['Using ',num2str(nfix),' fixed points.']);end
            if negfix > 0
                if max(obj.egfix(:)) > length(obj.pfix)
                    error('FATAL: egfix does index correcty into pfix.');
                end
                disp(['Using ',num2str(negfix),' fixed edges.']);
            end

            if ~isempty(obj.pfix); p = [obj.pfix; p]; end
            N = size(p,1); % Number of points N
            disp(['Number of initial points after rejection is ',num2str(N)]);
            %% Iterate
            pold = inf;                                                    % For first iteration
            if obj.plot_on >= 1
                clf,view(2),axis equal;
            end
            toc
            fprintf(1,' ------------------------------------------------------->\n') ;
            disp('Begin iterating...');
            while 1
                tic
                if ~mod(it,obj.nscreen) && delIT == 0
                    disp(['Iteration =' num2str(it)]) ;
                end

                % 3. Retriangulation by the Delaunay algorithm
                if max(sqrt(sum((p(1:size(pold,1),:)-pold).^2,2))/h0_l*111e3) > ttol         % Any large movement?
                    p = fixmesh(p);                                        % Ensure only unique points.
                    N = size(p,1); pold = p;                               % Save current positions
                    [t,p] = delaunay_elim(p,obj.fd,geps,0);                % Delaunay with elimination

                    if isempty(t)
                      disp('Exiting')
                      return
                    end

                    % Getting element quality and check "goodness"
                    if exist('pt','var'); clear pt; end
                    [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2));
                    tq = gettrimeshquan( pt, t);
                    mq_m = mean(tq.qm);
                    mq_l = min(tq.qm);
                    mq_s = std(tq.qm);
                    mq_l3sig = mq_m - 3*mq_s;
                    obj.qual(it,:) = [mq_m,mq_l3sig,mq_l];

                    % If mesh quality went down "significantly" since last iteration
                    % ..or..
                    % If not allowing improvements with reduction in quality
                    % then if the number of points significantly decreased
                    % due to a mesh improvement iteration, then rewind.
                    if ~mod(it,imp+1) && ((obj.qual(it,1) - obj.qual(it-1,1) < -0.10)  || ...
                        (~obj.improve_with_reduced_quality && ...
                        (N - length(p_before_improve))/length(p_before_improve) < -0.10))
                        disp('Mesh improvement was unsuccessful...rewinding...');
                        p = p_before_improve;
                        N = size(p,1);                                     % Number of points changed
                        pold = inf;
                        it = it + 1;
                        continue
                    else
                        N = size(p,1); pold = p;                           % Assign number of points and save current positions
                    end
                    % 4. Describe each bar by a unique pair of nodes.
                    bars = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])];           % Interior bars duplicated
                    bars = unique(sort(bars,2),'rows');                    % Bars as node pairs

                    % 5. Graphical output of the current mesh
                    if obj.plot_on >= 1 && (mod(it,obj.nscreen)==0 || it == 1)
                        cla,m_triplot(p(:,1),p(:,2),t)
                        m_grid
                        title(['Iteration = ',num2str(it)]);
                        if negfix > 0
                            m_plot(reshape(obj.pfix(obj.egfix,1),[],2)',...
                                 reshape(obj.pfix(obj.egfix,2),[],2)','r-')
                        end
                        if nfix > 0
                            m_plot(obj.pfix(:,1),obj.pfix(:,2),'b.')
                        end
                        plt = cell2mat(obj.boubox');
                        % reduce point spacing for asecthics
                        [plt2(:,2),plt2(:,1)] = my_interpm(plt(:,2),plt(:,1),0.1) ;
                        hold on ; axis manual
                        m_plot(plt2(:,1),plt2(:,2),'g','linewi',2)
                        drawnow
                    end
                end

                % Getting element quality and check goodness
                if exist('pt','var'); clear pt; end
                [pt(:,1),pt(:,2)] = m_ll2xy(p(:,1),p(:,2));
                tq = gettrimeshquan( pt, t);
                mq_m = mean(tq.qm);
                mq_l = min(tq.qm);
                mq_s = std(tq.qm);
                mq_l3sig = mq_m - 3*mq_s;
                obj.qual(it,:) = [mq_m,mq_l3sig,mq_l];

                % Improve the quality of triangles next to fixed edges by
                % deleting the point part of thin triangles without the fixed
                % point in it. Thin triangles have poor geometric quality <
                % 10%.
                if ~isempty(obj.egfix) && ~mod(it,delImp)
                    del = heal_fixed_edges(p,t,obj.egfix) ;
                    if ~isempty(del)
                        delIT = delIT + 1 ;
                        if delIT < 5
                            p(del,:)= [];
                            pold = inf;
                            disp(['Deleting ',num2str(length(del)),...
                                  ' points close to fixed edges']);
                            continue;
                        else
                            % Abandon strategy..if it will not terminate
                            disp('Moving to next iteration');
                        end
                    end
                    delIT = 0 ;
                 end

                % Termination quality, mesh quality reached is copacetic.
                qual_diff = mq_l3sig - obj.qual(max(1,it-imp),2);
                if ~mod(it,imp)
                    if abs(qual_diff) < obj.qual_tol
                        % Do the final elimination of small connectivity
                        if obj.delaunay_elim_on_exit
                            [t,p] = delaunay_elim(p,obj.fd,geps,1);
                        end
                        disp('Quality of mesh is good enough, exit')
                        close all;
                        break;
                    end
                end

                % Saving a temp mesh
                if ~mod(it,obj.nscreen) && delIT == 0
                    disp(['Number of nodes is ' num2str(length(p))])
                    disp(['Mean mesh quality is ' num2str(mq_m)])
                    disp(['Min mesh quality is ' num2str(mq_l)])
                    disp(['3rd sigma lower mesh quality is ' num2str(mq_l3sig)])
                    tempp = p; tempt = t;
                    save('Temp_grid.mat','it','tempp','tempt');
                    clearvars tempp tempt
                end

                % 6. Move mesh points based on bar lengths L and forces F
                barvec = pt(bars(:,1),:)- pt(bars(:,2),:);                 % List of bar vectors
                if strcmp(obj.grd.proj.name,'UTM')
                    % UTM is already in meters (useful for small domains)
                    L = sqrt(sum(barvec.^2,2))*Re;
                else
                    % Get spherical earth distances
                    long   = zeros(length(bars)*2,1);
                    lat    = zeros(length(bars)*2,1);
                    long(1:2:end) = p(bars(:,1),1);
                    long(2:2:end) = p(bars(:,2),1);
                    lat(1:2:end)  = p(bars(:,1),2);
                    lat(2:2:end)  = p(bars(:,2),2);
                    L = m_lldist(long,lat); L = L(1:2:end)*1e3;            % L = Bar lengths in meters
                end
                ideal_bars = 0.5*(pt(bars(:,1),:) + pt(bars(:,2),:));      % Used to determine what bars are in bbox
                [ideal_bars(:,1),ideal_bars(:,2)] = ...                    % needs to be in non-projected
                    m_xy2ll(ideal_bars(:,1),ideal_bars(:,2));              % coordinates
                hbars = 0*ideal_bars(:,1);

                for box_num = 1:length(obj.h0)                             % For each bbox, find the bars that are in and calculate
                    if ~iscell(obj.fh)                                     % their ideal lengths.
                        fh_l = obj.fh;
                    else
                        fh_l = obj.fh{box_num};
                    end
                    h0_l = obj.h0(box_num);
                    if box_num > 1
                        h0_l = h0_l/111e3;                                 % create buffer to evalulate fh between nests
                        iboubox = obj.boubox{box_num}(1:end-1,:) ;
                        inside = inpoly(ideal_bars,iboubox) ;
                    else
                        inside = true(size(hbars));
                    end
                    hbars(inside) = feval(fh_l,ideal_bars(inside,:));      % Ideal lengths
                end

                L0 = hbars*Fscale*median(L)/median(hbars);                 % L0 = Desired lengths using ratio of medians scale factor
                LN = L./L0;                                                % LN = Normalized bar lengths

                % Mesh improvements (deleting and addition)
                if ~mod(it,imp)
                    p_before_improve = p;
                    nn = []; pst = [];
                    if abs(qual_diff) < imp*obj.qual_tol && ...
                        (obj.improve_with_reduced_quality || qual_diff > 0)
                        % Remove elements with small connectivity
                        nn = get_small_connectivity(p,t);
                        disp(['Deleting ' num2str(length(nn)) ' due to small connectivity'])

                        % Remove points that are too close (< LN = 0.5)
                        if any(LN < 0.5)
                            % Do not delete pfix too close.
                            nn1 = setdiff(reshape(bars(LN < 0.5,:),[],1),[(1:nfix)']);
                            disp(['Deleting ' num2str(length(nn1)) ' points too close together'])
                            nn = unique([nn; nn1]);
                        end

                        % Split long edges however many times to
                        % better lead to LN of 1
                        if any(LN > 2)
                            nsplit = floor(LN);
                            nsplit(nsplit < 1) = 1;
                            adding = 0;
                            % Probably we can just split once?
                            for jj = 2:2
                                il = find(nsplit >= jj);
                                xadd = zeros(length(il),jj-1);
                                yadd = zeros(length(il),jj-1);
                                for jjj = 1 : length(il)
                                    deltax = (p(bars(il(jjj),2),1)- p(bars(il(jjj),1),1))/jj;
                                    deltay = (p(bars(il(jjj),2),2)- p(bars(il(jjj),1),2))/jj;
                                    xadd(jjj,:) = p(bars(il(jjj),1),1) + (1:jj-1)*deltax;
                                    yadd(jjj,:) = p(bars(il(jjj),1),2) + (1:jj-1)*deltay;
                                end
                                pst = [pst; xadd(:) yadd(:)];
                                adding = numel(xadd) + adding;
                            end
                            disp(['Adding ',num2str(adding) ,' points.'])
                        end
                    end
                    if ~isempty(nn) || ~isempty(pst)
                        % Doing the actual subtracting and add
                        p(nn,:)= [];
                        p = [p; pst];
                        pold = inf;
                        it = it + 1;
                        continue;
                    end
                end

                F    = (1-LN.^4).*exp(-LN.^4)./LN;                         % Bessens-Heckbert edge force
                Fvec = F*[1,1].*barvec;

                Ftot = full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));
                Ftot(1:nfix,:) = 0;                                        % Force = 0 at fixed points
                pt = pt + deltat*Ftot;                                     % Update node positions

                [p(:,1),p(:,2)] = m_xy2ll(pt(:,1),pt(:,2));

                %7. Bring outside points back to the boundary
                d = feval(obj.fd,p,obj,[],1); ix = d > 0;                  % Find points outside (d>0)
                ix(1:nfix) = 0;
                if sum(ix) > 0
                    pn = p(ix,:) + deps;
                    dgradx = (feval(obj.fd,[pn(:,1),p(ix,2)],obj,[])...%,1)...
                              -d(ix))/deps; % Numerical
                    dgrady = (feval(obj.fd,[p(ix,1),pn(:,2)],obj,[])...%,1)...
                              -d(ix))/deps; % gradient
                    dgrad2 = dgradx.^+2 + dgrady.^+2;
                    p(ix,:) = p(ix,:) - [d(ix).*dgradx./dgrad2,...
                                         d(ix).*dgrady./dgrad2];
                end

                % 8. Termination criterion: Exceed itmax
                it = it + 1 ;

                if ( it > obj.itmax )
                    % Do the final deletion of small connectivity
                    if obj.delaunay_elim_on_exit
                        [t,p] = delaunay_elim(p,obj.fd,geps,1);
                    end
                    disp('too many iterations, exit')
                    close all;
                    break ;
                end
                toc
            end
            %%
            warning('on','all')
            %%
            disp('Finished iterating...');
            fprintf(1,' ------------------------------------------------------->\n') ;

            %% Doing the final cleaning and fixing to the mesh...
            % Clean up the mesh if specified
            if ~strcmp(obj.cleanup,'none')
                % Put the mesh class into the grd part of meshgen and clean
                obj.grd.p = p; obj.grd.t = t;
                [obj.grd,qout] = clean(obj.grd,obj.cleanup,...
                                       'nscreen',obj.nscreen,'djc',obj.dj_cutoff,...
									    'pfix',obj.pfix);
                obj.grd.pfix = obj.pfix ;
				obj.grd.egfix= obj.egfix ;
                obj.qual(end+1,:) = qout;
            else
                % Fix mesh on the projected space
                [p(:,1),p(:,2)] = m_ll2xy(p(:,1),p(:,2));
                [p,t] = fixmesh(p,t);
                [p(:,1),p(:,2)] = m_xy2ll(p(:,1),p(:,2));
                % Put the mesh class into the grd part of meshgen
                obj.grd.p = p; obj.grd.t = t;
                obj.grd.pfix = obj.pfix ;
                obj.grd.egfix= obj.egfix ;
            end

            % Check element order, important for the global meshes crossing
            % -180/180 boundary
            obj.grd = CheckElementOrder(obj.grd);

            if obj.plot_on
                figure; plot(obj.qual,'linewi',2);
                hold on
                % plot the line dividing cleanup and distmesh
                plot([it it],[0 1],'--k')
                xticks(1:5:obj.itmax);
                xlabel('Iterations'); ylabel('Geometric element quality');
                title('Geometric element quality with iterations');
                set(gca,'FontSize',14);
                legend('q_{mean}','q_{mean}-q_{3\sigma}', 'q_{min}','Location','best');
                grid minor
            end
            return;
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            % Auxiliary subfunctions %
            %%%%%%%%%%%%%%%%%%%%%%%%%%

            function [t,p] = delaunay_elim(p,fd,geps,final)
                % Removing mean to reduce the magnitude of the points to
                % help the convex calc
                if exist('pt1','var'); clear pt1; end
                [pt1(:,1),pt1(:,2)] = m_ll2xy(p(:,1),p(:,2));
                if isempty(obj.egfix)
                    p_s  = pt1 - repmat(mean(pt1),[N,1]);
                    TR   = delaunayTriangulation(p_s);
                else
                    TR   = delaunayTriangulation(pt1(:,1),pt1(:,2),obj.egfix);
                    pt1  = TR.Points;
                end
                for kk = 1:final+1
                    if kk > 1
                        % Perform the following below upon exit from the mesh
                        % generation algorithm
                        nn = get_small_connectivity(pt1,t);
                        nn1 = heal_fixed_edges(pt1,t,obj.egfix) ;
                        nn = unique([nn; nn1]) ;
                        TR.Points(nn,:) = [];
                        pt1(nn,:) = [];
                    end
                    t = TR.ConnectivityList;
                    pmid = squeeze(mean(reshape(pt1(t,:),[],3,2),2));      % Compute centroids
                    [pmid(:,1),pmid(:,2)] = m_xy2ll(pmid(:,1),pmid(:,2));  % Change back to lat lon
                    t    = t(feval(fd,pmid,obj,[]) < -geps,:);             % Keep interior triangles
                    if kk == 1
                        % Deleting very straight triangles
                        tq_n = gettrimeshquan( pt1, t);
                        bad_ele = any(tq_n.vang < 1*pi/180 | ...
                                      tq_n.vang > 179*pi/180,2);
                        t(bad_ele,:) = [];
                    end
                end
                if length(pt1) ~= length(p)
                    clear p
                    [p(:,1),p(:,2)] = m_xy2ll(pt1(:,1),pt1(:,2));
                end
            end

            function nn = get_small_connectivity(p,t)
                % Get node connectivity (look for 4)
                [~, enum] = VertToEle(t);
                % Make sure they are not boundary nodes
                bdbars = extdom_edges2(t, p);
                bdnodes = unique(bdbars(:));
                I = find(enum <= 4);
                nn = setdiff(I',[(1:nfix)';bdnodes]);                      % and don't destroy pfix or egfix!
                return;
            end

            function del = heal_fixed_edges(p,t,egfix)
                % kjr april2019
                % if there's a triangle with a low geometric quality that
                % contains a fixed edge, remove the non-fixed vertex
                % perform this on every other iteration to allow non-fixed
                % points to create equilateral triangles nearby the locked
                % edge.
                % returns points IDs that should be deleted.
                TR = triangulation(t,p) ;
                elock = edgeAttachments(TR,egfix) ;
                tq = gettrimeshquan(p,t);
                elock = unique(cell2mat(elock'));
                dmy = elock(tq.qm(elock) < 0.25);
                badtria = t(dmy,:);
                del     = badtria(badtria > nfix) ;
            end

        end % end distmesh2d_plus

    end % end methods

end % end class

author: CHLNDDEV


text:
function varargout=rdmat(varargin)
%
% [tm,signal,Fs,siginfo]=rdmat(recordName)
%
% Import a signal in physical units from a *.mat file generated by WFDB2MAT.
% Required Parameters:
%
% recorName
%       String specifying the name of the *.mat file.
%
% Outputs are:
%
% tm
%       A Nx1 array of doubles specifying the time in seconds.
% signal
%       A NxM matrix of doubles contain the signals in physical units.
% Fs
%       A 1x1 integer specifying the sampling frequency in Hz for the entire record.
%siginfo
%       A LxN cell array specifying the signal siginfo. Currently it is a
%       structure with the following fields:
%
%        siginfo.Units
%        siginfo.Baseline
%        siginfo.Gain
%        siginfo.Description
%
% NOTE:
%       You can use the WFDB2MAT command in order to convert the record data into a *.mat file,
%       which can then be loaded into MATLAB/Octave's workspace using the LOAD command.
%       This sequence of procedures is quicker (by several orders of magnitude) than calling RDSAMP.
%       The LOAD command will load the signal data in raw units, use RDMAT to load the signal in physical units.
%
% KNOWN LIMITATIONS:
%       This function currently does support several of the features described
%       in the WFDB record format (such as multiresolution signals) :
%          http://www.physionet.org/physiotools/wag/header-5.htm
%       If you are not sure that the record (or database format) you are reading is
%       supported, you can do an integrity check by comparing the output with RDSAMP:
%
%       [tm,signal,Fs,siginfo]=rdmat('200m');
%       [tm2,signal2]=rdsamp('200m');
%       if(sum(abs(signal-signal2)) !=0);
%          error('Record not compatible with RDMAT');
%       end
%
%
% Written by Ikaro Silva, 2014
% Last Modified: November 26, 2014
% Version 1.2
%
% Since 0.9.7
%
% %Example:
% wfdb2mat('mitdb/200')
%tic;[tm,signal,Fs,siginfo]=rdmat('200m');toc
%tic;[tm2,signal2]=rdsamp('200m');toc
% sum(abs(signal-signal2))
%
%
% See also RDSAMP, WFDB2MAT

%endOfHelp

%Set default pararameter values
inputs={'recordName'};
defGain=200; %Default value for missing gains
wfdbNaN=-32768; %This should be the case for all WFDB signal format types currently supported by RDMAT

for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};'])
    end
end

outputs={'tm','val','Fs','siginfo'};
fid = fopen([recordName, '.hea'], 'rt');
if(fid==-1)
    error(['Could not open file: ' recordName '.hea !'])
end

%Following the documentation described in :
%http://www.physionet.org/physiotools/wag/header-5.htm
%to parse the header file

%Skip any comment lines
str=fgetl(fid);
while(strcmp(str(1),'#'))
    str=fgetl(fid);
end

%Process Record Line Info
info=textscan(str,'%s %u %f %u %s %s');
M=info{2}; %Number of signals present
Fs=info{3};

%Process Signal Specification lines. Assumes no comments between lines.
siginfo=[];
for m = 1:M
    str=fgetl(fid);
    info=textscan(str,'%s %s %s %u %u %f %u %u %s');
    gain=info{3}{:};

    %Get Signal Units if present
    ind=strfind(gain,'/');
    if(~isempty(ind))
        siginfo(m).Units=gain(ind+1:end);
        gain=gain(1:ind-1);
    end

    %Get Signal Baseline if present
    ind=strfind(gain,'(');
    if(~isempty(ind))
        ind2=strfind(gain,')');
        siginfo(m).Baseline=str2num(gain(ind+1:ind2-1));
        gain=gain(1:ind-1);
    else
        %If Baseline is missing, set it equal to ADC Zero
        adc_zero=info{5};
        if(~isempty(adc_zero))
            siginfo(m).Baseline=double(adc_zero);
        else
            error('Could not obtain signal baseline');
        end
    end

    %Get Signal Gain
    gain=str2num(gain);
    if(gain==0)
        %Set gain to default value in this case
        gain=defGain;
    end
    siginfo(m).Gain=double(gain);

    %Get Signal Descriptor
    siginfo(m).Description=info{9}{:};

end
fclose(fid);

load([recordName '.mat']);
val(val==wfdbNaN)= NaN;
for m = 1:M
    %Convert from digital units to physical units.
    % Mapping should be similar to that of rdsamp.c:
    % http://www.physionet.org/physiotools/wfdb/app/rdsamp.c
    val(m, :) = (val(m, :) - siginfo(m).Baseline ) / siginfo(m).Gain;
end

%Reshape to the Toolbox's standard format
val=val';

%Generate time vector
N=size(val,1);
tm =linspace(0,(N-1)/Fs,N);

for n=1:nargout
    eval(['varargout{n}=' outputs{n} ';'])
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [g,dgdx,dgdP] = g_ttest(x,P,u,in)

% transformations
%%% normal priors on phi_1
sP = P;
dsPdP = ones(numel(P),1);
%%% sparse laplace priors on phi_2
[sP(2),dsPdP(2)] = VBA_sparsifyPrior (P(2));

% prediction
g = [sP(1) ;
     sP(1) + sP(2) ];

% derivatives
dgdP = diag(dsPdP);
dgdP(1,2) = 1*dsPdP(1);
dgdx = [];

% % prediction
% g = [sP(1) - sP(2)/2;
%      sP(1) + sP(2)/2 ];
%
% % derivatives
% dgdP = zeros(2);
% dgdP(1,1) = dsPdP(1);
% dgdP(2,1) = -(1/2)*dsPdP(2);
% dgdP(1,2) = dsPdP(1);
% dgdP(2,2) = +(1/2)*dsPdP(2);
% dgdx = [];

end

author: MBB-team

repo: VBA-toolbox


text:
clear
addpath(genpath('code'));
addpath(genpath('manopt'));

author: geopavlakos


text:
function [A] = ft2A(ft)
% Convert length from feet to angstroms.
% Chad A. Greene 2012
A = ft*3048000000;

author: Sable

repo: mcbench-benchmarks


text:
function	[out,dout]=rodrigues(in)

% RODRIGUES	Transform rotation matrix into rotation vector and viceversa.
%
%		Sintax:  [OUT]=RODRIGUES(IN)
% 		If IN is a 3x3 rotation matrix then OUT is the
%		corresponding 3x1 rotation vector
% 		if IN is a rotation 3-vector then OUT is the
%		corresponding 3x3 rotation matrix
%

%%
%%		Copyright (c) March 1993 -- Pietro Perona
%%		California Institute of Technology
%%

%% ALL CHECKED BY JEAN-YVES BOUGUET, October 1995.
%% FOR ALL JACOBIAN MATRICES !!! LOOK AT THE TEST AT THE END !!

%% BUG when norm(om)=pi fixed -- April 6th, 1997;
%% Jean-Yves Bouguet

%% Add projection of the 3x3 matrix onto the set of special ortogonal matrices SO(3) by SVD -- February 7th, 2003;
%% Jean-Yves Bouguet

[m,n] = size(in);
%bigeps = 10e+4*eps;
bigeps = 10e+20*eps;

if ((m==1) && (n==3)) | ((m==3) && (n==1)) %% it is a rotation vector
   theta = norm(in);
   if theta < eps
      R = eye(3);

      %if nargout > 1,

      dRdin = [0 0 0;
	       0 0 1;
	       0 -1 0;
	       0 0 -1;
	       0 0 0;
	       1 0 0;
	       0 1 0;
	       -1 0 0;
          0 0 0];

       %end;

   else
      if n==length(in)  in=in'; end; 	%% make it a column vec. if necess.

	 %m3 = [in,theta]

	 dm3din = [eye(3);in'/theta];

	 omega = in/theta;

	 %m2 = [omega;theta]

	 dm2dm3 = [eye(3)/theta -in/theta^2; zeros(1,3) 1];

	 alpha = cos(theta);
	 beta = sin(theta);
	 gamma = 1-cos(theta);
	 omegav=[[0 -omega(3) omega(2)];[omega(3) 0 -omega(1)];[-omega(2) omega(1) 0 ]];
	 A = omega*omega';

	 %m1 = [alpha;beta;gamma;omegav;A];

	 dm1dm2 = zeros(21,4);
	 dm1dm2(1,4) = -sin(theta);
	 dm1dm2(2,4) = cos(theta);
	 dm1dm2(3,4) = sin(theta);
	 dm1dm2(4:12,1:3) = [0 0 0 0 0 1 0 -1 0;
	                     0 0 -1 0 0 0 1 0 0;
			     0 1 0 -1 0 0 0 0 0]';

         w1 = omega(1);
	 w2 = omega(2);
	 w3 = omega(3);

	 dm1dm2(13:21,1) = [2*w1;w2;w3;w2;0;0;w3;0;0];
	 dm1dm2(13: 21,2) = [0;w1;0;w1;2*w2;w3;0;w3;0];
	 dm1dm2(13:21,3) = [0;0;w1;0;0;w2;w1;w2;2*w3];

	 R = eye(3)*alpha + omegav*beta + A*gamma;

	 dRdm1 = zeros(9,21);

	 dRdm1([1 5 9],1) = ones(3,1);
	 dRdm1(:,2) = omegav(:);
	 dRdm1(:,4:12) = beta*eye(9);
	 dRdm1(:,3) = A(:);
	 dRdm1(:,13:21) = gamma*eye(9);

	 dRdin = dRdm1 * dm1dm2 * dm2dm3 * dm3din;

      end;
      out = R;
      dout = dRdin;

      %% it is prob. a rot matr.
   elseif ((m==n) & (m==3) & (norm(in' * in - eye(3)) < bigeps)...
	    & (abs(det(in)-1) < bigeps))
      R = in;

      % project the rotation matrix to SO(3);
      [U,S,V] = svd(R);
      R = U*V';

      tr = (trace(R)-1)/2;
      dtrdR = [1 0 0 0 1 0 0 0 1]/2;
      theta = real(acos(tr));

      if sin(theta) >= 1e-5,

	 dthetadtr = -1/sqrt(1-tr^2);

	 dthetadR = dthetadtr * dtrdR;
	 % var1 = [vth;theta];
	 vth = 1/(2*sin(theta));
	 dvthdtheta = -vth*cos(theta)/sin(theta);
	 dvar1dtheta = [dvthdtheta;1];

	 dvar1dR =  dvar1dtheta * dthetadR;

	 om1 = [R(3,2)-R(2,3), R(1,3)-R(3,1), R(2,1)-R(1,2)]';

	 dom1dR = [0 0 0 0 0 1 0 -1 0;
	       0 0 -1 0 0 0 1 0 0;
	       0 1 0 -1 0 0 0 0 0];

	 % var = [om1;vth;theta];
	 dvardR = [dom1dR;dvar1dR];

	 % var2 = [om;theta];
	 om = vth*om1;
	 domdvar = [vth*eye(3) om1 zeros(3,1)];
	 dthetadvar = [0 0 0 0 1];
	 dvar2dvar = [domdvar;dthetadvar];

	 out = om*theta;
	 domegadvar2 = [theta*eye(3) om];

	 dout = domegadvar2 * dvar2dvar * dvardR;

      else
	 if tr > 0; 			% case norm(om)=0;

	    out = [0 0 0]';

	    dout = [0 0 0 0 0 1/2 0 -1/2 0;
		  0 0 -1/2 0 0 0 1/2 0 0;
		  0 1/2 0 -1/2 0 0 0 0 0];
	 else 				% case norm(om)=pi; %% fixed April 6th

	    out = theta * (sqrt((diag(R)+1)/2).*[1;2*(R(1,2:3)>=0)'-1]);
	    %keyboard;

	    if nargout > 1,
	       fprintf(1,'WARNING!!!! Jacobian domdR undefined!!!\n');
		 	dout = NaN*ones(3,9);
	    end;
	 end;
      end;

   else
      error('Neither a rotation matrix nor a rotation vector were provided');
   end;

return;

%% test of the Jacobians:

%%%% TEST OF dRdom:
om = randn(3,1);
dom = randn(3,1)/1000000;

[R1,dR1] = rodrigues(om);
R2 = rodrigues(om+dom);

R2a = R1 + reshape(dR1 * dom,3,3);

gain = norm(R2 - R1)/norm(R2 - R2a)

%%% TEST OF dOmdR:
om = randn(3,1);
R = rodrigues(om);
dom = randn(3,1)/10000;
dR = rodrigues(om+dom) - R;

[omc,domdR] = rodrigues(R);
[om2] = rodrigues(R+dR);

om_app = omc + domdR*dR(:);

gain = norm(om2 - omc)/norm(om2 - om_app)

%%% OTHER BUG: (FIXED NOW!!!)

omu = randn(3,1);
omu = omu/norm(omu)
om = pi*omu;
[R,dR]= rodrigues(om);
[om2] = rodrigues(R);
[om om2]

%%% NORMAL OPERATION

om = randn(3,1);
[R,dR]= rodrigues(om);
[om2] = rodrigues(R);
[om om2]

author: strawlab

repo: MultiCamSelfCal


text:
% All the installation instructions are in a separate .md file named vonBertalanffy.md in docs/source/installation
% Setup the paths to the data, scripts and functions
% Check if this COBRA toolbox extension is in the Matlab path
stm_dir = which('setupThermoModel.m');
cc_dir = which('addThermoToModel.m');

% test if call to python2 works
[status,result] = system('python2 --version');
if status~=0
    disp(result)
    %check the current version of python that is running
    [status,result] = system('python --version');
    disp('Current version of python that is running:')
    disp(result)

    %find out what python is running and change to another version of
    %python (after a new version of python is installe)
    %     which python
    %     ls -al /usr/bin/python
    %     sudo rm -rf /usr/bin/python
    %     sudo ln -s /usr/bin/python3.7 /usr/bin/python

    error('python2 should be installed')
end

%test if call to cxcalc is installed
[status,result] = system('cxcalc');
if status ~= 0
    if status==127
        disp(result)
        error('Check that ChemAxon Marvin Beans is working licence is working and cxcalc is on the system path.')
    else
        disp(result)
        setenv('PATH', [getenv('PATH') ':/usr/local/bin/chemaxon/bin'])%RF
        setenv('PATH', [getenv('PATH') ':/opt/ChemAxon/MarvinBeans/bin/'])
        setenv('CHEMAXON_LICENSE_URL',[getenv('HOME') '/.chemaxon/license.cxl'])
        [status,result] = system('cxcalc');
        if status ~= 0
            disp(result)
            error('Check that ChemAxon Marvin Beans is installed, licence is working and cxcalc is on the system path.')
        end
    end
end

[FILEPATH,NAME,EXT] = fileparts(which('testVonBertalanffy.m'));
%test if call to cxcalc is working
[status,result] = system(['cxcalc pka -a 20 -b 20 -M false ' FILEPATH filesep 'test.inchi']);
if contains(result,'FAILED')
    disp(result)
    error('Check that ChemAxon Marvin Beans is installed, but may not have a working licence.')
else
    disp('ChemAxon Marvin Beans is installed and working.')
end

%check if call to babel works
[status,result] = system('ldd /usr/bin/babel');
if ~isempty(strfind(result,'MATLAB'))
    disp(result)
    fprintf('%s\n','babel must depend on the system libstdc++.so.6 not the one from MATLAB')
    fprintf('%s\n','Trying to edit the ''LD_LIBRARY_PATH'' to make sure that it has the correct system path before the Matlab path!')
    setenv('LD_LIBRARY_PATH',['/usr/lib/x86_64-linux-gnu:' getenv('LD_LIBRARY_PATH')]);
    fprintf('%s\n','The solution will be arch dependent');
end

[status,result] = system('babel');
if status ~= 0
    setenv('LD_LIBRARY_PATH',['/usr/lib/x86_64-linux-gnu:' getenv('LD_LIBRARY_PATH')]);
end

[status,result] = system('babel');
if status ~= 0
    error('Check that OpenBabel is installed and on the system path.')
end

% Add extension to path if it is not already there
if isempty(stm_dir) || isempty(cc_dir)
    vonB_dir = pwd;

    c = {'ComponentContribution'; 'SetupThermoModel'; 'testing'}; % Check that current directory contains all necessary subdirectories
    d = dir;
    c_test = {d.name}';

    if all(ismember(c,c_test))
        addpath(genpath(vonB_dir)); % Add to path
        savepath; % Save new path
    else
        error(['Current directory must be */vonBertalanffy, which should contain the following files and directories:\n' sprintf('%s\n',c{:})]);
    end
end

author: opencobra

repo: cobratoolbox


text:
function [ r, seed ] = r8_uniform_01 ( seed )

%*****************************************************************************80
%
%% R8_UNIFORM_01 returns a unit pseudorandom R8.
%
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2^31 - 1 )
%      r8_uniform_01 = seed / ( 2^31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R8_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.  SEED should not be 0.
%
%    Output, real R, a random value between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8_UNIFORM_01 - Fatal error!' );
  end

  seed = floor ( seed );

  seed = mod ( seed, i4_huge );

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  r = seed * 4.656612875E-10;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function str = showdetail(item)

% function str = showdetail(item)
% Generic showdetail function for cfg_item classes. It displays the
% name, tag, class and default function call for this item..
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

str{1,1} = sprintf('Name   : %s', item.name);
str{2,1} = sprintf('Tag    : %s', gettag(item));
if ~isempty(item.def)
    str{3} = sprintf('Default: %s', func2str(item.def));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function b= ImageDenoising(img, level)
     distorted_img = imnoise(img,'gaussian',0,level);

     yRGB = im2double(distorted_img);
     % Generate the same seed used in the experimental results of [1]
     randn('seed', 0);
     % Standard deviation of the noise --- corresponding to intensity
     %  range [0,255], despite that the input was scaled in [0,1]
     sigma = 25;
     % Add the AWGN with zero mean and standard deviation 'sigma'
     zRGB = yRGB + (sigma/255)*randn(size(yRGB));
     % Denoise 'zRGB'. The denoised image is 'yRGB_est', and 'NA = 1'
     %  because the true image was not provided
     [~, yRGB_est] = CBM3D(1, zRGB, sigma);
     % Compute the putput PSNR
%      PSNR = 10*log10(1/mean((yRGB(:)-yRGB_est(:)).^2))
     % show the noisy image 'zRGB' and the denoised 'yRGB_est'
%      figure; imshow(min(max(zRGB,0),1));
%      figure; imshow(min(max(yRGB_est,0),1));
     b = uint8(yRGB_est*255);

end

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
function [Xcs,Ycs,Zcs]=contourRegulariseAngle(varargin)

switch nargin
    case 2
        Vcs=varargin{1};
        np=varargin{2};
        thetaStart=0; %Default angle
        interpMethod='pchip'; %Default method
    case 3
        Vcs=varargin{1};
        np=varargin{2};
        thetaStart=varargin{3};
        interpMethod='pchip'; %Default method
            case 4
        Vcs=varargin{1};
        np=varargin{2};
        thetaStart=varargin{3};
        interpMethod=varargin{4}; %Default method
    otherwise
        error('Wrong number of input arguments');
end

%% CONTROL PARAMETERS
closeLoopOpt=1;

numContours=numel(Vcs);

%Allocate coordinate matrices
Xcs=nan(numContours,np); Ycs=nan(numContours,np); Zcs=nan(numContours,np);
startDefined=0;
for q=1:1:numContours

    %Join goups if present N.B. assumes they belong to the same curve!
    Vs=[];
    for qGroup=1:1:numel(Vcs{q})
        Vss=Vcs{q}{qGroup};
        if ~isempty(Vss)
            if isPolyClockwise(Vss)
                Vss=flipud(Vss);
            end
            Vs=[Vs; Vss];
        end
    end

    %Resample curve so they have the same number of points
    if ~isempty(Vs)

        numDigitsMerge=6-numOrder(mean(sum(diff(Vs,[],1).^2,2)));
        [~,ind1,ind2]=unique(pround(Vs,numDigitsMerge),'rows');

        Vs=Vs(ismember(1:size(Vs,1),ind1),:);

%         figure; plotV(Vs,'k.-'); hold on;
%         plotV(Vs(1,:),'r.');
%         plotV(Vs(end,:),'b.-');

        [Vs]=evenlySampleCurve(Vs,np,interpMethod,closeLoopOpt);

        %Reorder curve so that start points are close
         %%
%         if startDefined==0 %After 1 contour has been added
%             [thetaAll,~,~] = cart2pol(Vs(:,1),Vs(:,2),Vs(:,3));
%             thetaDiff=thetaAll-thetaStart;
%             [~,indMin]=min(thetaDiff);
%             V_start=Vs(indMin,:);
%             %Reorder first curve
%             Vs=[Vs(indMin:size(Vs,1),:); Vs(1:indMin-1,:)];
%             startDefined=1;
%         else
%             %Reorder curves so start and end points match closely
%             D=sqrt(sum((Vs-V_start(ones(1,size(Vs,1),1),:)).^2,2));
%             [~,indMin]=min(D);
%             if indMin~=1
%                 Vs=[Vs(indMin:size(Vs,1),:); Vs(1:indMin-1,:)];
%             end
%             V_start=Vs(1,:);
%         end

        %%
        meanVs=mean(Vs,1);
        Vsc=Vs-meanVs(ones(size(Vs,1),1),:);
        [thetaAll,~,~] = cart2pol(Vsc(:,1),Vsc(:,2),Vsc(:,3));
        thetaDiff=thetaAll-thetaStart;
        [~,indMin]=min(thetaDiff);

        if indMin~=1
            Vs=[Vs(indMin:size(Vs,1),:); Vs(1:indMin-1,:)];
        end

        %%

        %Store coordinates
        [Xcs(q,:),Ycs(q,:),Zcs(q,:)]=getColumns(Vs);
%         any(isnan(Zcs(q,:)))
    end
end

%%
% _*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:
x = [9 2.140641477955609996536345; ...
	2.5 0.7031566406452431872257; ...
	0.1 -10.42375494041107679516822; ...
	7e-4 -1429.147493371120205005198; ...
	7e-5 -14286.29138623969227538398; ...
	7e-6 -142857.7200612932791081972; ...
	2e-6 -500000.5772123750382073831; ...
	1e-6 -1000000.577214019968668068; ...
	7e-7 -1428572.005785942019703646; ...
	-0.5 .03648997397857652055902367; ...
	-1.1 10.15416395914385769902271 ...
	];
for i = 1:rows(x)
	actual = digamma(x(i,1));
	expected = x(i,2);
	e = abs(actual - expected)/expected;
	if e > 1e-12
		error(sprintf('digamma(%g) = %g should be %g', x(i,1), actual, expected));
	end
end
if digamma(-1) ~= -Inf
	error('digamma(-1) should be -Inf');
end
if digamma(0) ~= -Inf
  error('digamma(0) should be -Inf');
end
if ~isnan(digamma(-Inf))
  error('digamma(-Inf) should be NaN');
end

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed