text:
classdef ReasenbergDeclusterClass < ZmapFunction
    % Reasenberg Declustering codes
    % made from code originally from A.Allman that has chopped up and rebuilt

    properties
        taumin  duration    = days(1)   % look ahead time for not clustered events
        taumax  duration    = days(10)  % maximum look ahead time for clustered events
        P                   = 0.95      % confidence level that this is next event in sequence
        xk                  = 0.5       % is the factor used in xmeff

        % effective lower magnitude cutoff for catalog.
        % During clusteres, it is raised by a factor xk*cmag1
        xmeff               = 1.5

        rfact               = 10        % factor for interaction radius for dependent events
        err                 = 1.5       % epicenter error
        derr                = 2         % depth error, km
        %declustRoutine      = "ReasenbergDeclus"
        declusteredCatalog   ZmapCatalog
        replaceSequenceWithEquivMainshock   logical = false

         % if empty, clustering details will not be saved to workspace
        clusterDetailsVariableName          char    = 'cluster_details'

        % if empty, catalog will not be saved to workspace
        declusteredCatalogVariableName      char    = 'declustered_catalog'
        memorizeOriginalCatalog             logical = true
    end

    properties(Constant)
        PlotTag = "ReasenbergDecluster"

        ParameterableProperties = ["taumin", "taumax", "P",...
            "xk","xmeff","rfact","err","derr",..."declustRoutine",...
            "replaceSequenceWithEquivMainshock",...
            "clusterDetailsVariableName",...
            "declusteredCatalogVariableName",...
            "memorizeOriginalCatalog"];

        References = ['Paul Reasenberg (1985) ',...
            '"Second -order Moment of Central California Seismicity"',...
            ', JGR, Vol 90, P. 5479-5495.'];

        % interaction formula bundles assumptions about stress drop into a mag-dependent formula
        InteractFormula = struct('Reasenberg1985', @(m) 0.011 .* 10.^ (0.4 .* m),...
            'WellsCoppersmith1994', @(m) 0.01 * 10 .^ (0.5 * m))
    end

    methods
        function obj = ReasenbergDeclusterClass(catalog, varargin)
            % ReasenbergDeclusterClass

            obj@ZmapFunction(catalog);

            report_this_filefun();
            obj.parseParameters(varargin);
            obj.clusterDetailsVariableName = matlab.lang.makeValidName(obj.clusterDetailsVariableName);
            obj.declusteredCatalogVariableName = matlab.lang.makeValidName(obj.declusteredCatalogVariableName);
            obj.StartProcess();
        end

        function InteractiveSetup(obj)

            % make the interface

            zdlg = ZmapDialog();
            zdlg.AddHeader('Reasenberg Declustering parameters','FontSize',12);
            zdlg.AddHeader('look-ahead times');
            zdlg.AddDurationEdit('taumin', '(min) for UNclustered events' ,obj.taumin, 'TauMin look ahead time for not clustered events',@days);
            zdlg.AddDurationEdit('taumax', '(max) for   clustered events', obj.taumax,  'TauMax maximum look ahead time for clustered events',@days);
            zdlg.AddHeader('');
            zdlg.AddEdit('P',       'Confidence Level',             obj.P,          'P1 Confidence level : observing the next event in the sequence');
            zdlg.AddEdit('xk',      'XK factor',                    obj.xk,         'XK factor used in xmeff');
            zdlg.AddEdit('xmeff',   'Effective min mag cutoff',     obj.xmeff,   'XMEFF "effective" lower magnitude cutoff for catalog, during clusters, it is xmeff^{xk*cmag1}');
            zdlg.AddEdit('rfact',   'Interation radius factor:',    obj.rfact,      'RFACT>/b>factor for interaction radius for dependent events');
            zdlg.AddEdit('err',     'Epicenter error',              obj.err,        'Epicenter error');
            zdlg.AddEdit('derr',    'Depth error',                  obj.derr,       'derrDepth error');
            %zdlg.AddHeader('Output')
            zdlg.AddCheckbox('replaceSequenceWithEquivMainshock','Replace clusters with equivalent event',...
                obj.replaceSequenceWithEquivMainshock, {},...
                'Will replace each set of cluster earthquakes with a single event of equivalent Magnitude');
            zdlg.AddEdit('clusterDetailsVariableName',      'Save Clusters to workspace as', ...
                obj.clusterDetailsVariableName, 'if empty, then nothing will be separately saved');
            zdlg.AddEdit('declusteredCatalogVariableName',  'Save Declustered catalog to workspace as', ...
                obj.declusteredCatalogVariableName,'if empty, then nothing will be separately saved');
            zdlg.AddCheckbox('memorizeOriginalCatalog',     'Memorize Original catalog after sucessful decluster:', ...
                obj.memorizeOriginalCatalog, {}, 'Memorize original catalog prior to declustering');

            zdlg.Create('Name', 'Reasenberg Declustering','WriteToObj',obj,'OkFcn',@obj.Calculate);

        end

        function Results = Calculate(obj)
            calcFn = @obj.declus;
            [obj.declusteredCatalog, misc] = calcFn(obj);
            if nargout == 1
                Results = obj.declusteredCatalog;
            end

            obj.CalcFinishedFcn();
        end

        %{
        function [clustnum, eqs_in_clust] = decluster_from_python(obj)
              obj.verbose = config.get('verbose', obj.verbose);
            % Get relevant parameters
            neq = catalog.get_number_events()  % Number of earthquakes

            min_lookahead_days = obj.taumin;
            max_lookahead_days = obj.taumax;

            % Get elapsed days
            elapsed = days_from_first_event(catalog);

            assert(all(elapsed(2:end) >= elapsed(1:end-1)), "catalog needs to be in ascending date order")

            % easy-access variables
            dmethod='p2p';
            switch dmethod
                case 'gc'
                    surf_pos = [catalog.latitude, catalog.longitude];
                    event_distance = event_gc_distance;
                case'p2p'
                    surf_pos = geodetic_to_ecef(catalog.latitude, catalog.longitude);  % assumes at surface
                    event_distance = event_p2p_distance;
                otherwise
                    except(ValueError("unknown configuration dmethod. it should be 'gc' or 'p2p'"))
            end

            mags = catalog.magnitude;
            deps = catalog.depth;

            if isempty(obj.err)
                horiz_error = catalog.data.get('horizError', 0);
            else
                horiz_error = obj.err;
            end
            if isempty(obj.derr)
                depth_error = catalog.data.get('depthError', 0);
            else
                depth_error = obj.derr;
            end
            % Pre-allocate cluster index vectors
            vcl = zeros(neq, 1);

            % set the interaction zones, in km
            % Reasenberg 1987 or alternate version: Wells & Coppersmith 1994 / Helmstetter (SRL) 2007
            zone_noclust, zone_clust = obj.get_zone_distances_per_mag(mags, obj.rfact,...
                obj.interaction_formulas.(obj.interaction_formula), obj.taumax)

            k = 0  % clusterindex

            % variable to store information whether earthquake is already clustered
            clusmaxmag = ones(neq,1) * -inf;
            clus_biggest_idx = zeros(neq,1);

            % for every earthquake in catalog, main loop
            for i = 0 : neq-1 % in range(neq - 1)
                my_mag = mags(i);

                % variable needed for distance and timediff
                my_cluster = vcl(i);
                not_classified = my_cluster == 0;

                % attach interaction time

                if not_classified
                    obj.debug_print(i, ' is not in a cluster')
                    % this event is not associated with a cluster, yet
                    look_ahead_days = min_lookahead_days;

                elseif my_mag >= clusmaxmag(my_cluster)
                    % note, if this is now the biggest, then the cluster range collapses into its radius
                    printf('%d is the biggest event of cluster M=%g\n', i, my_mag);
                    % this is the biggest event  in this cluster, so far (or equal to it).
                    clusmaxmag(my_cluster) = my_mag;
                    clus_biggest_idx(my_cluster) = i;
                    look_ahead_days = min_lookahead_days;
                else
                    printf('%d is already in cluster, but not biggest', i);
                    % this event is already tied to a cluster, but is not the largest
                    idx_biggest = clus_biggest_idx(my_cluster);
                    days_since_biggest = elapsed(i) - elapsed(idx_biggest);
                    look_ahead_days = obj.clust_look_ahead_time(clusmaxmag(my_cluster),...
                        days_since_biggest, obj.xk, obj.xmeff, obj.P);

                    look_ahead_days(look_ahead_days<min_lookahead_days) = min_lookahead_days;
                    look_ahead_days(look_ahead_days>max_lookahead_days) = max_lookahead_days;
                end
                % extract eqs that fit interaction time window --------------

                max_elapsed = elapsed(i) + look_ahead_days;
                next_event = i + 1;
                last_event = bisect_left(elapsed, max_elapsed, next_event);
                temporal_evs = np.arange(next_event, last_event)
                if my_cluster ~= 0
                    temporal_evs = temporal_evs(vcl(temporal_evs) ~= my_cluster);
                end
                if len(temporal_evs) == 0
                    continue
                end
                % ------------------------------------
                % one or more events have now passed the time window test. Now compare
                % this subcatalog in space to A) most recent and B) largest event in cluster
                % ------------------------------------

                obj.debug_print('temporal_evs:', temporal_evs)
                my_biggest_idx = clus_biggest_idx(my_cluster)
                if not_classified
                    bg_ev_for_dist = i;
                else
                    bg_ev_for_dist = my_biggest_idx;
                end

                obj.debug_print('bg_ev_for_dist:', bg_ev_for_dist)
                % noinspection PyTypeChecker
                dist_to_recent = event_distance(surf_pos, deps, i, temporal_evs, horiz_error, depth_error)
                dist_to_biggest = event_distance(surf_pos, deps, bg_ev_for_dist, temporal_evs, horiz_error, depth_error)
                printf('dist_to_recent', dist_to_recent)
                printf('dist_to_biggest', dist_to_biggest)
                % extract eqs that fit the spatial interaction
                if look_ahead_days == min_lookahead_days
                    l_big = dist_to_biggest == 0;  % all false
                    l_recent = dist_to_recent <= zone_noclust(my_mag);
                    printf('Connecting those near to this event [dist <= {zone_noclust[my_mag]}]')
                else
                    l_big = dist_to_biggest <= zone_clust(clusmaxmag(my_cluster))
                    l_recent = dist_to_recent <= zone_clust(clusmaxmag(my_cluster))
                    printf('Connecting those near to this OR largest event [dist <= {zone_clust[clusmaxmag(my_cluster)]}]')
                end
                spatial_evs = l_recent | l_big;

                if ~any(spatial_evs)
                    continue
                end
                % ------------------------------------
                % one or more events have now passed both spatial and temporal window tests
                %
                % if there are events in this cluster that are already related to another
                % cluster, figure out the smallest cluster number. Then, assign all events
                % (both previously clustered and unclustered) to this new cluster number.
                % ------------------------------------

                % spatial events only include events AFTER i, not i itself
                % so vcl(events_in_any_cluster) is independent from vcl(i)

                candidates = temporal_evs(spatial_evs) ; % eqs that fit spatial and temporal criterion
                events_in_any_cluster = candidates(vcl(candidates) ~= 0);  % eqs which are already related with a cluster
                events_in_no_cluster = candidates(vcl(candidates) == 0);  % eqs that are not already in a cluster

                % if this cluster overlaps with any other cluster, then merge them
                % assign every event in all related clusters to the same (lowest) cluster number
                % set this cluster's maximum magnitude "clusmaxmag" to the largest magnitude of all combined events
                % set this cluster's clus_biggest_idx to the most recent largest event of all combined events

                if len(events_in_any_cluster) > 0
                    if not_classified
                        related_clust_nums = unique(vcl(events_in_any_cluster));
                    else
                        % include this cluster number in the reckoning
                        related_clust_nums = unique(np.hstack((vcl(events_in_any_cluster), my_cluster,)))
                    end
                    % associate all related events with my cluster
                    my_cluster = related_clust_nums(0);
                    vcl(i) = my_cluster;
                    vcl(candidates) = my_cluster;

                    for clustnum  = related_clust_nums
                        vcl(vcl == clustnum) = my_cluster;
                    end
                    events_in_my_cluster = vcl == my_cluster;
                    biggest_mag = np.max(mags(events_in_my_cluster));
                    biggest_mag_idx = find(mags == biggest_mag & events_in_my_cluster, 1, 'last');

                    % reset values for other clusters
                    clusmaxmag(related_clust_nums) = -inf;
                    clus_biggest_idx(related_clust_nums) = 0;

                    % assign values for this cluster
                    clusmaxmag(my_cluster) = biggest_mag;
                    clus_biggest_idx(my_cluster) = biggest_mag_idx;

                elseif my_cluster == 0
                    k = k + 1;
                    vcl(i) = k;
                    my_cluster = k;
                    clusmaxmag(my_cluster) = my_mag;
                    clus_biggest_idx(my_cluster) = i;
                else
                    pass  % no events found, and attached to existing cluster
                end
                % attach clustnumber to catalog yet unrelated to a cluster
                vcl(events_in_no_cluster) = my_cluster;

            end
            clustnum = vcl;
            eqs_in_clust = vcl > 0;
        end
        %}

        function [outputcatalog, details] = declus(obj, vals)
            % DECLUS main decluster algorithm
            % A.Allmann
            % main decluster algorithm
            % modified version, uses two different circles for already related events
            % works on catalog
            % different clusters stored with respective numbers in clus
            % Program is based on Raesenberg paper JGR;Vol90;Pages5479-5495;06/10/85
            %
            %basic variables used in the program
            %
            % interactzone_main_km  interaction zone for not clustered events
            % interactzone_in_clust_km     interaction zone for clustered events
            % rtest  radius in which the program looks for clusters
            % tau    look ahead time
            % tdiff  time difference between jth event and biggest eq
            % mbg    index of earthquake with biggest magnitude in a cluster
            % k      index of the cluster
            % k1     working index for cluster
            %
            % modified by Celso Reyes, 2017

            max_mag_in_cluster=[];
            idx_biggest_event_in_cluster=[];

            %% calculate interaction_zone  (1 value per event)

            interactzone_main_km = obj.InteractFormula.Reasenberg1985(obj.RawCatalog.Magnitude); %interaction zone for mainshock
            interactzone_in_clust_km = obj.rfact * interactzone_main_km;      %interaction zone if included in a cluster

            tau_min = days(obj.taumin);
            tau_max = days(obj.taumax);

            %calculation of the eq-time relative to 1902
            eqtime = days( obj.RawCatalog.Date - min(obj.RawCatalog.Date) );

            %variable to store information whether earthquake is already clustered
            clus = zeros(1,obj.RawCatalog.Count);

            k = 0;                                %clusterindex

            wai = waitbar(0,' Please Wait ...  Declustering the catalog');
            set(wai,'NumberTitle','off','Name','Declustering Progress');
            drawnow
            declustering_start = tic;
            %for every earthquake in catalog, main loop
            confine_value = @(value, min_val, max_val) max( min( value, max_val), min_val);
            for i = 1: (obj.RawCatalog.Count-1)

                % "myXXX" refers to the XXX for this event

                my_mag = obj.RawCatalog.Magnitude(i);

                if rem(i,50)==0
                    waitbar(i/(obj.RawCatalog.Count-1));
                end

                % variable needed for distance and timediff
                my_cluster = clus(i);
                alreadyInCluster = my_cluster~=0;
                not_classified = my_cluster==0;
                assert(not_classified~=alreadyInCluster)

                % attach interaction time
                if not_classified
                    look_ahead_days = tau_min;
                else
                    if my_mag >= max_mag_in_cluster(my_cluster)
                        max_mag_in_cluster(my_cluster) = my_mag;
                        idx_biggest_event_in_cluster(my_cluster) = i;
                        look_ahead_days = tau_min;
                    else
                        bgdiff = eqtime(i) - eqtime(idx_biggest_event_in_cluster(my_cluster));
                        look_ahead_days = clustLookAheadTime(obj.xk, max_mag_in_cluster(my_cluster), obj.xmeff, bgdiff, obj.P);
                        look_ahead_days = confine_value(look_ahead_days, tau_min, tau_max);
                    end
                end

                %extract eqs that fit interation time window
                temporal_evs = timediff(i, look_ahead_days, clus, eqtime); %local version

                if isempty(temporal_evs)
                    continue;
                end

                % ---------------------------
                % only continue if events passed the time test
                % ---------------------------

                rtest1 = interactzone_in_clust_km(i);
                if look_ahead_days == obj.taumin
                    rtest2 = 0;
                else
                    rtest2 = interactzone_main_km(idx_biggest_event_in_cluster(my_cluster));
                end

                if alreadyInCluster                % if i is already related with a cluster
                    tm1 = clus(temporal_evs) ~= my_cluster;  % eqs with a clustnumber different than i
                    if any(tm1)
                        temporal_evs = temporal_evs(tm1);
                    end
                    bg_ev_for_dist = idx_biggest_event_in_cluster(my_cluster);
                else
                    bg_ev_for_dist = i;
                end

                %calculate distances from the epicenter of biggest and most recent eq
                [dist1,dist2]=distance2(i,bg_ev_for_dist,temporal_evs, obj.RawCatalog);

                %extract eqs that fit the spatial interaction time
                sl0 = dist1<= rtest1 | dist2<= rtest2;

                if ~any(sl0)
                    continue
                end

                % ----------
                % only continue if events passed the distance test
                % ----------

                ll = temporal_evs(sl0);            %eqs that fit spatial and temporal criterion
                lla = ll(clus(ll)~=0);   %eqs which are already related with a cluster
                llb = ll(clus(ll)==0);   %eqs that are not already in a cluster
                if ~isempty(lla)         %find smallest clustnumber in the case several
                    sl1 = min(clus(lla));     %numbers are possible
                    if alreadyInCluster
                        my_cluster = min([sl1, my_cluster]);
                    else
                        my_cluster = sl1;
                    end
                    if clus(i)==0
                        clus(i) = my_cluster;
                    end
                    % merge related clusters together into cluster with the smallest number
                    sl2 = lla(clus(lla) ~= my_cluster);
                    if clus(i) ~= my_cluster
                        clus(clus==clus(i)) = my_cluster;
                    end

                    for j1 = sl2
                        if clus(j1) ~= my_cluster
                            clus(clus==clus(i)) = my_cluster;
                        end
                    end
                end

                if my_cluster==0   %if there was neither an event in the interaction zone nor i, already related to cluster
                    k = k+1;                         %
                    my_cluster = k;
                    clus(i) = my_cluster;
                    max_mag_in_cluster(my_cluster) = my_mag;
                    idx_biggest_event_in_cluster(my_cluster) = i;
                end

                if size(llb)>0     % attach clustnumber to events yet unrelated to a cluster
                    clus(llb) = my_cluster;  %
                end

            end

            close(wai);
            msg.dbfprintf('Declustering complete. It took %g seconds\n',toc(declustering_start));

            %% this table contains all we need to know about the clusters. maybe.
            details = table;
            details.Properties.UserData = struct;
            for j = 1 : numel(obj.ParameterableProperties)
                details.Properties.UserData.(obj.ParameterableProperties(j)) = obj.(obj.ParameterableProperties(j));
            end

            details.Properties.Description  = 'Details for cluster, from reasenberg declustering';
            details.eventNumber             = (1:obj.RawCatalog.Count)';
            details.clusterNumber           = clus(:);
            details.clusterNumber(details.clusterNumber==0) = missing;
            details.isBiggest               = false(size(details.clusterNumber));
            details.isBiggest(idx_biggest_event_in_cluster) = true;

            details.Latitude                = obj.RawCatalog.Y;
            details.Properties.VariableUnits(width(details)) = {'degrees'};

            details.Longitude               = obj.RawCatalog.X;
            details.Properties.VariableUnits(width(details)) = {'degrees'};

            details.Depth                   = obj.RawCatalog.Z;
            details.Properties.VariableUnits(width(details)) = {'kilometers'};

            details.Magnitude               = obj.RawCatalog.Magnitude;

            details.MagnitudeType           = obj.RawCatalog.MagnitudeType;

            details.Date                    = obj.RawCatalog.Date;

            details.InteractionZoneIfMain   = interactzone_main_km;
            details.Properties.VariableUnits(width(details)) = {'kilometers'};

            details.InteractionZoneIfInClust = interactzone_in_clust_km;
            details.Properties.VariableUnits(width(details)) = {'kilometers'};

            clusterFreeCatalog = obj.RawCatalog.subset(ismissing(details.clusterNumber));
            %biggest_events_in_cluster = obj.RawCatalog.subset(details.isBiggest);

            outputcatalog = clusterFreeCatalog;

            if ~any(clus)
                return
            end

            %build a matrix clust that stored clusters
            [~, biggest_events_in_cluster, max_mag_in_cluster,~,~] = funBuildclu(obj.RawCatalog,idx_biggest_event_in_cluster,clus,max_mag_in_cluster);

            % replace cluster sequences with summary events
            if obj.replaceSequenceWithEquivMainshock
                equi = obj.equevent(details(~ismissing(details.clusterNumber),:));  % calculates equivalent events
                if isempty(equi)
                    disp('No clusters in the catalog with this input parameters');
                    return;
                end
                tmpcat = cat(clusterFreeCatalog, equi);  %new, unsorted catalog

            else
                tmpcat = cat(clusterFreeCatalog, biggest_events_in_cluster); % builds catalog with biggest events instead
                disp('Original mainshocks kept');
            end

            tmpcat.sort('Date');

            tmpcat.Name = string(tmpcat.Name) + " (declust)";

            % save clustering details to workspace
            if ~isempty(obj.clusterDetailsVariableName)
                assignin('base',matlab.lang.makeValidName(obj.clusterDetailsVariableName),details);
            end

            if obj.memorizeOriginalCatalog
                mm = MemorizedCatalogManager();
                mm.memorize(obj.RawCatalog,'predeclust')
            end

            ZG              = ZmapGlobal.Data;
            ZG.original     = obj.RawCatalog;    %save catalog in variable original
            ZG.cluscat      = ZG.original.subset(clus(clus~=0));

            % save declustered catalog to workspace
            if ~isempty(obj.declusteredCatalogVariableName)
                assignin('base', obj.declusteredCatalogVariableName, tmpcat);
            end

            st1 = sprintf([' The declustering found %d clusters of earthquakes, a total of %d'...
                ' events (out of %d). The map window now displays the declustered catalog containing %d events.'], ...
                biggest_events_in_cluster.Count, ZG.cluscat.Count, ZG.original.Count , ZG.primeCatalog.Count);

            msgbox(st1,'Declustering Information')

            watchoff
            outputcatalog = tmpcat;

            obj.Result(1).values.cluster_details  = details;

        end

        function plot(obj, varargin)
            f = figure('Name','Reasenberg Deslustering Results');
            ax = subplot(2,2,1);
            ZG = ZmapGlobal.Data;
            biggest = obj.Result.values.cluster_details(obj.Result.values.cluster_details.isBiggest,:);
            non_cluster = obj.Result.values.cluster_details(ismissing(obj.Result.values.cluster_details.clusterNumber),:);
            msf = str2func(ZG.MainEventOpts.MarkerSizeFcn);
            scatter3(biggest.Longitude,biggest.Latitude,biggest.Depth,msf(biggest.Magnitude));
            ax.ZDir='reverse';
            title(ax,'Biggest events in each cluster');
            hold on
            ax.XLabel.String = 'Longitude';
            ax.YLabel.String = 'Latitude';
            ax.ZLabel.String = 'Depth [km]';
            feats = findobj(allchild(findobj('Tag','mainmap_ax')),'-regexp','Tag','mainmap_.+');
            copyobj(feats,ax); %copy features

            ax = subplot(2,2,2);

            ax = subplot(2,1,2);

            isInClust = ~ismissing(obj.Result.values.cluster_details.clusterNumber);
            isNotBig = ~obj.Result.values.cluster_details.isBiggest;
            clust = obj.Result.values.cluster_details(isInClust&isNotBig, :);
            scatter(clust.Date,clust.Depth,[],[.8 .8 .8],'Marker','.','DisplayName','other events in each cluster');
            ax.YDir='reverse';
            hold on;
            scatter(biggest.Date,biggest.Depth,msf(biggest.Magnitude),categorical(biggest.clusterNumber),'DisplayName','primary events in each cluster');
            cb = colorbar;
            cb.Label.String = 'Cluster #';
            ax.YLabel.String = 'Depth [km]';
            ax.XLabel.String = 'Date';
            title(ax,'Clusters through time');
            legend('show')
        end
    end

    methods(Static)
        function tau = clust_look_ahead_time(mag_big, dt_big, xk, xmeff, P)
            % CLUSTLOOKAHEAD calculate look ahead time for clustered events (days)
            deltam = (1-xk) .* mag_big - xmeff;
            if deltam < 0
                deltam = 0;
            end
            denom = 10.0 .^ ((deltam - 1) * (2/3));
            top = log(1 - P) * dt_big;
            tau = top / denom;
        end

        function h = AddMenuItem(parent, catalog, varargin)
            % create a menu item
            label = 'Reasenberg Decluster';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)ReasenbergDeclusterClass(catalog),...
                varargin{:});
        end

        function equi = equevent(tb)
            % equevent calc equivalent event to cluster
            % equi = equevent(catalog, cluster, bg)
            %   catalog : earthquake catalog
            %   cluster :
            %   bg : index of a big event (?)
            %  equevent.m                        A.Allmann
            % calculates equivalent event to a cluster
            % weight according to seismic moment
            % time for equivalent event is time of first biggest event
            %
            report_this_filefun();

            equi = ZmapCatalog;
            equi.Name='clusters';

            if isempty(tb)
                return
            end
            j = 0;
            nClusts = max(tb.clusterNumber);
            [elat, elon, edep, emag] = deal(nan(nClusts,1));
            edate(nClusts,1) = datetime(missing);
            emagtype(nClusts,1) = categorical(missing);

            for n = 1 : max(tb.clusterNumber)
                clust_events = tb(tb.clusterNumber==n,:);
                if isempty(clust_events)
                    continue;
                end
                j = j + 1;

                eqmoment = 10.^(clust_events.Magnitude .* 1.2);
                emoment = sum(eqmoment);         %moment

                weights = eqmoment./emoment;      %weightfactor
                elat(j)     = sum(clust_events.Latitude .* weights);
                elon(j)     = sum(clust_events.Longitude .* weights); %longitude
                edep(j)     = sum(clust_events.Depth .* weights); %depth
                emag(j)     = (log10(emoment))/1.2;
                theBiggest  = find(clust_events.isBiggest,1,'first');
                edate(j)    = clust_events.Date(theBiggest);
                emagtype(j) = clust_events.MagnitudeType(theBiggest);

            end

            %equivalent events for each cluster
            equi.Latitude = elat(1:j);
            equi.Longitude = elon(1:j);
            equi.Date = edate(1:j);
            equi.Magnitude = emag(1:j);
            equi.MagnitudeType = emagtype(1:j);
            equi.Depth = edep(1:j);
            [equi.Dip, equi.DipDirection, equi.Rake]=deal(nan(size(equi.Date)));
        end
    end

end

%% helper
function ac  = timediff(clus_idx, look_ahead_time, clus, eqtimes)
    % TIMEDIFF calculates the time difference between the ith and jth event
    % works with variable eqtime from function CLUSTIME
    % gives the indices ac of the eqs not already related to cluster k1
    % eqtimes should be sorted!
    %
    % clus_idx : ith cluster (ci)
    % look_ahead : look-ahead time (tau)
    % clus: clusters (length of catalog)
    % eqtimes: datetimes for event catalog, in days  [did not use duration because of overhead]
    %
    % tdiff: is time between jth event and eqtimes(clus_idx)
    % ac: index of each event within the cluster

    %assert(clus_idx <100, 'testing. remove me')

    comparetime = eqtimes(clus_idx);

    first_event = clus_idx + 1; % in cluster
    last_event = numel(eqtimes);
    max_elapsed = comparetime + look_ahead_time;

    if eqtimes(end) >= max_elapsed
        last_event = find(eqtimes(first_event : last_event) < max_elapsed, 1, 'last') + clus_idx;
    end

    if first_event == last_event
        % no additional events were found.
        ac = [];
        return
    end

    this_clusternum = clus(clus_idx);

    range = first_event : last_event;

    if this_clusternum == 0
        ac = range;
    else
        % indices of eqs not already related to this cluster
        ac = (find(clus(range) ~= this_clusternum)) + clus_idx;
    end
    ac = ac(:);
end

author: CelsoReyes


text:
function gf = gradddivHexp(Hkt_r, varargin)
% function gf = gradddivH(Hkt_r, varargin)
% gf is 1 x k*t row vector
% Hkt_r = reshape(Hkt,1,k*t) -> row vector
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %W matrix
% Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
% Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

alphaH=1; %for now....
Vxt = varargin{1};  %data
Wxk_tmp = varargin{2};  %W matrix
Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)
peval = varargin{5}; %parameters

if ~isfield(peval, 'w_lambda') peval.w_lambda=0; end

t=size(Vxt,2);
k=length(peval.h_dovec);

Hkt_tmp = exp(reshape(Hkt_r,k,t));

Wxk = zeros(peval.numpix, peval.ncomp+1);
Hkt = zeros(peval.ncomp+1, peval.nt);

Wxk(:,peval.w_dovec)=Wxk_tmp;
Hkt(peval.h_dovec,:)=Hkt_tmp;

Wxk(:,peval.w_fixvec)=Wxk_fix;
Hkt(peval.h_fixvec,:)=Hkt_fix;

deltasum=sum(sum(Wxk_tmp,1))-k;
if and(~isempty(Wxk_tmp), abs(deltasum)>10^-6)
    error('Wxk is not correctly normalized! (sum(Wxk_tmp,1)<>1)\n sum(Wxk_tmp,1)=%f',deltasum)
end

gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt; %d/dh(d-divergence)
% one is tehre because Wxt is normalized: sum(Wxt,1)=1
gf=reshape(gfkt(1:k,:),1,k*t); %making row vector

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


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

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

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

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

sigma = X' * X / m;
[U,S,V] = svd(sigma);

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

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
function [h_ctx,h_iskl,h_oskl,h_slp] = spm_eeg_inv_checkmeshes(varargin)
% Display tesselated surfaces of cortex, skull and scalp
% FORMAT [h_ctx,h_iskl,h_oskl,h_slp] = spm_eeg_inv_checkmeshes(S)
% S      - input data struct (optional)
%
% h_ctx  - handle to cortex patch
% h_iskl - handle to inner skull patch
% h_oskl - handle to outer skull patch
% h_slp  - handle to scalp patch
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Jeremie Mattout
% $Id: spm_eeg_inv_checkmeshes.m 6182 2014-09-18 12:03:18Z guillaume $

%-Initialisation
%--------------------------------------------------------------------------
[D,val] = spm_eeg_inv_check(varargin{:});
try
    %disp(D.inv{val}.mesh);
    mesh  = spm_eeg_inv_transform_mesh(eye(4), D.inv{val}.mesh);
    Mctx  = mesh.tess_ctx;
    Miskl = mesh.tess_iskull;
    Moskl = mesh.tess_oskull;
    Mslp  = mesh.tess_scalp;
catch
    spm('alert!','Please create meshes',mfilename);
    return
end

if spm('CmdLine'), [h_ctx,h_iskl,h_oskl,h_slp] = deal([]); return; end

%-SPM Graphics figure
%--------------------------------------------------------------------------
Fgraph  = spm_figure('GetWin','Graphics'); spm_figure('Clear',Fgraph)

%-Cortex Mesh Display
%--------------------------------------------------------------------------
h_ctx   = patch('vertices',Mctx.vert,'faces',Mctx.face,...
            'EdgeColor','b','FaceColor','b');
hold on

%-Inner-skull Mesh Display
%--------------------------------------------------------------------------
h_iskl  = patch('vertices',Miskl.vert,'faces',Miskl.face,...
            'EdgeColor','r','FaceColor','none');

%-Outer-skull Mesh Display
%--------------------------------------------------------------------------
h_oskl  = patch('vertices',Moskl.vert,'faces',Moskl.face,...
            'EdgeColor',[1 .5 .35],'FaceColor','none');

%-Inner Scalp Mesh Display
%--------------------------------------------------------------------------
h_slp   = patch('vertices',Mslp.vert,'faces',Mslp.face,...
            'EdgeColor',[1 .7 .55],'FaceColor','none');

%-Slices
%--------------------------------------------------------------------------
pls = 0.05:0.2:0.9;
N   = nifti(mesh.sMRI);
d   = size(N.dat);
pls = round(pls.*d(3));
for i=1:numel(pls)
    [x,y,z] = ndgrid(1:d(1),1:d(2),pls(i));
    f1 = N.dat(:,:,pls(i));
    M  = N.mat;
    x1 = M(1,1)*x+M(1,2)*y+M(1,3)*z+M(1,4);
    y1 = M(2,1)*x+M(2,2)*y+M(2,3)*z+M(2,4);
    z1 = M(3,1)*x+M(3,2)*y+M(3,3)*z+M(3,4);

    s  = surf(x1,y1,z1,f1);
    set(s,'EdgeColor','none')
end

%-Fiducials
%--------------------------------------------------------------------------
pnt   = mesh.fid.fid.pnt;
h_fid = plot3(pnt(:,1), pnt(:,2), pnt(:,3), '.c', 'MarkerSize', 30);
h_fid_txt = text(pnt(:,1), pnt(:,2), pnt(:,3), mesh.fid.fid.label);

axis image off;
colormap('gray');
view(-135,45);
rotate3d on
drawnow
hold off

author: spm


text:
function [i, j] = toindex(x, y)
global F;

j = ceil(x);
i = size(F,1) - floor(y);

author: Sable

repo: mcbench-benchmarks


text:
function [O, DzDw  ] = convo ( I , gamma , DzDy1 ,  DzDy2)
%% Convolution layer
%% This layer has a parameter: D_{l} = sum_{m=1}^{s} gamma_{l,m}B_{m};
%% Copyright (c) 2017 Yan Yang
%% All rights reserved.

%% network setting
config;
fN = nnconfig.FilterNumber;
fS = nnconfig.FilterSize;
pad = nnconfig.Padding;
gp = nnconfig.EnableGPU;
s=fS*fS-1;
[m,n] = size(I);

D = zeros(fS, fS, fN);
B = filter_base( );
for i = 1:fN
    D(:,:,i) = reshape(B*gamma(:,i),fS,fS);
end
DT = rot90(D,2);

if nargin == 2
    for i=1:fN
        O(:,:,i) = imfilter( double(I) ,double(D(:,:,i)),'same','circular','conv');
    end
end

if nargin == 4
    DzDy = DzDy1 + DzDy2 ;
    if gp
        Dp = gpuArray(D);
        E = gpuArray(zeros(m,m));
        D1p = gpuArray(DT);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  O
        for i = 1:fN
            d1 = padImage_g(DzDy(:,:,i),[pad,pad],'circular');
            d1 = gpuArray(d1);
            dd = conv2(d1,D1p(:,:,i),'valid');
            E = E + dd;
        end
        O = gather(E);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DzDw
        DzDy = gpuArray(DzDy);
        Bj1p = FBASE(B,1);
        Bjp = rot90(Bj1p,2);
        DzDw = gpuArray(zeros(s,fN));
        d1 = padImage_g(I,[pad,pad],'circular');
        d1 = gpuArray(d1);

        for j = 1:fN
            for k=1:s  % m
                dd = conv2(d1,Bjp(:,:,k),'valid');
                tp = DzDy(:,:,j).* dd;
                DzDw (k,j)=sum(tp(:));
            end
        end

        DzDw=gather(DzDw);

    else

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  O
        O = zeros(m,n);
        DT = rot90(D,2);
        for i=1:fN
            O = O + imfilter(double(DzDy(:,:,i)),double(DT(:,:,i)),'same','circular','conv');
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DzDw
        for j = 1:fN
            for k=1:s
                Bj = reshape(B(:,k),fS,fS);
                tp = DzDy(:,:,j).* imfilter(double(I),double(Bj),'same','circular','conv');
                DzDw (k,j)=sum(tp(:));
            end
        end
    end
end
end

repo: Deep-ADMM-Net


text:
% PLOTPROJ - plot projections of one or more ICA components along with
%              the original data (returns the data plotted)
%
% Usage:
%   >> [projdata] = plotproj(data,weights,compnums);
%   >> [projdata] = plotproj(data,weights,compnums, ...
%                                 title,limits,chanlist,channames,colors);
%
% Inputs:
%   data        = single epoch of RUNICA input data (chans,frames)
%   weights     = unmixing matrix (=weights*sphere)
%   compnums    = vector of component numbers to project and plot
%
% Optional inputs:
%   title       = 'fairly short plot title' {0 -> 'PLOTPROJ'}
%   limits      = [xmin xmax ymin ymax]  (x's in msec)
%                          {0, or both y's 0 -> data limits}
%   chanlist    = list of data channels to plot {0 -> all}
%   channames   = channel location file or structure (see READLOCS)
%   colors      = file of color codes, 3 chars per line  ('.' = space)
%                          {0 -> default color order (black/white first)}
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 05-01-96
%
% See also: PLOTDATA

% Without color arg, reads filename for PROJCOLORS from icadefs.m

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

% 05-25-96 added chanlist, nargin tests, rearranged variable order -sm
% 07-29-96 debugged, added chanlist channames option -sm
% 10-26-96 added test for column of compnums  -sm
% 02-18-97 improved usage message -sm
% 02-19-97 merged versions -sm
% 03-19-97 changed VAR to diag(COV), use datamean arg instead of frames/baseframes -sm
% 04-24-97 tested datamean for 1-epoch; replaced COV with mean-squares -sm
% 05-20-97 read icadefs for PROJCOLORS & MAXPLOTDATACHANS -sm
% 06-05-97 use arbitrary chanlist as default channames -sm
% 06-07-97 changed order of args to conform to runica -sm
% 06-12-97 made sumdata(chanlist in line 159 below -sm
% 07-23-97 removed datamean from args; let mean distribut4e over components -sm
% 09-09-97 corrected write out line " summing " -sm
% 10-31-97 removed errcode var -sm
% 11-05-97 added test for channames when chanlist ~= 1:length(chanlist) -sm & ch
% 12-19-00 adjusted new ICAPROJ args -sm
% 01-12-01 removed sphere arg -sm
% 01-25-02 reformated help & license, added links -ad

function [projdata] = plotproj(data,weights,compnums,titl,limits,chanlist,channels,colors);

icadefs       % read default PROJCOLORS & MAXPLOTDATACHANS variables from icadefs.m
DEFAULT_TITLE = 'plotproj()';

%
% Substitute for missing arguments
%

if nargin < 8,
    colors = 'white1st.col';
elseif colors==0,
    colors = 'white1st.col';
end

if nargin < 7,
    channels = 0;
end
if nargin < 6
    chanlist = 0;
end
if nargin < 5,
    limits = 0;
end
if nargin < 4,
    titl = 0;
end
if titl==0,
    titl = DEFAULT_TITLE;
end

if nargin < 3,
    fprintf('plotproj(): must have at least four arguments.\n\n');
    help plotproj
    return
end
%
% Test data size
%
[chans,framestot] = size(data);

frames = framestot; % assume one epoch

[wr,wc]           = size(weights);

if wc ~= chans
    fprintf('plotproj(): sizes of weights and data incompatible.\n\n');
    return
end
%
% Substitute for 0 arguments
%
if chanlist == 0,
    chanlist = [1:chans];
end
if compnums == 0,
    compnums = [1:wr];
end
if size(compnums,1)>1,        % handle column of compnums !
    compnums = compnums';
end
if length(compnums) > 256,
    fprintf('plotproj(): cannot plot more than %d channels of data at once.\n',256);
    return
end

if channels ~= 0  % if chan name file given
   if ~all(chanlist == [1:length(chanlist)])
     fprintf('plotproj(): Cannot read an arbitrary chanlist of channel names.\n');
     return
   end
end
if channels==0,
    channels = chanlist;
end

if max(compnums)>wr,
   fprintf(...
 '\n    plotproj(): Component index (%d) > number of components (%d).\n', ...
                                 max(compnums),wr);
   return
end

fprintf('Reconstructing (%d chan, %d frame) data summing %d components.\n', ...
                   chans,frames,length(compnums));
%
% Compute projected data for single components
%
projdata = data(chanlist,:);
fprintf('plotproj(): Projecting component(s) ');
for s=compnums,            % for each component
   fprintf('%d ',s);
   proj = icaproj(data,weights,s);   % let offsets distribute
   projdata = [projdata proj(chanlist,:)];  % append projected data onto projdata
   % size(projdata) = [length(chanlist)  framestot*(length(compnums)+1)]
end
fprintf('\n');
%
% Compute percentage of variance accounted for
%
sumdata = icaproj(data,weights,compnums);% let offsets distribute
sigmssq = mean(sum(data(chanlist,:).*data(chanlist,:)));
   data(chanlist,:) = data(chanlist,:) - sumdata(chanlist,:);
difmssq = mean(sum(data(chanlist,:).*data(chanlist,:)));

pvaf   = round(100.0*(1.0-difmssq/sigmssq)); % percent variance accounted for

rtitl = ['(p.v.a.f. ' int2str(pvaf) '%)'];
%
% Make the plot
%
plotdata(projdata,length(data),limits,titl,channels,colors,rtitl);
                                                % make the plot

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Test file for Scaling Functions
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef test_scales < matlab.unittest.TestCase

    properties
        scaling_function
    end

    methods(TestMethodSetup)
        function setScalingFunction(testCase)
            testCase.scaling_function = MelScaling();
        end
    end

    methods (Test)
        function test_scales_invertible(testCase)
             for hertz = 200:100:2000
                scale = testCase.scaling_function.hertz_to_scale(hertz);
                testCase.verifyEqual(hertz, testCase.scaling_function.scale_to_hertz(scale),'relTol',sqrt(eps));
             end
        end
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [np_info] = np_readfileinfo (filename, option)
%
% function [np_info] = np_readfileinfo (filename, option)
%
% Purpose:
%
% np_readfileinfo reads out header information from a NEURO PRAX
% data file (*.EEG).
%
% Syntax:
%
%   [np_info] = np_readfileinfo(filename);
%   [np_info] = np_readfileinfo(filename,'NO_MINMAX');
%
% Input data:
%
%   filename        -   the complete filename with path
%                       (e. g.  C:\Document...\20030716103637.EEG)
%   option          -   if option = 'NO_MINMAX' then physical minima
%   (optional)          and maxima off all channels will not be calculated
%                       (faster for long recordings)
%
% Output data:
%
%   np_info         -   structure
%
%      np_info.filename       -   filename of *.EEG file
%      np_info.pathname       -   pathname of *.EEG file
%      np_info.name           -   patient's name
%      np_info.firstname      -   patient's firstname
%      np_info.birthday       -   patient's birthday
%      np_info.ID             -   identification number
%      np_info.date           -   start date of the recording (format: 'dd-mmm-yyyy')
%      np_info.time           -   start time of the recording (format: 'hh:mm:ss')
%      np_info.duration       -   duration of the recording (global for splitted EEG files!)
%      np_info.setup          -   electrode setup
%      np_info.pmtype         -   type of primary montage
%      np_info.algorithm      -   used algorithm for measurement
%      np_info.channels       -   cell-array with channel names
%      np_info.channeltypes   -   cell-array with channel types
%      np_info.units          -   cell-array with channel units
%      np_info.PhysMin        -   physical minimum for each channel (global for splitted EEG files!)
%      np_info.PhysMax        -   physical maximum for each channel (global for splitted EEG files!)
%      np_info.N              -   number of samples per channel (global for splitted EEG files!)
%      np_info.K              -   number of channels
%      np_info.fa             -   sampling frequency
%      np_info.fp_data        -   filepointer to data
% --- additional SPLITTING information ---
%      np_info.SPLIT_Z        -   number of splitted EEG files (=Z; Z = 1 : no splitting)
%      np_info.SPLIT_filename -   all filenames of splitted EEG file (Zx1 cell-array)
%      np_info.SPLIT_N        -   samples of splitted EEG file (Zx1 array)
%      np_info.SPLIT_fp_data  -   file pointers of splitted EEG file (Zx1 array)
%      np_info.SPLIT_PhysMin  -   physical minimum for each channel and each file (ZxK array)
%      np_info.SPLIT_PhysMax  -   physical maximum for each channel and each file (ZxK array)
%
% Version:  1.3. (2005-09-19)
%
% (1) The field 'units' will be read correctly from the channel header.
%
% (2) The sampling frequency (fa) is identical for all channels and is set
%     to the value of the first channel.
%
% (3) The channel types will be read from the channel header directly.
%
% (4) The field 'time' is set to the correct recording time.
%
% (5) Additional structure fields:
%     setup      -   the electrode setup (previously "primmon")
%     pmtype     -   the type of the primary montage
%     algorithm  -   the feedback algorithm
%
% (6) No longer available: the structure field "primmon".
%
% (7) Splitted EEG files will be supported.
%
% See also: np_readdata, np_readmarker
%
% eldith GmbH
% Gustav-Kirchhoff-Str. 5
% D-98693 Ilmenau
% Germany
% 02.02.2005

if (1==nargin),
    option='';
end

% -------------------------------------------------------------------------
% Initialisierung der Struktur
% -------------------------------------------------------------------------
np_info=struct('filename','','pathname','','name','','firstname','','birthday','','ID',0,'date','00.00.00',...
               'time','00:00:00','duration','','setup','','pmtype','','algorithm','','channels','','channeltypes','',...
               'units','','PhysMin',0,'PhysMax',0,'N',0,'K',0,'fa',0,'fp_data',0,...
               'SPLIT_Z',1,'SPLIT_filename','','SPLIT_N',[],'SPLIT_fp_data',[],...
               'SPLIT_PhysMin',[],'SPLIT_PhysMax',[]);

% -------------------------------------------------------------------------
% alle zusammengeh�renden EEG-Dateien ermitteln (z.B. bei SPLITTING)
% -------------------------------------------------------------------------
[pa,fn,ex]=fileparts(filename);
if strcmp(pa,''),
    pa=pwd;
end
MaskEEGFileNames=[pa filesep fn(1:14) '*.EEG'];
d=dir(MaskEEGFileNames);
e=cell(length(d),1);
for i=1:length(d)
    e{i}=[pa filesep d(i).name];
end
filenames=sort(e);

% -------------------------------------------------------------------------
% SPLIT-Procedures
% -------------------------------------------------------------------------
np_info.SPLIT_Z=length(filenames);
np_info.SPLIT_filename=cell(np_info.SPLIT_Z,1);
np_info.SPLIT_N=zeros(np_info.SPLIT_Z,1);
np_info.SPLIT_fp_data=zeros(np_info.SPLIT_Z,1);

for Z=1:np_info.SPLIT_Z
    % Infos f�r n�chste EEG-Datei lesen
    np_splitted_info=np_read_splitted_fileinfo(filenames{Z},option);

    % fixe Daten in np_info eintragen
    if (1==Z),
        np_info.filename=np_splitted_info.filename;
        np_info.pathname=np_splitted_info.pathname;
        np_info.name=np_splitted_info.name;
        np_info.firstname=np_splitted_info.firstname;
        np_info.birthday=np_splitted_info.birthday;
        np_info.ID=np_splitted_info.ID;
        np_info.date=np_splitted_info.date;
        np_info.time=np_splitted_info.time;
        np_info.setup=np_splitted_info.setup;
        np_info.pmtype=np_splitted_info.pmtype;
        np_info.algorithm=np_splitted_info.algorithm;
        np_info.K=np_splitted_info.K;
        np_info.fa=np_splitted_info.fa;
        np_info.fp_data=np_splitted_info.fp_data;
        np_info.channels=np_splitted_info.channels;
        np_info.channeltypes=np_splitted_info.channeltypes;
        np_info.units=np_splitted_info.units;
        np_info.PhysMin=zeros(1,np_splitted_info.K);
        np_info.PhysMax=zeros(1,np_splitted_info.K);
        np_info.SPLIT_PhysMin=zeros(np_info.SPLIT_Z,np_splitted_info.K);
        np_info.SPLIT_PhysMax=zeros(np_info.SPLIT_Z,np_splitted_info.K);
    end

    % SPLIT-Infos zusammenstellen
    np_info.SPLIT_filename{Z}=np_splitted_info.filename;
    np_info.SPLIT_N(Z)=np_splitted_info.N;
    np_info.SPLIT_fp_data(Z)=np_splitted_info.fp_data;
    np_info.SPLIT_PhysMin(Z,:)=np_splitted_info.PhysMin;
    np_info.SPLIT_PhysMax(Z,:)=np_splitted_info.PhysMax;

    % Samplez�hler erh�hen
    np_info.N=np_info.N+np_info.SPLIT_N(Z);

end

% global duration berechnen
np_info.fa;
np_info.duration=num2str(np_info.N./np_info.fa);

% physikalische Minima und Maxima berechnen
if (1~=strcmp(upper(option),'NO_MINMAX')),
    for k=1:np_info.K
        np_info.PhysMin(k)=+1e38;
        np_info.PhysMax(k)=-1e38;
        for Z=1:np_info.SPLIT_Z
            np_info.PhysMin(k)=min([np_info.PhysMin(k) np_info.SPLIT_PhysMin(Z,k)]);
            np_info.PhysMax(k)=max([np_info.PhysMax(k) np_info.SPLIT_PhysMax(Z,k)]);
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ c, s ] = spherical_harmonic ( l, m, theta, phi )

%*****************************************************************************80
%
%% SPHERICAL_HARMONIC evaluates spherical harmonic functions.
%
%  Discussion:
%
%    The spherical harmonic function Y(L,M,THETA,PHI)(X) is the
%    angular part of the solution to Laplace's equation in spherical
%    coordinates.
%
%    Y(L,M,THETA,PHI)(X) is related to the associated Legendre
%    function as follows:
%
%      Y(L,M,THETA,PHI)(X) = FACTOR * P(L,M)(cos(THETA)) * exp ( i * M * PHI )
%
%    Here, FACTOR is a normalization factor:
%
%      FACTOR = sqrt ( ( 2 * L + 1 ) * ( L - M )! / ( 4 * PI * ( L + M )! ) )
%
%    In Mathematica, a spherical harmonic function can be evaluated by
%
%      SphericalHarmonicY [ l, m, theta, phi ]
%
%    Note that notational tradition in physics requires that THETA
%    and PHI represent the reverse of what they would normally mean
%    in mathematical notation; that is, THETA goes up and down, and
%    PHI goes around.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Eric Weisstein,
%    CRC Concise Encyclopedia of Mathematics,
%    CRC Press, 1999.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input, integer L, the first index of the spherical harmonic function.
%    Normally, 0 <= L.
%
%    Input, integer M, the second index of the spherical harmonic function.
%    Normally, -L <= M <= L.
%
%    Input, real THETA, the polar angle, for which
%    0 <= THETA <= PI.
%
%    Input, real PHI, the longitudinal angle, for which
%    0 <= PHI <= 2*PI.
%
%    Output, real C(1:L+1), S(1:L+1), the real and imaginary
%    parts of the functions Y(L,0:L,THETA,PHI).
%
  m_abs = abs ( m );

  plm(1:l+1) = legendre_associated_normalized ( l, m_abs, cos ( theta ) );

  c(1:l+1) = plm(1:l+1) * cos ( m * phi );
  s(1:l+1) = plm(1:l+1) * sin ( m * phi );

  if ( m < 0 )
    c(1:l+1) = -c(1:l+1);
    s(1:l+1) = -s(1:l+1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% testing MDP
function metrics = MDP_test_hungarian(seq_idx, seq_set, tracker)

is_show = 0;
is_save = 1;
is_text = 0;
is_pause = 0;

opt = globals();
opt.is_text = is_text;
opt.exit_threshold = 0.7;

if is_show
    close all;
end

if strcmp(seq_set, 'train') == 1
    seq_name = opt.mot2d_train_seqs{seq_idx};
    seq_num = opt.mot2d_train_nums(seq_idx);
else
    seq_name = opt.mot2d_test_seqs{seq_idx};
    seq_num = opt.mot2d_test_nums(seq_idx);
end

% build the dres structure for images
filename = sprintf('%s/%s_dres_image.mat', opt.results, seq_name);
if exist(filename, 'file') ~= 0
    object = load(filename);
    dres_image = object.dres_image;
    fprintf('load images from file %s done\n', filename);
else
    dres_image = read_dres_image(opt, seq_set, seq_name, seq_num);
    fprintf('read images done\n');
    save(filename, 'dres_image', '-v7.3');
end

% read detections
filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'det', 'det.txt');
dres_det = read_mot2dres(filename);

if strcmp(seq_set, 'train') == 1
    % read ground truth
    filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'gt', 'gt.txt');
    dres_gt = read_mot2dres(filename);
    dres_gt = fix_groundtruth(seq_name, dres_gt);
end

% load the trained model
if nargin < 3
    object = load('tracker.mat');
    tracker = object.tracker;
end

% intialize tracker
I = dres_image.I{1};
tracker = MDP_initialize_test(tracker, size(I,2), size(I,1), dres_det, is_show);

% for each frame
trackers = [];
id = 0;
for fr = 1:seq_num
    if is_text
        fprintf('frame %d\n', fr);
    else
        fprintf('.');
        if mod(fr, 100) == 0
            fprintf('\n');
        end
    end

    % extract detection
    index = find(dres_det.fr == fr);
    dres = sub(dres_det, index);

    % nms
%     boxes = [dres.x dres.y dres.x+dres.w dres.y+dres.h dres.r];
%     index = nms_new(boxes, 0.6);
%     dres = sub(dres, index);

    dres = MDP_crop_image_box(dres, dres_image.Igray{fr}, tracker);

    if is_show
        figure(1);

        % show ground truth
        if strcmp(seq_set, 'train') == 1
            subplot(2, 2, 1);
            show_dres(fr, dres_image.I{fr}, 'GT', dres_gt);
        end

        % show detections
        subplot(2, 2, 2);
        show_dres(fr, dres_image.I{fr}, 'Detections', dres_det);
    end

    % separate trackers into the first and the second class
    [index1, index2] = sort_trackers(trackers);

    for k = 1:2
        % process trackers in the first class or the second class
        if k == 1
            index_track = index1;
        else
            index_track = index2;
        end

        num_track = numel(index_track);
        flags = zeros(num_track, 1);

        % process tracked targets
        for i = 1:num_track
            ind = index_track(i);
            if trackers{ind}.state == 0
                flags(i) = 1;
            elseif trackers{ind}.state == 2
                % track target
                trackers{ind} = track(fr, dres_image, dres, trackers{ind}, opt);
                % connect target
                if trackers{ind}.state == 3
                    if k == 1
                        index_tmp = index_track(1:i-1);
                    else
                        index_tmp = [index_track(1:i-1); index1];
                    end
                    [dres_tmp, index] = generate_initial_index(trackers(index_tmp), dres);
                    dres_associate = sub(dres_tmp, index);
                    trackers{ind} = associate(fr, dres_image,  dres_associate, trackers{ind}, opt);
                end

                if trackers{ind}.state == 2 || trackers{ind}.state == 0
                    flags(i) = 1;
                end
            end
        end

        % process lost targets
        if k == 1
            index_tmp = index_track(flags == 1);
        else
            index_tmp = [index_track(flags == 1); index1];
        end
        [dres_tmp, index] = generate_initial_index(trackers(index_tmp), dres);
        dres_associate = sub(dres_tmp, index);
        num_det = numel(index);

        % compute distance matrix
        index_track = index_track(flags == 0);
        num_track = numel(index_track);
        dist = zeros(num_track, num_det);
        for i = 1:num_track
            % lost target
            ind = index_track(i);
            dist(i,:) = compute_distance(fr, dres_image, dres_associate, trackers{ind});
        end

        % Hungarian algorithm
        assignment = assignmentoptimal(dist);

        % process the assignment
        for i = 1:numel(assignment)
            det_id = assignment(i);
            ind = index_track(i);
            if det_id == 0
                % no association
                trackers{ind}.state = 3;
                dres_one = sub(trackers{ind}.dres, numel(trackers{ind}.dres.fr));
                dres_one.fr = fr;
                dres_one.id = trackers{ind}.target_id;
                dres_one.state = 3;

                if trackers{ind}.dres.fr(end) == fr
                    dres_tmp = trackers{ind}.dres;
                    index_tmp = 1:numel(dres_tmp.fr)-1;
                    trackers{ind}.dres = sub(dres_tmp, index_tmp);
                end
                trackers{ind}.dres = concatenate_dres(trackers{ind}.dres, dres_one);
            else
                % association
                dres_one = sub(dres_associate, det_id);
                trackers{ind} = LK_associate(fr, dres_image, dres_one, trackers{ind});

                trackers{ind}.state = 2;
                % build the dres structure
                dres_one = [];
                dres_one.fr = fr;
                dres_one.id = trackers{ind}.target_id;
                dres_one.x = trackers{ind}.bb(1);
                dres_one.y = trackers{ind}.bb(2);
                dres_one.w = trackers{ind}.bb(3) - trackers{ind}.bb(1);
                dres_one.h = trackers{ind}.bb(4) - trackers{ind}.bb(2);
                dres_one.r = 1;
                dres_one.state = 2;

                if trackers{ind}.dres.fr(end) == fr
                    dres_tmp = trackers{ind}.dres;
                    index_tmp = 1:numel(dres_tmp.fr)-1;
                    trackers{ind}.dres = sub(dres_tmp, index_tmp);
                end
                trackers{ind}.dres = interpolate_dres(trackers{ind}.dres, dres_one);
                % update LK tracker
                trackers{ind} = LK_update(fr, trackers{ind}, dres_image.Igray{fr}, dres_associate, 1);
            end
        end
    end

    % find detections for initialization
    [dres, index] = generate_initial_index(trackers, dres);
    for i = 1:numel(index)
        % extract features
        dres_one = sub(dres, index(i));
        f = MDP_feature_active(tracker, dres_one);
        % prediction
        label = svmpredict(1, f, tracker.w_active, '-q');
        % make a decision
        if label < 0
            continue;
        end

        % reset tracker
        tracker.prev_state = 1;
        tracker.state = 1;
        id = id + 1;

        trackers{end+1} = initialize(fr, dres_image, id, dres, index(i), tracker);
    end

    % resolve tracker conflict
    trackers = resolve(trackers, dres, opt);

    dres_track = generate_results(trackers);
    if is_show
        figure(1);

        % show tracking results
        subplot(2, 2, 3);
        show_dres(fr, dres_image.I{fr}, 'Tracking', dres_track, 2);

        % show lost targets
        subplot(2, 2, 4);
        show_dres(fr, dres_image.I{fr}, 'Lost', dres_track, 3);

        if is_pause
            pause();
        else
            pause(0.01);
        end
    end
end

% write tracking results
filename = sprintf('%s/%s.txt', opt.results, seq_name);
fprintf('write results: %s\n', filename);
write_tracking_results(filename, dres_track, opt.tracked);

% evaluation
if strcmp(seq_set, 'train') == 1
    benchmark_dir = fullfile(opt.mot, opt.mot2d, seq_set, filesep);
    metrics = evaluateTracking({seq_name}, opt.results, benchmark_dir);
else
    metrics = [];
end

% save results
if is_save
    filename = sprintf('%s/%s_results.mat', opt.results, seq_name);
    save(filename, 'dres_track', 'metrics');
end

% sort trackers according to number of tracked frames
function [index1, index2] = sort_trackers(trackers)

sep = 10;
num = numel(trackers);
len = zeros(num, 1);
state = zeros(num, 1);
for i = 1:num
    len(i) = trackers{i}.streak_tracked;
    state(i) = trackers{i}.state;
end

index1 = find(len > sep);
[~, ind] = sort(state(index1));
index1 = index1(ind);

index2 = find(len <= sep);
[~, ind] = sort(state(index2));
index2 = index2(ind);

% initialize a tracker
% dres: detections
function tracker = initialize(fr, dres_image, id, dres, ind, tracker)

if tracker.state ~= 1
    return;
else  % active

    % initialize the LK tracker
    tracker = LK_initialize(tracker, fr, id, dres, ind, dres_image);
    tracker.state = 2;
    tracker.streak_occluded = 0;
    tracker.streak_tracked = 0;

    % build the dres structure
    dres_one.fr = dres.fr(ind);
    dres_one.id = tracker.target_id;
    dres_one.x = dres.x(ind);
    dres_one.y = dres.y(ind);
    dres_one.w = dres.w(ind);
    dres_one.h = dres.h(ind);
    dres_one.r = dres.r(ind);
    dres_one.state = tracker.state;
    tracker.dres = dres_one;
end

% track a target
function tracker = track(fr, dres_image, dres, tracker, opt)

% tracked
if tracker.state == 2
    tracker.streak_occluded = 0;
    tracker.streak_tracked = tracker.streak_tracked + 1;
    tracker = MDP_value(tracker, fr, dres_image, dres, []);

    % check if target outside image
    [~, ov] = calc_overlap(tracker.dres, numel(tracker.dres.fr), dres_image, fr);
    if ov < opt.exit_threshold
        if opt.is_text
            fprintf('target outside image by checking boarders\n');
        end
        tracker.state = 0;
    end
end

% associate a lost target
function tracker = associate(fr, dres_image, dres_associate, tracker, opt)

% occluded
if tracker.state == 3
    tracker.streak_occluded = tracker.streak_occluded + 1;
    % find a set of detections for association
    [dres_associate, index_det] = generate_association_index(tracker, fr, dres_associate);
    tracker = MDP_value(tracker, fr, dres_image, dres_associate, index_det);
    if tracker.state == 2
        tracker.streak_occluded = 0;
    end

    if tracker.streak_occluded > opt.max_occlusion
        tracker.state = 0;
        if opt.is_text
            fprintf('target %d exits due to long time occlusion\n', tracker.target_id);
        end
    end

    % check if target outside image
    [~, ov] = calc_overlap(tracker.dres, numel(tracker.dres.fr), dres_image, fr);
    if ov < opt.exit_threshold
        if opt.is_text
            fprintf('target outside image by checking boarders\n');
        end
        tracker.state = 0;
    end
end

% associate a lost target
function dist = compute_distance(fr, dres_image, dres_associate, tracker)

% occluded
if tracker.state == 3
    % find a set of detections for association
    [dres_associate, index_det] = generate_association_index(tracker, fr, dres_associate);

    % compute association scores
    num = numel(dres_associate.fr);
    dist = Inf(1, num);
    if isempty(index_det) == 0
        % extract features with LK association
        dres = sub(dres_associate, index_det);
        [features, flag] = MDP_feature_occluded(fr, dres_image, dres, tracker);

        m = size(features, 1);
        labels = -1 * ones(m, 1);
        [~, ~, probs] = svmpredict(labels, features, tracker.w_occluded, '-b 1 -q');

        dist(index_det(flag == 1)) = probs(flag == 1, 2);
        dist(dist > 0.5) = inf;
    end
end

% resolve conflict between trackers
function trackers = resolve(trackers, dres_det, opt)

% collect dres from trackers
dres_track = [];
for i = 1:numel(trackers)
    tracker = trackers{i};
    dres = sub(tracker.dres, numel(tracker.dres.fr));

    if tracker.state == 2
        if isempty(dres_track)
            dres_track = dres;
        else
            dres_track = concatenate_dres(dres_track, dres);
        end
    end
end

% compute overlaps
num_det = numel(dres_det.fr);
if isempty(dres_track)
    num_track = 0;
else
    num_track = numel(dres_track.fr);
end

flag = zeros(num_track, 1);
for i = 1:num_track
    [~, o] = calc_overlap(dres_track, i, dres_track, 1:num_track);
    o(i) = 0;
    o(flag == 1) = 0;
    [mo, ind] = max(o);
    if mo > opt.overlap_sup
        num1 = trackers{dres_track.id(i)}.streak_tracked;
        num2 = trackers{dres_track.id(ind)}.streak_tracked;
        o1 = max(calc_overlap(dres_track, i, dres_det, 1:num_det));
        o2 = max(calc_overlap(dres_track, ind, dres_det, 1:num_det));

        if num1 > num2
            sup = ind;
        elseif num1 < num2
            sup = i;
        else
            if o1 > o2
                sup = ind;
            else
                sup = i;
            end
        end

        trackers{dres_track.id(sup)}.state = 3;
        trackers{dres_track.id(sup)}.dres.state(end) = 3;
        if opt.is_text
            fprintf('target %d suppressed\n', dres_track.id(sup));
        end
        flag(sup) = 1;
    end
end

author: yuxng


text:
function pdf = beta_pdf(x, a, b)
% PURPOSE: pdf of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: pdf = beta_pdf(x,a,b)
% where:   x = vector of components
%          a = beta distribution parameter, a = scalar
%          b = beta distribution parameter  b = scalar
% NOTE: mean[(beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: pdf at each element of x of the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox

if (nargin ~=3)
    error('Wrong # of arguments to beta_pdf');
end

if any(any((a<=0)|(b<=0)))
   error('Parameter a or b is nonpositive');
end

I = find((x<0)|(x>1));

pdf = x.^(a-1) .* (1-x).^(b-1) ./ beta(a,b);
pdf(I) = 0*I;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function sphere_llt_grid_display_test ( )

%*****************************************************************************80
%
%% SPHERE_LLT_GRID_DISPLAY_TEST tests SPHERE_LLT_GRID_DISPLAY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2015
%
%  Author:
%
%    John Burkardt
%
  lat_num = 10;
  long_num = 12;
  pc(1,1:3) = [ 0.0, 0.0, 0.0 ];
  r = 10.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_LLT_GRID_DISPLAY_TEST\n' );
  fprintf ( 1, '  SPHERE_LLT_GRID_DISPLAY displays an LLT grid on a sphere.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of latitudes is  %d\n', lat_num );
  fprintf ( 1, '  Number of longitudes is %d\n', long_num );
%
%  Get points.
%
  node_num = sphere_llt_grid_point_count ( lat_num, long_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of grid points is %d\n', node_num );

  node_xyz = sphere_llt_grid_points ( r, pc, lat_num, long_num, node_num );
%
%  Get lines.
%
  line_num = sphere_llt_grid_line_count ( lat_num, long_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of line segments is %d\n', line_num );

  line_data = sphere_llt_grid_lines ( lat_num, long_num, line_num );

  filename = 'sphere_llt_grid.png';

  sphere_llt_grid_display ( r, pc, lat_num, long_num, node_num, node_xyz, line_num, ...
    line_data, filename );

  return
end


text:
function plot_DR30_fa(plot_obj,plot_args,legend_string)
% Plots the point on the minDCF curve at which the number of false
% alarms falls below 30 (see Doddington's Rule of 30).
% Inputs:
%   plot_args: A cell array of arguments to be passed to 'plot' that control
%     the appearance of the curve. See Matlab's help on 'plot' for information.
%   legend_string: Optional. A string to describe this curve in the legend.
%     This string will be prepended with the system name.

if ischar(plot_args)
    plot_args = {plot_args};
end

figure(plot_obj.fh);
assert(iscell(plot_args))
lh = plot(plot_obj.plo(plot_obj.dr30FA),plot_obj.minDCF(plot_obj.dr30FA),plot_args{:});

if exist('legend_string','var') && ~isempty(legend_string)
    assert(ischar(legend_string))
    plot_obj.add_legend_entry(lh,legend_string,true);
end

end

author: nesl


text:
function xs = searchSeries(x,gpstruct,LB,UB,optimState,options)
%SEARCHSERIES Search step by time series prediction (UNSUPPORTED).

if nargin < 1
    xs = 'srs';
    return;
end

MeshSize = optimState.meshsize;
SearchFactor = optimState.searchfactor;

xsuccess = bsxfun(@rdivide,optimState.xsuccess,optimState.scale);
fsuccess = optimState.fsuccess;

D = size(x,2);
t = size(xsuccess,1);

% Random rotation matrix
[R,~] = qr(randn(D));
R = eye(D);

fmu = zeros(1,D);
fs = zeros(1,D);

xsuccess = xsuccess*R';

for i = 1:D
    index = max(1,t-10):t;

    stdy = sqrt(std(xsuccess(index,i)).^2 + optimState.TolMesh.^2);
    gptemp = gpinit(stdy,options,optimState);
    gptemp.x = (index)';
    % gptemp.x = -fsuccess(index);
    gptemp.y = xsuccess(index,i);
    gptemp.hyp(1).lik = log(MeshSize);
    gptemp.prior.lik{1}{2} = gptemp.hyp(1).lik;

    if numel(index) < 3 || 1
        lastindex = 1:numel(index);
    else
        lastindex = ceil(numel(index)/2):numel(index);
    end
    if numel(gptemp.hyp(1).mean) == 1
        gptemp.hyp(1).mean = mean(gptemp.y(lastindex));
        gptemp.prior.mean{1}{2} = gptemp.hyp(1).mean;
    else
        gptemp.hyp(1).mean = polyfit(gptemp.x(lastindex),gptemp.y(lastindex),1)';
        gptemp.prior.mean{1}{2} = gptemp.hyp(1).mean(1);
        gptemp.prior.mean{2}{2} = gptemp.hyp(1).mean(2);
    end

    bounds = unwrap2vec(gptemp.bounds);
    lb = bounds(1:2:end-1);
    ub = bounds(2:2:end);

    prior = @(theta_) independent_prior(gptemp.prior,theta_);
    gptemp.inf = {@inference_with_prior, @exact_inference, prior};

    % hyp0(1) = gptemp.hyp;
    %s = min(max(unwrap2vec(independent_prior(gptemp.prior)),lb),ub);
    %hyp0(2) = rewrap(gptemp.hyp,s);

    % Quick-and-dirty grid optimization
    % ell = 0:0.1:2;
    ell = 1;
    sf = log(stdy) + (-2:0.5:2);
    theta = combvec(ell,sf)';
    nlZ = zeros(1,size(theta,1));
    for j = 1:size(theta,1)
        nlZ(j) = gp_optimizer(theta(j,:), gptemp.hyp, gptemp.inf, gptemp.mean, gptemp.cov, gptemp.lik, gptemp.x, gptemp.y);
    end
    [~,index] = min(nlZ);
    gptemp.hyp.cov = theta(index,:)';

    cov(i,:) = gptemp.hyp.cov(:)';

    % optoptions = optimset('TolFun',0.01,'TolX',1e-4,'MaxFunEval',100);
    % optoptions.Hessian = 'user-supplied';
    % gptemp.hyp = gpHyperOptimize(hyp0,gptemp,optoptions);

%   [~,~,fmu(i),fs2] = gpgrad(0,gptemp);
    dt = exp(SearchFactor-1);
    [~,~,fmu(i),fs2] = gpgrad(t+dt,gptemp);
    fs(i) = sqrt(fs2);
end

% cov

fs = sqrt(fs.^2 + (MeshSize*options.PollMeshMultiplier^(-options.SearchGridNumber)).^2);

xs = bsxfun(@plus,bsxfun(@times,randn(options.Nsearch,D),fs),fmu);

xs = xs*R;

if 1
    figure(5);
    hold off;
    %plot(-fsuccess,xsuccess(:,1:end-1));
    plot(1:t,xsuccess*R);
    hold on;
    scatter((t+dt)*ones(size(fmu)),fmu);
    %scatter(0*ones(size(fmu)),fmu);
    drawnow;
    set(gcf,'Color','w');
    set(gca,'TickDir','out');
end

end

%--------------------------------------------------------------------------
function gptemp = gpinit(stdy,options,optimState)
%BGA_GPINIT

TolFun = options.TolFun;
TolMesh = optimState.TolMesh;

%% gp covariance function

% Base covariance function properties: smoothness and length scales
% gptemp.cov = {@isotropic_matern_covariance, 5};
gptemp.cov = {@isotropic_sqdexp_covariance};
gptemp.covbase = gptemp.cov;

ncov = str2double(feval(gptemp.cov{1}));
gptemp.hyp.cov = zeros(ncov,1);

% Covariance priors
for i = 1:ncov-1
     gptemp.prior.cov{i} = {@gaussian_prior, 1, 2^2};
     gptemp.bounds.cov{i} = [-5;Inf];
end
gptemp.hyp.cov = [log(3);log(stdy)];
gptemp.prior.cov{ncov} = {@gaussian_prior, log(stdy), 2^2};
gptemp.bounds.cov{ncov} = [log(TolFun);Inf];

if 0
    gptemp.cov = {@covSum,{gptemp.covbase,{@covPERiso,{@covMaterniso,Smoothness*2+1}}}};
    gptemp.hyp.cov = [gptemp.hyp.cov;zeros(ncov+1,1)];
    for i = ncov+(1:ncov)
         gptemp.prior.cov{i} = {@gaussian_prior,-1,2^2};
         gptemp.bounds.cov{i} = [-5;Inf];
    end
    gptemp.prior.cov{ncov*2+1} = {@gaussian_prior,0,2^2};
    gptemp.bounds.cov{ncov*2+1} = [log(TolFun);Inf];
elseif 0
    gptemp.cov = {@covariance_product,{gptemp.covbase,@NNone_covariance}};
    gptemp.hyp.cov = [gptemp.hyp.cov;zeros(ncov,1)];
    for i = ncov+(1:ncov-1)
         gptemp.prior.cov{i} = {@gaussian_prior,-1,2^2};
         gptemp.bounds.cov{i} = [-5;Inf];
    end
    gptemp.prior.cov{ncov*2} = {@gaussian_prior,0,2^2};
    gptemp.bounds.cov{ncov*2} = [log(TolFun);Inf];
end

%% gp likelihood function

gptemp.lik = @likGauss;   % Gaussian likelihood
likmu = log(TolMesh); liks2 = 1^2;
gptemp.hyp.lik = likmu;
gptemp.prior.lik = {{@delta_prior, likmu}};
% gptemp.prior.lik = {{@gaussian_prior, likmu, liks2}};

gptemp.hyp.lik = gptemp.hyp.lik(:);
gptemp.bounds.lik = {[log(TolMesh); 10]};

%% gp mean function
if 1
    gptemp.mean = @constant_mean;             % Constant mean function
    gptemp.hyp.mean = NaN;
    gptemp.prior.mean = {{@delta_prior, 0}};      % Fixed mean
else
    gptemp.mean = {@meanSum,{@meanLinear,@meanConst}};
    gptemp.hyp.mean = NaN(2,1);
    gptemp.prior.mean{1} = {@delta_prior, 0};      % Fixed prior
    gptemp.prior.mean{2} = {@delta_prior, 0};
end

%% gp sampling weight
gptemp.hypweight = 1;
gptemp.hypmean = [];

%% gp inference method
gpinf = @inference_with_prior;              % Exact inference
gptemp.infMethod = 'exact';

prior = @(theta_) independent_prior(gptemp.prior,theta_);
gptemp.inf = {@inference_with_prior, @exact_inference, prior};

% gptemp.inf = {@infPrior,gpinf,gptemp.prior};

gptemp = gpstruct_check(gptemp, 1);

end

%--------------------------------------------------------------------------
function nlZ = gp_optimizer(theta, hyp, inf, mean, cov, lik, x, y)
%GP_OPTIMIZER Wrapper function for GP optimization

hyp.cov = theta(:);
f = @() (feval(inf{:}, hyp, mean, cov, lik, x, y));
[~, nlZ] = f();

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
% compare geodesic with euclidean distance inside a 2D shape

rep = 'images/';
name = 'bird';
name = 'giraffe';
name = 'camel';
name = 'chicken';

test = 'meangeodesic';
test = 'eccentricity';

% put s=-1 for Linf eccentricity
s = -1;

repsrc = 'data/';

rep = 'results/geodesic-vs-euclidean/';
if ~exist(rep)
    mkdir(rep);
end

n = 256;
M = rescale( load_image([repsrc name],n), 0,1 );
M = double(M>0.5);

% make sure pixels on the boundary are black
if M(1)==1
    M = 1-M;
end

% compare the geodesic ditance to the euclidean distance
clf;
imagesc(M); axis image; axis off;
[y,x] = ginput(1);
start_points = round([x y]');
W = ones(n);
L = zeros(n)-Inf; L(M==1) = +Inf;
options.constraint_map = [];
[D0,S,Q] = perform_fast_marching(W, start_points, options);
options.constraint_map = L;
[D1,S,Q] = perform_fast_marching(W, start_points, options);
D0(M==0) = Inf; D1(M==0) = Inf;

% compute a nice color image
D0a = convert_distance_color(D0);
D1a = convert_distance_color(D1);

% place a dot at the center
[Y,X] = meshgrid(1:n,1:n);
r = 0.02*n;
I = find( (X-x).^2 + (Y-y).^2 <= r^2 );
for i=1:3
    A = D0a(:,:,i); A(I) = i==1;
    D0a(:,:,i) = A;
    A = D1a(:,:,i); A(I) = i==1;
    D1a(:,:,i) = A;
end

% display
clf;
subplot(1,2,1)
imagesc(D0a); axis image; axis off; title('Euclidean');
subplot(1,2,2)
imagesc(D1a); axis image; axis off; title('Geodesic');

warning off;
imwrite(rescale(D0a), [rep name '-euclidean.png' ], 'png');
imwrite(rescale(D1a), [rep name '-geodesic.png' ], 'png');
warning on;

author: alecjacobson

repo: gptoolbox


text:
classdef BackgroundSubtractorMOG < handle
    %BACKGROUNDSUBTRACTORMOG  Gaussian Mixture-based Background/Foreground Segmentation Algorithm
    %
    % The class implements the algorithm described in [KB2001].
    %
    % ## References
    % [KB2001]:
    % > Pakorn KaewTraKulPong and Richard Bowden. "An improved adaptive
    % > background mixture model for real-time tracking with shadow detection".
    % > In Video-Based Surveillance Systems, pages 135-144. Springer, 2002.
    % > [PDF](http://personal.ee.surrey.ac.uk/Personal/R.Bowden/publications/avbs01/avbs01.pdf)
    %
    % See also: cv.BackgroundSubtractorMOG.BackgroundSubtractorMOG,
    %  cv.BackgroundSubtractorMOG.apply,
    %  cv.BackgroundSubtractorMOG.getBackgroundImage
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % The number of last frames that affect the background model.
        History
        % The number of gaussian components in the background model.
        % The model needs to be reinitalized to reserve memory.
        NMixtures
        % The "background ratio" parameter of the algorithm.
        % If a foreground pixel keeps semi-constant value for about
        % `BackgroundRatio * History` frames, it's considered background and
        % added to the model as a center of a new component. It corresponds to
        % `TB` parameter in the paper.
        BackgroundRatio
        % Sigma of noise
        NoiseSigma
    end

    %% BackgroundSubtractor
    methods
        function this = BackgroundSubtractorMOG(varargin)
            %BACKGROUNDSUBTRACTORMOG  Creates mixture-of-gaussian background subtractor
            %
            %     bs = cv.BackgroundSubtractorMOG()
            %     bs = cv.BackgroundSubtractorMOG('OptionName', optionValue, ...)
            %
            % ## Options
            % * __History__ Length of the history. default 200
            % * __NMixtures__ Number of Gaussian mixtures. default 5
            % * __BackgroundRatio__ Background ratio. default 0.7
            % * __NoiseSigma__ Noise strength (standard deviation of the
            %   brightness or each color channel). 0 means some automatic
            %   value. default 0
            %
            % Default constructor sets all parameters to default values.
            %
            % See also: cv.BackgroundSubtractorMOG
            %
            this.id = BackgroundSubtractorMOG_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     bs.delete()
            %
            % See also: cv.BackgroundSubtractorMOG
            %
            if isempty(this.id), return; end
            BackgroundSubtractorMOG_(this.id, 'delete');
        end

        function fgmask = apply(this, im, varargin)
            %APPLY  Updates the background model and computes the foreground mask
            %
            %     fgmask = bs.apply(im)
            %     fgmask = bs.apply(im, 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __im__ Next video frame.
            %
            % ## Output
            % * __fgmask__ The output foreground mask as an 8-bit binary image
            %   (0 for background, 255 for foregound).
            %
            % ## Options
            % * __LearningRate__ The value between 0 and 1 that indicates how
            %   fast the background model is learnt. Negative parameter value
            %   makes the algorithm to use some automatically chosen learning
            %   rate. 0 means that the background model is not updated at all,
            %   1 means that the background model is completely reinitialized
            %   from the last frame. default -1
            %
            % See also: cv.BackgroundSubtractorMOG.getBackgroundImage
            %
            fgmask = BackgroundSubtractorMOG_(this.id, 'apply', im, varargin{:});
        end

        function bgImg = getBackgroundImage(this)
            %GETBACKGROUNDIMAGE  Computes a background image
            %
            %     bgImg = bs.getBackgroundImage()
            %
            % ## Output
            % * __bgImg__ The output background image.
            %
            % ### Note
            % Method not implemented for this class, throws exception.
            %
            % See also: cv.BackgroundSubtractorMOG.apply
            %
            bgImg = BackgroundSubtractorMOG_(this.id, 'getBackgroundImage');
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.BackgroundSubtractorMOG.empty
            %
            BackgroundSubtractorMOG_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.BackgroundSubtractorMOG.clear
            %
            b = BackgroundSubtractorMOG_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BackgroundSubtractorMOG.save, cv.BackgroundSubtractorMOG.load
            %
            name = BackgroundSubtractorMOG_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.BackgroundSubtractorMOG.load
            %
            BackgroundSubtractorMOG_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.BackgroundSubtractorMOG.save
            %
            BackgroundSubtractorMOG_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.History(this)
            value = BackgroundSubtractorMOG_(this.id, 'get', 'History');
        end
        function set.History(this, value)
            BackgroundSubtractorMOG_(this.id, 'set', 'History', value);
        end

        function value = get.NMixtures(this)
            value = BackgroundSubtractorMOG_(this.id, 'get', 'NMixtures');
        end
        function set.NMixtures(this, value)
            BackgroundSubtractorMOG_(this.id, 'set', 'NMixtures', value);
        end

        function value = get.BackgroundRatio(this)
            value = BackgroundSubtractorMOG_(this.id, 'get', 'BackgroundRatio');
        end
        function set.BackgroundRatio(this, value)
            BackgroundSubtractorMOG_(this.id, 'set', 'BackgroundRatio', value);
        end

        function value = get.NoiseSigma(this)
            value = BackgroundSubtractorMOG_(this.id, 'get', 'NoiseSigma');
        end
        function set.NoiseSigma(this, value)
            BackgroundSubtractorMOG_(this.id, 'set', 'NoiseSigma', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Computes the probability density function of the multivariate gaussian distribution.
function probabilities = multivariate_gaussian(X, mu, sigma2)
    % Get number of training sets and features.
    [m n] = size(X);

    % Init probabilities matrix.
    probabilities = ones(m, 1);

    % Go through all training examples and through all features.
    for i=1:m
        for j=1:n
            p = (1 / sqrt(2 * pi * sigma2(j))) * exp(-(X(i, j) - mu(j)) .^ 2 / (2 * sigma2(j)));
            probabilities(i) = probabilities(i) * p;
        end
    end
end

author: trekhleb

repo: machine-learning-octave


text:
% revised on September 21, 2010 to make parameters consistent

% September 25, 2009
% written by Shuiwang Ji and Jieping Ye

% This function implements the accelerated gradient algorithm for
% multi-task learning regularized by trace norm as
% described in Ji and Ye (ICML 2009).

% References:
%Ji, S. and Ye, J. 2009. An accelerated gradient method for trace norm minimization.
%In Proceedings of the 26th Annual international Conference on Machine Learning
%(Montreal, Quebec, Canada, June 14 - 18, 2009). ICML '09, vol. 382. ACM, New York,
%NY, 457-464.

%[W, fval_vec,itr_counter] =
%accel_grad_mtl(A,Y,lambda,opt)

% required inputs:
% A: K x 1 cell in which each cell is N x D where N is the sample size
% and D is data dimensionality and K is the number of tasks
% Y: K x 1 cell in which each cell contains the output of the
% corresponding task
% lambda: regularization parameter

% optional inputs:
% opt.L0: Initial guess for the Lipschitz constant
% opt.gamma: the multiplicative factor for Lipschitz constant
% opt.W_init: initial weight matrix
% opt.epsilon: precision for termination
% opt.max_itr: maximum number of iterations

% outputs:
% W: the computed weight matrix
% fval_vec: a vector for the sequence of function values
% itr_counter: number of iterations executed

function [W,fval_vec,itr_counter] = accel_grad_mtl(A,Y,lambda,opt)

Xtrain = A;
Ytrain = Y;
clear A;
clear Y;

if nargin<4
    opt = [];
end

if isfield(opt, 'L0')
    L0 = opt.L0;
else
    L0 = 100;
end

if isfield(opt, 'gamma')
    gamma = opt.gamma;
else
    gamma = 1.1;
end

if isfield(opt, 'W_init')
    W_init = opt.W_init;
else
    W_init = zeros(size(Xtrain{1},2),length(Ytrain));
end

if isfield(opt, 'epsilon')
    epsilon = opt.epsilon;
else
    epsilon = 10^-5;
end

if isfield(opt, 'max_itr')
    max_itr = opt.max_itr;
else
    max_itr = 100;
end

W_old = W_init;
alpha = 1;
fval_vec = [];
L = L0;
fval_old = rand(1,1);
fval = rand(1,1);
itr_counter = 0;
Z_old = W_old;
while abs((fval_old-fval)/fval_old)>epsilon
    fval_old = fval;
    [Wp,P,sval] = ComputeQP(Xtrain,Ytrain,Z_old,L,lambda);
    f = 0;
    for i = 1:length(Xtrain)
        f = f+0.5*norm(Ytrain{i}-Xtrain{i}*Wp(:,i))^2;
    end
    fval = f+lambda*sval;
    Q = P+lambda*sval;
    while fval>Q
        fprintf('Searching step size (fval = %f, Q = %f)...\n',fval,Q);
        L = L*gamma;
        [Wp,P,sval] = ComputeQP(Xtrain,Ytrain,Z_old,L,lambda);
        f = 0;
        for i = 1:length(Xtrain)
            f = f+0.5*norm(Ytrain{i}-Xtrain{i}*Wp(:,i))^2;
        end
        fval = f+lambda*sval;
        Q = P+lambda*sval;
    end

    alpha_old = alpha;
    alpha = (1+sqrt(1+4*alpha_old^2))/2;
    Z_old = Wp+((alpha_old-1)/alpha)*(Wp-W_old);

    W_old = Wp;
    fval_vec = [fval_vec,fval];
    itr_counter = itr_counter+1;
    if itr_counter>max_itr
        break;
    end
    if mod(itr_counter,100)==0
        fprintf('Iteration = %8d,  objective = %f\n',itr_counter, fval);
    end

end
W = Wp;
return;

function [Wp,P,sval] = ComputeQP(X,Y,W,L,lambda)

[W1,delta_W] = ComputeGradStep(X,Y,W,L);
[U,D,V] = svd(W1,0);
D = diag(D);
D = D-(lambda/L);
idx = find(D>0);
sval = sum(D(idx));
Wp = U(:,idx)*diag(D(idx))*V(:,idx)';
P = 0;
for i = 1:length(X)
    P = P+0.5*norm(Y{i}-X{i}*W(:,i))^2;
end
P = P+trace(delta_W'*(Wp-W))+0.5*L*norm(Wp-W,'fro')^2;
return;

function [W1,delta_W] = ComputeGradStep(X,Y,W,L)
delta_W = ComputeDerivative(X,Y,W);
W1 = W-(1/L)*delta_W;
return;

function dev = ComputeDerivative(X,Y,W)

dev = [];
for i = 1:length(X)
    dev = [dev, X{i}'*(X{i}*W(:,i)-Y{i})];
end
return;

author: thomaskuestner


text:
function [vertices2, faces2] = subdivideMesh(vertices, faces, n)
%SUBDIVIDEMESH Subdivides each face of the mesh
%
%   [V2 F2] = subdivideMesh(V, F, N)
%   Subdivides the mesh specified by (V,F) such that each face F is divided
%   into N^2 smaller faces.
%
%   Example
%     [v f] = createOctahedron;
%     figure; drawMesh(v, f); view(3);
%     [v2 f2] = subdivideMesh(v, f, 4);
%     figure; drawMesh(v2, f2); view(3)
%
%   See also
%     meshes3d, drawMesh
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2013-08-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

%% Initialisations

if ~isnumeric(faces) || size(faces, 2) ~= 3
    error('Requires a triangular mesh');
end

% compute the edge array
% edgeVertexIndices = computeMeshEdges(faces);
edges = meshEdges(faces);
nEdges = size(edges, 1);

% index of faces around each edge
% edgeFaceIndices = meshEdgeFaces(vertices, edges, faces);

% index of edges around each face
faceEdgeIndices = meshFaceEdges(vertices, edges, faces);

%% Create new vertices on edges

% several interpolated positions
t = linspace(0, 1, n + 1)';
coef2 = t(2:end-1);
coef1 = 1 - t(2:end-1);

% initialise the array of new vertices
vertices2 = vertices;

% keep an array containing index of new vertices for each original edge
edgeNewVertexIndices = zeros(nEdges, n-1);

% create new vertices on each edge
for iEdge = 1:nEdges
    % extract each extremity as a point
    v1 = vertices(edges(iEdge, 1), :);
    v2 = vertices(edges(iEdge, 2), :);

    % compute new points
    newPoints = coef1 * v1 + coef2 * v2;

    % add new vertices, and keep their indices
    edgeNewVertexIndices(iEdge,:) = size(vertices2, 1) + (1:n-1);
    vertices2 = [vertices2 ; newPoints]; %#ok<AGROW>
end

%% Process each face

faces2 = zeros(0, 3);

nFaces = size(faces, 1);
for iFace = 1:nFaces
    % compute index of each corner vertex
    face = faces(iFace, :);
    iv1 = face(1);
    iv2 = face(2);
    iv3 = face(3);

    % compute index of each edge
    faceEdges = faceEdgeIndices{iFace};
    ie1 = faceEdges(1);
    ie2 = faceEdges(2);
    ie3 = faceEdges(3);

    % indices of new vertices on edges
    edge1NewVertexIndices = edgeNewVertexIndices(ie1, :);
    edge2NewVertexIndices = edgeNewVertexIndices(ie2, :);
    edge3NewVertexIndices = edgeNewVertexIndices(ie3, :);

    % keep vertex 1 as reference for edges 1 and 3
    if edges(ie1, 1) ~= iv1
        edge1NewVertexIndices = edge1NewVertexIndices(end:-1:1);
    end
    if edges(ie3, 1) ~= iv1
        edge3NewVertexIndices = edge3NewVertexIndices(end:-1:1);
    end

    % create the first new face, on 'top' of the original face
    topVertexInds = [edge1NewVertexIndices(1) edge3NewVertexIndices(1)];
    newFace = [iv1 topVertexInds];
    faces2 = [faces2; newFace]; %#ok<AGROW>

    % iterate over middle strips
    for iStrip = 2:n-1
        % index of extreme vertices of current row
        ivr1 = edge1NewVertexIndices(iStrip);
        ivr2 = edge3NewVertexIndices(iStrip);

        % extreme vertices as points
        v1 = vertices2(ivr1, :);
        v2 = vertices2(ivr2, :);

        % create additional vertices within the bottom row of the strip
        t = linspace(0, 1, iStrip+1)';
        coef2 = t(2:end-1);
        coef1 = 1 - t(2:end-1);
        newPoints = coef1 * v1 + coef2 * v2;

        % compute indices of new vertices in result array
        newInds = size(vertices2, 1) + (1:iStrip-1);
        botVertexInds = [ivr1 newInds ivr2];

        % add new vertices
        vertices2 = [vertices2 ; newPoints]; %#ok<AGROW>

        % create top faces of current strip
        for k = 1:iStrip-1
            newFace = [topVertexInds(k) botVertexInds(k+1) topVertexInds(k+1)];
            faces2 = [faces2; newFace]; %#ok<AGROW>
        end

        % create bottom faces of current strip
        for k = 1:iStrip
            newFace = [topVertexInds(k) botVertexInds(k) botVertexInds(k+1)];
            faces2 = [faces2; newFace]; %#ok<AGROW>
        end

        % bottom vertices of current strip are top vertices of next strip
        topVertexInds = botVertexInds;
    end

    % for edge 2, keep vertex 2 of the current face as reference
    if edges(ie2, 1) ~= iv2
        edge2NewVertexIndices = edge2NewVertexIndices(end:-1:1);
    end

    % consider new vertices together with extremities
    botVertexInds = [iv2 edge2NewVertexIndices iv3];

    % create top faces for last strip
    for k = 1:n-1
        newFace = [topVertexInds(k) botVertexInds(k+1) topVertexInds(k+1)];
        faces2 = [faces2; newFace]; %#ok<AGROW>
    end

    % create bottom faces for last strip
    for k = 1:n
        newFace = [topVertexInds(k) botVertexInds(k) botVertexInds(k+1)];
        faces2 = [faces2; newFace]; %#ok<AGROW>
    end

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips