text:
%--- help for read_identification_restrictions ---
%
%  INTERNAL FUNCTION: Formats identification restrictions
%
%  See also:
%     :func:`rfvar.identification`
%
%

author: jmaih


text:
% Question No: 6

% Consider an image composed of small, no overlapping blobs. Segmenting the
% blobs based on region growing.

function regrow(x)
f1=imread(x);
f=double(f1);
s=255;
t=65;
if numel(s)==1
    si=f==s;
    s1=s;
else
    si=bwmorph(s,'shrink',Inf);
    j=find(si);
    s1=f(j);
end
ti=false(size(f));
for k=1:length(s1)
    sv=s1(k);
    s=abs(f-sv)<=t;
    ti=ti|s;
end
[g nr]=bwlabel(imreconstruct(si,ti));
figure,imshow(f1),title('Original Image');
figure,imshow(g),title('Segmented Image - Region Growing');
display('No: of regions');
nr
end

author: Sable

repo: mcbench-benchmarks


text:
classdef TestBriefDescriptorExtractor
    %TestBriefDescriptorExtractor

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
    end

    methods (Static)
        function test_compute_img
            obj = cv.BriefDescriptorExtractor('Bytes',32);
            typename = obj.typeid();
            ntype = obj.defaultNorm();

            img = imread(TestBriefDescriptorExtractor.im);
            kpts = cv.FAST(img, 'Threshold',20);
            [desc, kpts2] = obj.compute(img, kpts);
            validateattributes(kpts2, {'struct'}, {'vector'});
            assert(all(ismember(fieldnames(kpts), fieldnames(kpts2))));
            validateattributes(desc, {obj.descriptorType()}, ...
                {'size',[numel(kpts2) obj.descriptorSize()]});
        end

        function test_compute_imgset
            img = imread(TestBriefDescriptorExtractor.im);
            imgs = {img, img};
            kpts = cv.FAST(img, 'Threshold',20);
            kpts = {kpts, kpts};

            obj = cv.BriefDescriptorExtractor();
            [descs, kpts] = obj.compute(imgs, kpts);
            validateattributes(kpts, {'cell'}, {'vector'});
            validateattributes(descs, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(d,k) validateattributes(d, {obj.descriptorType()}, ...
                {'size',[numel(k) obj.descriptorSize()]}), descs, kpts);
        end

        function test_error_1
            try
                cv.BriefDescriptorExtractor('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function varargout = plotss(varargin)
% VL_PLOTSS Plot scale space
%   VL_PLOTSS(SS) plots the scale space SS. SS is a structure
%   with the following members:
%
%   ss.firstOctave::
%     The index of the first octave in the scale space.
%
%   ss.lastOctave::
%     The index of the last octave in the scale space.
%
%   ss.octaveResolution::
%     The octave resolution, i.e. the nubmer of subdivisions
%     per octave.
%
%   ss.octaveFirstLevel::
%     The index of the first level of subdivisions for each octave.
%
%   ss.octaveLastLevel::
%     The iundex of last leve of subdivisions for each cotave.
%
%   ss.data::
%     A cell array of 3D arrays representing the scale space data.
%     The cell array has a length equal to the nubmer of octaves
%     contained in the scale space. Each entry is a 3D array, the
%     first two dimensions of which correspond to image rows and
%     columns respectively, and the third to scale levels.
%
%   ss.sigma0::
%     Base smoothing.
%
%   A scale space is a representation of a 2D signal (image) at
%   multiple scales. In the simplest case, a scale SIGMA is defined as
%   the input image I(x,y) convolved by a Gaussian kernel of isotropic
%   standard deviation SIGMA:
%
%     I(x,y;sigma) = (g_sigma * I)(x,y)
%
%   where scales are sampled as follows:
%
%     sigma(o,s) = sigma0 2^{o + s / ss.octaveResolution),
%     ss.firstOctave <= o <= ss.lastOctave,
%     ss.octaveFirstLeve <= s <= ss.octaveLastLevel.
%
%   Moving from one octave to the next, the size of the kernel
%   doubles. Hence the effective bandwith of the signal halves, and
%   resolution can be reduced by half. Typically, for example, sigma0
%   = 1.6, so at octave 0 the image can be effectively sampled with a
%   step of 1, and the resolution of octave 0 is the same as the one
%   at which the input image is presented. Then at octave o the
%   sampling step is 2^o.
%
%   ss.octaveResolution is the number of scale subdivisions per
%   octave. ss.firstOctave and ss.lastOctave give the additional
%   flexibility of specifying a range for the level index s to exceed
%   the standard setting [0, ss.octaveResolution-1]. In this manner
%   the same scales can be represented twice, at two sampling
%   rates. This is often convenient in feature computation (e.g. to
%   find local maxima in scale of a function).
%
%   VL_PLOTSS(SS, 'Option', value) supports the following options:
%
%   Uniform:: false
%     If TRUE then use a fixed gray scale for all the levels.
[varargout{1:nargout}] = vl_plotss(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function outputImage = rms(this, varargin)
% Computes root mean square along specified dimension, i.e. sqrt(mean(Y.^2))
%
%
%   Y = MrImage()
%   Y.rms(applicationDimension)
%
% This is a method of class MrImage.
%
% IN
%   applicationDimension    image dimension along which operation is
%                           performed (e.g. 4 = time, 3 = slices)
%                           default: The last dimension with more than one
%                           value is chosen
%                           (i.e. 3 for 3D image, 4 for 4D image)
%
% OUT
%   outputImage             rms of all images along application dimension
%
% EXAMPLE
%   rms
%
%   See also MrImage MrImage.perform_unary_operation

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-12-23
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin > 1
    applicationDimension = varargin{1};
    outputImage = mean(this.^2, applicationDimension).^(1/2);
else
    outputImage = mean(this.^2).^(1/2);
end

author: translationalneuromodeling

repo: tapas


text:
function x = is_distributed(x)
x = is_1d_distributed(x) || is_Nd_distributed(x) || is_discrete(x) || is_custom_distributed(x);

repo: REST


text:
function [suspicious_index suspicious_score u] = OD_onlinePCA(A, beta)
%
% Outlier detection via over-sampling online PCA
%
% This function is used for outlier detection. The main idea is using
% the variation of the first principal direction detect the outlierness
% of each instance(event) in the leave one out procedure. Here the
% over-sampling on target instance is also used for enlarge the
% outlierness
%
%
% input
%   A: the data matrix, each row represents an instance
%   beta: forgetting factor
%         For example, beta=0.9 means we decrease the influence of previous
%         data by a factor of 0.9
% output
%   suspicious_score: the suspicious score for each instance
%   suspicious_index: the ranking of instances according to their
%                     suspicious score
%                     For example, suspicious_index(i)=j means the ith
%                     instance is in jth position in the ranking.
%
% References
% Anomaly Detection via Over-Sampling Principal Component Analysis
% Send your comment and inquiry to Dr. Yi-Ren Yeh (yirenyeh@gmail.com)
%

[n p]= size(A);

A_m = mean(A);
d = 0.0001;
u = ones(p,1);
for i = 1:n
    [u d] = Track_w(A(i,:)-A_m, u, d, 1);
end
u = u/norm(u);

sim_pool = zeros(n,1);
ratio = 1/(n*beta);
for i = 1:n
    temp_mu = (A_m+ratio*A(i,:))/(1+ratio);
    x =A(i,:)-temp_mu;
    [w] = Track_w(x, u, d, beta);
    w = w/norm(w);
    sim_pool(i,:) = abs(diag(u'*w));
    if (~mod(i,10000))
        display(['Iteration ',num2str(i)])
    end
end

[non,suspicious_index]=sort(sim_pool);
suspicious_score = 1-sim_pool;

%==========================================================================
function [w d] = Track_w(x, w, d, beta)
y = x*w;
d = beta*d+y^2;
e = x'-w*y;
w = w + (y*e)/d;
% function [w d] = Track_w(x, w, d, n, ratio)
% y = x*w;
% d = (1/n)*d+ratio*y^2;
% e = x'-w*y;
% w = w + (y*e)/d;

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function [ x, seed ] = pearson_05_sample ( a, b, c, seed )

%*****************************************************************************80
%
%% PEARSON_05_SAMPLE samples the Pearson 5 PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < A, 0.0 < B.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a2 = 0.0;
  b2 = b;
  c2 = 1.0 / a;

  [ x2, seed ] = gamma_sample ( a2, b2, c2, seed );

  x = c + 1.0 / x2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Bubble sort algorithm:
function list = bubble_sort(list)
% function to sort vector 'list' with using the 'Bubble sort' algorithm
% INPUT: 'list' array
% OUTPUT: sorted array
changed = true;
count = numel(list);
while(changed)
    changed = false;
    count = count - 1;
    for index = (1:count)
        if(list(index) > list(index+1))
            list([index index+1]) = list([index+1 index]); %swap
            changed = true;
        end

    end
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function Q=getbasematrix(X,ind)
%GETBASEMATRIX Internal function to extract basematrix for variable IND

if ind==0
    base = X.basis(:,1);
    Q = reshape(base,X.dim);
    return;
end

here = find(X.lmi_variables==ind);
if isempty(here)
    error
else
    base = X.basis(:,here+1);
    Q = reshape(base,X.dim);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef BoundaryConditions < handle

    properties (GetAccess = public)
        dirichlet
        dirichlet_values
        free
        freeFields
        masterSlave
        periodic_free
        periodic_constrained
        neumann
        neumann_values
    end

    properties (Access = private)
        ndimf
        ndofs
        scale
        dirichletInput
        pointloadInput
    end

    methods (Access = public)

        function obj = BoundaryConditions(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            [dirID, dirVals]     = obj.formatInputData(obj.ndimf, obj.dirichletInput);
            [neuID, neuVals]     = obj.formatInputData(obj.ndimf, obj.pointloadInput);
            obj.dirichlet        = dirID;
            obj.dirichlet_values = dirVals;
            obj.neumann          = neuID;
            obj.neumann_values   = neuVals;
            obj.free             = obj.computeFreeDOF();
        end

        function red = fullToReducedMatrix(obj, mat)
            switch obj.scale
                case 'MACRO'
                    red = obj.reduceMatrixDirichlet(mat);
                case 'MICRO'
                    red = obj.reduceMatrixPeriodic(mat);
            end
        end

        function red = fullToReducedVector(obj, vec)
            switch obj.scale
                case 'MACRO'
                    red = obj.reduceVectorDirichlet(vec);
                case 'MICRO'
                    red = obj.reduceVectorPeriodic(vec);
            end
        end

        function full = reducedToFullVector(obj, vec)
            switch obj.scale
                case 'MACRO'
                    full = obj.expandVectorDirichlet(vec);
                case 'MICRO'
                    full = obj.expandVectorPeriodic(vec);
            end
        end

        function changeBoundaryConditions(obj,neumannDOFs,neumannValues)
            obj.neumann        = neumannDOFs;
            obj.neumann_values = neumannValues;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.scale          = cParams.scale;
            obj.ndofs          = cParams.ndofs; % Stokes
            obj.ndimf          = cParams.bc{1}.ndimf; % Elastic
            obj.initPeriodicMasterSlave(cParams);
            obj.initDirichletInput(cParams);
        end

        function initDirichletInput(obj, s)
            nfields = numel(s.bc);
            dirich  = [];
            neumann    = [];
            dirichVals = [];
            neumnnVals = [];
            free       = [];
            globalNdof = 0;
            for i = 1:nfields
                bc = s.bc{i};
                obj.dirichletInput = bc.dirichlet;
                obj.pointloadInput = bc.pointload;

                inD = bc.dirichlet;
                inN = bc.pointload;
                [idxD, valD] = obj.formatInputData(bc.ndimf,inD);
                [idxN, valN] = obj.formatInputData(bc.ndimf,inN);
                idxD = idxD + globalNdof;
                idxN = idxN + globalNdof;

                dirich     = [dirich; idxD];
                dirichVals = [dirichVals; valD];
                neumann    = [neumann; idxN];
                neumnnVals = [neumnnVals; valN];

                firstDof = globalNdof + 1;
                lastDof  = firstDof + bc.ndofs - 1;
                obj.freeFields{i} = setdiff(firstDof:lastDof,idxD);
                globalNdof = globalNdof+ bc.ndofs;
            end
            obj.dirichlet        = dirich;
            obj.dirichlet_values = dirichVals;
            obj.neumann          = neumann;
            obj.neumann_values   = neumnnVals;
            obj.free             = obj.computeFreeDOF();
        end

        function initPeriodicMasterSlave(obj, cParams)
            switch obj.scale
                case 'MICRO'
                    if isfield(cParams.bc, 'masterSlave')
                        obj.masterSlave = cParams.bc.masterSlave;
                    end
                    MS = obj.masterSlave;
                    if isempty(MS)
                        mesh = cParams.mesh;
                        MS = obj.computeMasterSlave(mesh.coord);
                        obj.masterSlave = MS;
                    end
                    obj.periodic_free = obj.computePeriodicNodes(MS(:,1));
                    obj.periodic_constrained = obj.computePeriodicNodes(MS(:,2));
            end
        end

        function perDof = computePeriodicNodes(obj,perNodes)
            nunkn = obj.ndimf;
            nlib = size(perNodes,1);
            perDof = zeros(nlib*nunkn,1);
            for iunkn = 1:nunkn
                indDof = nlib*(iunkn - 1) + [1:nlib];
                perDof(indDof,1) = obj.nod2dof(obj.ndimf, perNodes,iunkn);
            end
        end

        function free = computeFieldFree(obj, ndofs, dirich)
            free  = setdiff(1:ndofs,dirich);
        end

        function free = computeFreeDOF(obj)
            ndof  = obj.ndofs;
            cnstr = [obj.periodic_constrained;obj.dirichlet];
            free  = setdiff(1:ndof,cnstr);
        end

        function [dofs, vals] = formatInputData(obj, ndimf, data)
            dofs = [];
            vals = [];
            if ~isempty(data)
                inod = data(:,1);
                iunk = data(:,2);
                vals = data(:,3);
                dofs = obj.nod2dof(ndimf, inod,iunk);
            end
        end

        function idof = nod2dof(obj, ndimf, inode, iunkn)
            idof(:,1)= ndimf*(inode - 1) + iunkn;
        end

        function Ared = reduceMatrixDirichlet(obj,A)
%             fr = obj.computeGlobalFree();
            fr = obj.free';
            Ared = A(fr,fr);
        end

        function b_red = reduceVectorDirichlet(obj,b)
            fr = obj.free';
            b_red = b(fr);
        end

        function Ared = reduceMatrixPeriodic(obj,A)
            MS = obj.masterSlave;
            vF = obj.free;
            vP = obj.computePeriodicNodes(MS(:,1));
            vQ = obj.computePeriodicNodes(MS(:,2));
            vI = setdiff(vF,vP);

            A_II = A(vI,vI);
            A_IP = A(vI,vP) + A(vI,vQ); %Grouping P and Q nodal values
            A_PI = A(vP,vI) + A(vQ,vI); % Adding P  and Q equation
            A_PP = A(vP,vP) + A(vP,vQ) + A(vQ,vP) + A(vQ,vQ); % Adding and grouping

            Ared = [A_II, A_IP; A_PI, A_PP];
        end

        function b_red = reduceVectorPeriodic(obj,b)
            vF = obj.free;
            vP = obj.periodic_free;
            vQ = obj.periodic_constrained;
            vI = setdiff(vF,vP);
            b_I = b(vI);
            b_P = b(vP) + b(vQ);
            b_red = [b_I; b_P];
        end

        function b = expandVectorDirichlet(obj,bfree)
            dir = obj.dirichlet;
            uD  = obj.dirichlet_values;
            fr  = obj.free;
            nsteps = length(bfree(1,:));
            ndof = sum(obj.ndofs);
            uD = repmat(uD,1,nsteps);

            b = zeros(ndof,nsteps);
            b(fr,:) = bfree;
            if ~isempty(dir)
                b(dir,:) = uD;
            end
        end

        function b = expandVectorPeriodic(obj,bfree)
            vF = obj.free;
            vP = obj.periodic_free;
            vC = obj.periodic_constrained;
            vI = setdiff(vF,vP);
            b = zeros(obj.ndofs,1);
            b(vI) = bfree(1:1:size(vI,2));
            b(vP) = bfree(size(vI,2)+1:1:size(bfree,1));
            b(vC) = b(vP);
        end

        function MS = computeMasterSlave(obj, coord)
           mR = MasterSlaveRelator(coord);
           MS = mR.getRelation();
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [params] = sv_calcMcNodes_2methods(params,nNodeStart, nNodeEnd)
% function [params] = sv_calcMcNodes_2methods(params,nNodeStart, nNodeEnd);
% ----------------------------------------------------------------
% Calculation of magnitude of completeness specifying the nodes to be calculated for distributing on different CPUs
% Same as sv_calcMcNodes but just using McEMR and MAXC
%
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bRandom            Perform random simulation (=1) or real calculation (=0)
%   params.nCalculation       Number of random simulations
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fMaxRadius         Maximum Radius using a constant number of events; works only with bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
%   params.fMinMag            Lower limit of magnitude range for testing
%   params.fMaxMag            Upper limit of magnitude range for testing
%   params.bTimePeriod        Calculate seismicity difference for 2 periods (0) until start and end of catalog or
%                             a specific time period before and after fSplitTime (1)
%   params.fTimePeriod        Length of time periods
%   params.bTstart            Check for starting time of temporal mapping
%   params.fTstart            Starting time for temporal mapping
%   params.bBstnum            Check for boostrap sampling
%   params.fBstnum            Number of bootstrap samples
%   params.fBinning           Bin size for magnitude binning
%   params.sComment           Comment on calculation

% Output parameters:
%   Same as input parameters including
%   params.mValueGrid         Matrix of calculated values
%   params.vcsGridNames       Names of parameters calculated
%   Check sv_NodeCalcMc.m for a list of variables!!
%
% J. Woessner; j.woessner@sed.ethz.ch
% last update: 27.07.04

global bDebug;
if bDebug
    report_this_filefun(mfilename('fullpath'));
end

% Check calculation for splitting of nodes
if nargin < 2
    nNodeStart = 1;
    nNodeEnd = length(params.mPolygon(:,1));
end

% Initialize
vResults = [];
params.sComment = [];
if isempty(params.fBinning)
    params.fBinning = 0.1;
end

% Determine time period of catalog
params.fTminCat = min(params.mCatalog(:,3));
params.fTmaxCat = max(params.mCatalog(:,3));
% Adjust to decimal years
fTimePeriod =params.fTimePeriod/365;

% Init result matrix
mValueGrid_ = [];

% Temporary saving the original catalog
mCatalog = params.mCatalog;

% Force saving all 500 nodes
fDivide = length(params.mPolygon(:,1))/500
fForceSave = length(params.mPolygon(:,1))/fDivide;
% Check for bootstrapping or not
% ------------------------------
% Case of calculations with bootstrapping
if (params.bBstnum == 1)
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog(:,3) & params.mCatalog(:,3) < fTstart+fTimePeriod);
        params.mCatalog = params.mCatalog(vSel,:);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = nNodeStart:nNodeEnd
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
                mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN nX NaN NaN NaN NaN];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc_2method(params,mNodeCatalog_);
                % Store the results
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_Bst...
                        rCalcNodeResult_.fStd_Mc rCalcNodeResult_.fBvalue_Bst rCalcNodeResult_.fStd_B...
                        rCalcNodeResult_.fAvalue_Bst rCalcNodeResult_.fStd_A nX...
                        rCalcNodeResult_.bH_EMR rCalcNodeResult_.fPval rCalcNodeResult_.bH_Bst rCalcNodeResult_.fPval_Bst];
            end; % End of if on nNode_
            if rem(nNode_,floor(fForceSave)) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                params.vcsGridNames = cellstr(char('Mc max. curvature', 'Mc EMR-method', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
                    'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature', 'Mc EMR-method', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
            'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)...
                '_Nodes_' num2str(nNodeStart) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+fTimePeriod;
    end; % End of while fTstart

    % Case of no bootstrapping
else
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog(:,3) & params.mCatalog(:,3) < fTstart+fTimePeriod);
        params.mCatalog = params.mCatalog(vSel,:);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = nNodeStart:nNodeEnd
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
                mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN nX];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc(params,mNodeCatalog_);
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_90 rCalcNodeResult_.fMc_95 rCalcNodeResult_.fMc_com...
                        rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_shi nX];
            end; % End of if on nX
            if rem(nNode_,500) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)','Number of events'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['tmp_result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
            'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)...
                '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '_Nodes_' num2str(nNodeStart) '_' num2str(nNodeEnd) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+fTimePeriod;
    end; % End of while fTstart
end; % END of if params.bBst

author: CelsoReyes


text:
% ZMAP script show_map.m. Creates Dialog boxes for Z-map calculation
% does the calculation and makes displays the map
% stefan wiemer 11/94
%
% make dialog interface and call maxzlta
%
% This is the info window text
%

report_this_filefun(mfilename('fullpath'));

watchon
think
ttlStr='The Alarm Cube Window                                ';
hlpStr1= ...
    ['  To be implemented                             '
    ' corners with the mouse                         '];
% Find out of figure already exists
watchon
if ~exist('iala') ; iala = ZG.compare_window_yrs; end
if ~exist('abo2') ; errordlg('No alarms with z >= Zmin detected!');return; end
if isempty(abo2)  == 1 ; errordlg('No alarms with z >= Zmin detected!');return; end

abo = abo2;

[existFlag,figNumber]=figure_exists('Alarm Display',1);
newCubeWindowFlag=~existFlag;

% Set up the Seismicity Map window Enviroment
%
if newCubeWindowFlag
    cube = figure_w_normalized_uicontrolunits( ...
        'Name','Alarm Display',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'Visible','off', ...
        'Position',[  200 200 400 600]);

    matdraw
    ter2 = 7.5;
    tresh = max(loc(:,3));

    uicontrol('Units','normal',...
        'Position',[.0 .65 .12 .06],'String','Refresh ',...
         'Callback','close,plala_cu')

    tre2 = max(abo(:,4)) - 0.5;
    new = uicontrol('style','edit','value',ZG.compare_window_yrs,...
        'string',num2str(tre2,3), 'background','y',...
        'Callback','tre2=str2num(get(new,''String''));''String'',num2str(tre2,3);',...
        'units','norm','pos',[.80 .01 .08 .06],'min',2.65,'max',10);

    newlabel = uicontrol('style','text','units','norm','pos',[.40 .00 .40 .08]);
    set(newlabel,'string','Alarm Threshold:','background',color_fbg);

    mamo1 = uicontrol('Units','normal',...
        'Position',[.90 .01 .08 .06],'String','Go',...
         'Callback','abo = abo2;plala_cu ');

    mamo = uicontrol('Units','normal',...
        'Position',[.02 .01 .27 .10],'String','Make Movie',...
         'Callback','delete(mamo);delete(mamo1); delete(newlabel); mamovie ')

    nilabel2 = uicontrol('style','text','units','norm','pos',[.50 .92 .25 .06]);
    set(nilabel2,'string','MinRad (in km):','background',color_fbg);
    set_ni2 = uicontrol('style','edit','value',tresh,'string',num2str(tresh,3),...
        'background','y');
    set(set_ni2,'Callback','tresh=str2double(get(set_ni2,''String'')); set(set_ni2,''String'',num2str(tresh,3))');
    set(set_ni2,'units','norm','pos',[.80 .92 .13 .06],'min',0.01,'max',10000);

    uicontrol('Units','normal',...
        'Position',[.93 .93 .07 .05],'String','Go ',...
         'Callback','think;pause(1); plala_cu')

    op3 = uimenu('Label','Tools');
    uimenu(op3,'Label','Find Anomalie Groups  ',...
        'Callback','agroupc,cian');
    uimenu(op3,'Label','Display one Anomalie Group ',...
        'Callback','cian2');
    uimenu(op3,'Label','Determine Valarm/Vtotal(Zalarm) ',...
        'Callback','sucrac');
    uimenu(op3,'Label','Determine # Alarmgroups (Zalarm) ',...
        'Callback','agzc');

end   % if exist newCube

report_this_filefun(mfilename('fullpath'));

figure_w_normalized_uicontrolunits(cube)
delete(gca)
abo = abo2;
if isempty(abo)
    zmap_message_center.set_info(' ','No data above threshold');
    return;
end
rect= [0.2 0.2 0.6 0.6];
axes('pos',rect)
set(gca,'visible','off')
abo = abo2;
abo(:,5) = abo(:,5)* par1/365 + a(1,3);
l = abo(:,4) > tre2;
abo = abo(l,:);
if length(abo)  < 1
    errordlg('No alarms with z >= Zmin detected!');
    return;
end
l = abo(:,3) < tresh;
abo = abo(l,:);
if length(abo)  < 1  ; errordlg('No alarms with z >= Zmin detected!');
    return;
end
hold on

if isempty(abo) == 0
    figure_w_normalized_uicontrolunits(map)
    update(mainmap())
    plot(abo(:,1),abo(:,2),'o',...
        'MarkerFaceColor','r','MarkerEdgeColor','y');

    figure_w_normalized_uicontrolunits(cube)
    plo  = plot3(abo(:,1),abo(:,2),abo(:,5),'ro');
    set(plo,'MarkerSize',6,'LineWidth',1.0)
    for i = 1:length(abo(:,1))
        li = [abo(i,1) abo(i,2) abo(i,5) ; abo(i,1) abo(i,2) abo(i,5)+iala];
        plot3(li(:,1),li(:,2),li(:,3),'b');
    end
end
view(3);

grid
hold on
%end

strib4 = [  ' Alarm Cube of '  name '; wl =  '  num2str(ZG.compare_window_yrs,3) '; Zcut = ' num2str(tre2,3)  ];
title2(strib4,'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',2.0,'visible','on')
%set(gca,'Color',[0.7 0.7 0.7])

viewer
watchoff
vie = gcf;
figure_w_normalized_uicontrolunits(cube)

watchoff
done;
figure_w_normalized_uicontrolunits(cube)

rotate3d on


text:
% mff_exportcategories - export MFF EEG event from EEGLAB structure to
%                    'categories.xml' file
% Usage:
%   mff_exportcategories(EEG, mffFile);
%
% Inputs:
%  EEG     - EEGLAB structure
%  mffFile - filename/foldername for the MFF file (MFF file/folder must
%            already exist)

% This file is part of mffmatlabio.
%
% mffmatlabio is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% mffmatlabio is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with mffmatlabio.  If not, see <https://www.gnu.org/licenses/>.

function indtle = mff_exportcategories(EEG, mffFile)

indtle = []; % time locking event indices
if EEG.trials == 1
    return;
end

p = fileparts(which('mff_importsignal.m'));
warning('off', 'MATLAB:Java:DuplicateClass');
javaaddpath(fullfile(p, 'MFF-1.2.2-jar-with-dependencies.jar'));
warning('on', 'MATLAB:Java:DuplicateClass');

% Create a factory.
mfffactorydelegate = javaObject('com.egi.services.mff.api.LocalMFFFactoryDelegate');
mfffactory         = javaObject('com.egi.services.mff.api.MFFFactory', mfffactorydelegate);

categoriesRType = javaObject('com.egi.services.mff.api.MFFResourceType', javaMethod('valueOf', 'com.egi.services.mff.api.MFFResourceType$MFFResourceTypes', 'kMFF_RT_Categories'));
if mfffactory.createResourceAtURI(fullfile(mffFile, 'categories.xml'), categoriesRType)
    fprintf('categories.xml file created successfully\n');
else
    fprintf('categories.xml ressource already exist, overwriting\n');
end
catsResource = mfffactory.openResourceAtURI( fullfile(mffFile, 'categories.xml'), categoriesRType);

% find time-locking events
% ------------------------
events   = EEG.event;
eventlat = abs(eeg_point2lat( [ events.latency ], [ events.epoch ], EEG.srate, [EEG.xmin EEG.xmax]));
indtle    = find(eventlat == 0);
if length(indtle) < EEG.trials
    indtle    = find(eventlat < 0.02);
    if length(indtle) ~= EEG.trials
        if isfield(EEG.event, 'begintime')
            indtle = find(cellfun(@isempty, { EEG.event.begintime }));
            if length(indtle) ~= EEG.trials
                disp('WARNING: NOT THE RIGHT NUMBER OF TIME-LOCKING EVENTS');
            end
        end
    end
end

events = events(indtle);
uniqueType = unique( { events.type } );
jCatList = javaObject('java.util.ArrayList');
epochLen = [];
for iType = 1:length(uniqueType)
    trials = strmatch( uniqueType{iType}, { events.type }, 'exact'); % get time locking events

    catObj = javaObject('com.egi.services.mff.api.Category');
    catObj.setName( events(trials(1)).type );
    jList = javaObject('java.util.ArrayList');
    for iTrial = 1:length(trials)
        segmentObj = javaObject('com.egi.services.mff.api.Segment');
        epoch      = events(trials(iTrial)).epoch;

        segmentObj.setName( events(trials(iTrial)).type );
        segmentObj.setBeginTime(  round((epoch-1)*EEG.pnts/round(EEG.srate)*1000000) );
        segmentObj.setEndTime(    round((epoch  )*EEG.pnts/round(EEG.srate)*1000000) );
        segmentObj.setEventBegin( round((events(trials(iTrial)).latency-1)/round(EEG.srate)*1000000) );
        if isfield(events(trials(iTrial)), 'duration') && ~isempty(events(trials(iTrial)).duration)
             duration = events(trials(iTrial)).duration;
        else duration = 0;
        end
        segmentObj.setEventEnd(   round((events(trials(iTrial)).latency-1+duration)/round(EEG.srate)*1000000) );
        if isfield(events, 'status')
            segmentObj.setStatus(  events(trials(iTrial)).status );
        end
        if isfield(events, 'mffkeysbackup')
            % if additional keys are present -> average (not sure this heuristic is always correct)
            segmentObj.setName( 'Average' );

            if ~isempty(events(trials(iTrial)).mffkeysbackup)
                jListKeys = javaObject('java.util.ArrayList');
                tmpKeys = eval(events(trials(iTrial)).mffkeysbackup);
                for iKey = 1:length(tmpKeys);
                    keyObj = javaObject('com.egi.services.mff.api.Key');
                    keyObj.setCode(tmpKeys(iKey).code);
                    keyObj.setData(tmpKeys(iKey).data);
                    keyObj.setDataType(tmpKeys(iKey).datatype);
                    keyObj.setDescription(tmpKeys(iKey).description);
                    jListKeys.add(keyObj);
                end
                segmentObj.setKeys(jListKeys);
            end
        end

        % check epoch length
        epochLenTmp = segmentObj.getEndTime() - segmentObj.getBeginTime();
        if isempty(epochLen)
            epochLen = epochLenTmp;
        elseif epochLen ~= epochLenTmp
            error('Epoch length export error')
        end

        jList.add(segmentObj);
    end
    catObj.setSegments(jList);
    jCatList.add(catObj);
end

catsResource.setCategories(jCatList);
catsResource.saveResource();

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Fig. 5.38  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
clf
numG = 160*conv ([1 2.5],[1 0.7]);
denG = conv([1 5 40],[1  .03  .06]);
sysG = tf(numG,denG);
sysD=tf([1   3],[1   20]);
sysDG=sysD*sysG;
K = 0.3;
sysH=tf(1,1);
sysT = feedback (K*sysG,sysH);
sysTD=feedback(sysDG,sysH);
step(sysT)
hold on
step(sysTD)
grid on
title('Step responses of auto-pilot with P and lead control')
hold off


text:
function [Ain, Cin,  bin, fin, center, res] = greedyROI_corr(Y, K, options, sn, debug_on, save_avi)
%% a greedy method for detecting ROIs and initializing CNMF. in each iteration,
% it searches the one with large (peak-median)/noise level and large local
% correlation
%% Input:
%   Y:  d X T matrx, imaging data
%   K:  scalar, maximum number of neurons to be detected.
%   options: struct data of paramters/options
%       d1:     number of rows
%       d2:     number of columns
%       gSiz:   maximum size of a neuron
%       nb:     number of background
%       min_corr: minimum threshold of correlation for segementing neurons
%   sn:     d X 1 vector, noise level of each pixel
%   debug_on: options for showing procedure of detecting neurons
%% Output:
%       Ain:  d X K' matrix, estimated spatial component
%       Cin:  K'X T matrix, estimated temporal component
%       bin:  d X nb matrix/vector, spatial components of the background
%       Cin:  nb X T matrix/vector, temporal components of the background
%       center: K' X 2, coordinate of each neuron's center
%       res:  d X T, residual after initializing Ain, Cin, bin, fin

%% Author: Pengcheng Zhou, Carnegie Mellon University.
% the method is an modification of greedyROI method used in Neuron paper of Eftychios
% Pnevmatikakis et.al. https://github.com/epnev/ca_source_extraction/blob/master/utilities/greedyROI2d.m
%% In each iteration of initializing neurons, it searchs the one with maximum
% value of (max-median)/noise * Cn, which selects pixels with high SNR and
% local correlation.

%% parameters
if exist('sn', 'var')&& ~(isempty(sn))
    Y_std = sn;
else
    Y_std = std(Y, 0, ndims(Y));
end
Y_std = Y_std(:);
if ~exist('debug_on', 'var'); debug_on = false; end
if ~exist('save_avi', 'var'); save_avi=false; end
d1 = options.d1;
d2 = options.d2;
gSig = options.gSig;
gSiz = options.gSiz;
if and(isempty(gSiz), isempty(gSig)); gSig = 3; gSiz = 10; end
if isempty(gSiz); gSiz=3*gSig; end
if isempty(gSig); gSig=gSiz/3; end
min_corr = options.min_corr;    % minimum local correaltion value to start one neuron
nb = options.nb;        % number of the background
pSiz = 1;       % after selecting one pixel, take the mean of square box
%near the pixel as temporal activity. the box size is (2*pSiz+1)
psf = ones(gSig)/(gSig^2);

min_snr = 3;        % minimum value of (peak-median)/sig
maxIter = 5;            % iterations for refining results
sz = 4;            %distance of neighbouring pixels for computing local correlation

if ~ismatrix(Y); Y = reshape(Y, d1*d2, []); end;
[~, T] = size(Y);       % number of frames
Ain = zeros(d1*d2, K);  % spatial components
Cin = zeros(K, T);      % temporal components
center = zeros(K, 2);   % center of the initialized components

%% compute correlation image and (max-median)/std ratio
ind_frame = round(linspace(1, T, min(T, 1000)));    % select few frames for speed issues
%tmp_noise = randn(d1*d2, length(ind_frame));
C1 = correlation_image(full(Y(:, ind_frame)), sz, d1, d2);
Cb =  zeros(size(C1)); %correlation_image(full(Y(:, ind_frame(1:3:end)))+tmp_noise(:, 1:3:end), [gSiz, gSiz+1], d1, d2);  %backgroung correlatin
Cn = C1-Cb; %here Cb is the background correlation. for 2photon imaging results. It might be useful when the background signal is large
Y_median = median(Y(:, ind_frame), 2);
Y = bsxfun(@minus, Y, Y_median);
% Y_std = sqrt(mean(Y.*Y, 2));

%% find local maximum
k = 0;      %number of found components
min_pixel = floor(gSig^2/2);  % minimum number of peaks to be a neuron
peak_ratio = full(max(Y, [], 2))./Y_std; %(max-median)/std
peak_ratio(isinf(peak_ratio)) = 0;  % avoid constant values
% save_avi = false;   %save procedures for demo
if debug_on
    figure('position', [100, 100, 800, 650]); %#ok<*UNRCH>
    subplot(331);
    imagesc(Cn, [0,1]); colorbar;
    axis equal off tight; hold on;
    title('correlation image');
    if save_avi
        avi_file = VideoWriter('greedyROI_example.avi');
        avi_file.open();
    end
end

max_thresh = min_snr * (min_corr);
while k<K
    %% find the pixel with the maximum ratio
    [max_v, ind_p] = max(peak_ratio.*(Cn(:)));
    peak_ratio(ind_p) = 0;  % no longer visit this pixel any more
    if max_v<max_thresh; break; end
    if Cn(ind_p)<min_corr; continue; end % ignore this local maximum due to small local correlation
    if  max_v/(min_corr)< min_snr;     continue;    end

    [r, c] = ind2sub([d1,d2], ind_p);

    % select its neighbours for computing correlation
    rsub = max(1, -gSiz+r):min(d1, gSiz+r);
    csub = max(1, -gSiz+c):min(d2, gSiz+c);
    [cind, rind] = meshgrid(csub, rsub);
    nr = length(rsub);  %size of the neighboring matrix
    nc = length(csub);
    ind_nhood = sub2ind([d1, d2], rind(:), cind(:));
    Y_box = Y(ind_nhood, :);

    % draw a small area near the peak and extract the mean activities
    r0 = rsub(1); c0 = csub(1);
    rsub = (max(1, -pSiz+r):min(d1, pSiz+r)) - r0+1;
    csub = (max(1, -pSiz+c):min(d2, pSiz+c)) -c0+1;
    [cind, rind] = meshgrid(csub, rsub);
    ind_peak = sub2ind([nr, nc], rind(:), cind(:));
    y0 = mean(Y_box(ind_peak, :), 1);
    y0(y0<0) = 0;

    % compute the correlation between the peak and its neighbours
    temp = reshape(corr(y0', Y_box'), nr, nc);
    active_pixel = full(temp>min_corr/2);
    l = bwlabel(active_pixel, 8);   % remove disconnected components
    active_pixel(l~=mode(l(ind_peak))) = false;
    tmp_v = sum(active_pixel(:));    %number of pixels with above-threshold correlation
    if debug_on
        subplot(332); cla;
        imagesc(reshape(peak_ratio.*Cn(:), d1, d2), [0, max_v]); colorbar;
        title(sprintf('neuron %d', k+1));
        axis equal off tight; hold on;
        plot(c,r, 'om');
        subplot(333);
        imagesc(temp, [min_corr, 1]);
        axis equal off tight;
        title('corr. with neighbours');
        subplot(3,3,4:6); cla;
        plot(y0); title('activity in the center');
        subplot(3,3,7:9); cla;
        if ~save_avi; pause; end
    end
    if tmp_v<min_pixel;         continue;  end % neuron is too small

    %% save neuron
    %   nonzero area
    data = Y_box(active_pixel(:), :);
    ind_active = ind_nhood(active_pixel(:));  %indices of active pixels within the whole frame
    peak_ratio(ind_nhood(ind_peak)) = 0;    % the small area near the peak is not able to initialize neuron anymore

    % do a rank-1 matrix factorization in this small area
    [ai, ci] = finetune2d(data, y0, maxIter);
    %     data(data<0) = 0;
    %     ai = (data*y0')/(y0*y0');
    %     ci = y0;
    if norm(ai)==0;        continue;   end
    k = k+1;
    Ain(ind_active, k) = ai;
    Cin(k, :) = ci;
    Y(ind_active, :) = data-ai*ci;
    center(k, :) = [r, c];

    if debug_on
        subplot(331);
        plot(c, r, '.r');
        subplot(332);
        plot(c,r, 'or');
        subplot(333);
        temp = zeros(nr, nc); temp(active_pixel) = ai;
        imagesc(temp);
        axis equal off tight;
        title('spatial component');
        subplot(3,3,7:9); cla;
        plot(ci); title('temporal component');
        if save_avi;
            temp = getframe(gcf);
            temp.cdata = imresize(temp.cdata, [800, 640]);
            avi_file.writeVideo(temp);
        else pause; end
    end

    if mod(k, 10)==0
        fprintf('%d/%d neurons have been detected\n', k, K);
    end

    if k==K;   break; end

    %% udpate peak_ratio and correlation image
    tmp_old = peak_ratio(ind_active);
    tmp_new = max(Y(ind_active, :), [], 2)./Y_std(ind_active);
    temp = zeros(nr, nc);
    temp(active_pixel) = max(0, tmp_old-tmp_new); % after each iteration, the peak ratio can not be increased
    peak_ratio(ind_nhood) = max(0, peak_ratio(ind_nhood) - reshape(imfilter(temp, psf), [], 1)); % update peak_ratio, results are smoothed
    Cn(ind_nhood) = correlation_image(full(Y(ind_nhood, ind_frame)), sz, nr, nc)-reshape(Cb(ind_nhood), nr, nc);  % update local correlation
end

center = center(1:k, :);
Ain = sparse(Ain(:, 1:k));
Cin = Cin(1:k, :);
Cin(Cin<0) = 0;
if save_avi; avi_file.close(); end
res = bsxfun(@plus, Y, Y_median);

%% initialize background
tsub = max(1, round(T/1000));
[bin, f] = nnmf(max(res(:, 1:tsub:T), 0), nb);
fin = imresize(f, [nb, T]);
fin = HALS_temporal(max(res, 0), bin, fin, maxIter);
bin = HALS_spatial(max(res, 0), bin, fin, [], maxIter);
end

function [ai, ci] = finetune2d(data, ci, nIter)
%do matrix factorization given the model data = ai*ci, where ai>=0
%
%Input:
%   data:   d x T matrix, small patch containing one neuron
%   ci:     initial value for trace
%   nIter  number of coordinate descent steps
%
%Output:
%   ai  M x N matrix, result of the fine-tuned neuron shape
%   ci  1 x T matrix, result of the neuron
%% copied from greedyROI.m

if ~exist('nIter', 'var'), nIter = 1; end
data(data<0)= 0;
%do block coordinate descent
for iter = 1:nIter,
    %update basis
    ai = max(0, (data*ci')/(ci*ci'));
    norm_ai = norm(ai, 2);
    if norm_ai==0; break;     end
    ai = ai/norm_ai;
    ci =  (ai'*data);
    %     ci(ci<0) = 0;
end
temp = (median(ci)-2*std(ci));
ci(ci<temp) = temp;
end

author: zhoupc