text:
function gpcf = gpcf_constant(varargin)
%GPCF_CONSTANT  Create a constant covariance function
%
%  Description
%    GPCF = GPCF_CONSTANT('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a squared exponential covariance function structure in
%    which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    GPCF = GPCF_CONSTANT(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for constant covariance function [default]
%      constSigma2       - magnitude (squared) [0.1]
%      constSigma2_prior - prior for constSigma2 [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010      Jaakko Riihimaki, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPCF_CONSTANT';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('constSigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('constSigma2_prior',prior_logunif, @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    init=true;
    gpcf.type = 'gpcf_constant';
  else
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_constant')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  % Initialize parameter
  if init || ~ismember('constSigma2',ip.UsingDefaults)
    gpcf.constSigma2=ip.Results.constSigma2;
  end

  % Initialize prior structure
  if init
    gpcf.p=[];
  end
  if init || ~ismember('constSigma2_prior',ip.UsingDefaults)
    gpcf.p.constSigma2=ip.Results.constSigma2_prior;
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_constant_pak;
    gpcf.fh.unpak = @gpcf_constant_unpak;
    gpcf.fh.lp = @gpcf_constant_lp;
    gpcf.fh.lpg = @gpcf_constant_lpg;
    gpcf.fh.cfg = @gpcf_constant_cfg;
    gpcf.fh.ginput = @gpcf_constant_ginput;
    gpcf.fh.cov = @gpcf_constant_cov;
    gpcf.fh.trcov  = @gpcf_constant_trcov;
    gpcf.fh.trvar  = @gpcf_constant_trvar;
    gpcf.fh.recappend = @gpcf_constant_recappend;
  end

end

function [w, s] = gpcf_constant_pak(gpcf, w)
%GPCF_CONSTANT_PAK  Combine GP covariance function parameters into
%                   one vector.
%
%  Description
%    W = GPCF_CONSTANT_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for example
%    in energy and gradient computations.
%
%       w = [ log(gpcf.constSigma2)
%             (hyperparameters of gpcf.constSigma2)]'
%
%  See also
%    GPCF_CONSTANT_UNPAK

  w = []; s = {};

  if ~isempty(gpcf.p.constSigma2)
    w = log(gpcf.constSigma2);
    s = [s 'log(constant.constSigma2)'];
    % Hyperparameters of constSigma2
    [wh sh] = gpcf.p.constSigma2.fh.pak(gpcf.p.constSigma2);
    w = [w wh];
    s = [s sh];
  end
end

function [gpcf, w] = gpcf_constant_unpak(gpcf, w)
%GPCF_CONSTANT_UNPAK  Sets the covariance function parameters
%                     into the structure
%
%  Description
%    [GPCF, W] = GPCF_CONSTANT_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance parameters have been set
%    to the values in W. Deletes the values set to GPCF from W
%    and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(gpcf.constSigma2)
%             (hyperparameters of gpcf.constSigma2)]'
%
%  See also
%   GPCF_CONSTANT_PAK

  gpp=gpcf.p;
  if ~isempty(gpp.constSigma2)
    gpcf.constSigma2 = exp(w(1));
    w = w(2:end);
    % Hyperparameters of magnSigma2
    [p, w] = gpcf.p.constSigma2.fh.unpak(gpcf.p.constSigma2, w);
    gpcf.p.constSigma2 = p;
  end
end

function lp = gpcf_constant_lp(gpcf)
%GPCF_CONSTANT_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_CONSTANT_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LPG, GP_E

% Evaluate the prior contribution to the error. The parameters that
% are sampled are from space W = log(w) where w is all the
% "real" samples. On the other hand errors are evaluated in the
% W-space so we need take into account also the Jacobian of
% transformation W -> w = exp(W). See Gelman et.al., 2004,
% Bayesian data Analysis, second edition, p24.

  lp = 0;
  gpp=gpcf.p;
  if ~isempty(gpp.constSigma2)
    lp = gpp.constSigma2.fh.lp(gpcf.constSigma2, gpp.constSigma2) +log(gpcf.constSigma2);
  end
end

function lpg = gpcf_constant_lpg(gpcf)
%GPCF_CONSTANT_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_CONSTANT_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  lpg = [];
  gpp=gpcf.p;

  if ~isempty(gpcf.p.constSigma2)
    lpgs = gpp.constSigma2.fh.lpg(gpcf.constSigma2, gpp.constSigma2);
    lpg = [lpg lpgs(1).*gpcf.constSigma2+1 lpgs(2:end)];
  end
end

function DKff = gpcf_constant_cfg(gpcf, x, x2, mask, i1)
%GPCF_CONSTANT_CFG  Evaluate gradient of covariance function
%                   with respect to the parameters
%
%  Description
%    DKff = GPCF_CONSTANT_CFG(GPCF, X) takes a
%    covariance function structure GPCF, a matrix X of input
%    vectors and returns DKff, the gradients of covariance matrix
%    Kff = k(X,X) with respect to th (cell array with matrix
%    elements). This is a mandatory subfunction used in gradient
%    computations.
%
%    DKff = GPCF_CONSTANT_CFG(GPCF, X, X2) takes a
%    covariance function structure GPCF, a matrix X of input
%    vectors and returns DKff, the gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_CONSTANT_CFG(GPCF, X, [], MASK)
%    takes a covariance function structure GPCF, a matrix X of
%    input vectors and returns DKff, the diagonal of gradients of
%    covariance matrix Kff = k(X,X2) with respect to th (cell
%    array with matrix elements). This subfunction is needed when
%    using sparse approximations (e.g. FIC).
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  [n, m] =size(x);

  DKff = {};

  if nargin==5
    % Use memory save option
    if i1==0
      % Return number of hyperparameters
      if ~isempty(gpcf.p.constSigma2)
        DKff=1;
      else
        DKff=0;
      end
      return
    end
  end

  % Evaluate: DKff{1} = d Kff / d constSigma2
  %           DKff{2} = d Kff / d coeffSigma2
  % NOTE! Here we have already taken into account that the parameters are transformed
  % through log() and thus dK/dlog(p) = p * dK/dp

  % evaluate the gradient for training covariance
  if nargin == 2 || (isempty(x2) && isempty(mask))

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones(n)*gpcf.constSigma2;
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_constant -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones([n size(x2,1)])*gpcf.constSigma2;
    end

    % Evaluate: DKff{1}    = d mask(Kff,I) / d constSigma2
    %           DKff{2...} = d mask(Kff,I) / d coeffSigma2
  elseif nargin == 4 || nargin == 5

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones(n,1)*gpcf.constSigma2; % d mask(Kff,I) / d constSigma2
    end
  end
  if nargin==5
    DKff=DKff{1};
  end

end

function DKff = gpcf_constant_ginput(gpcf, x, x2, i1)
%GPCF_CONSTANT_GINPUT  Evaluate gradient of covariance function with
%                      respect to x.
%
%  Description
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  [n, m] =size(x);
  if nargin==4
    % Use memory save option
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  end

  if nargin == 2 || isempty(x2)
    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n);
      end
    end

  elseif nargin == 3 || nargin == 4
    %K = feval(gpcf.fh.cov, gpcf, x, x2);

    ii1 = 0;
    for i=1:m
      for j = 1:n
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n, size(x2,1));
        gprior(ii1) = 0;
      end
    end
  end
  if nargin==5
    DKff=DKff{1};
  end
end

function C = gpcf_constant_cov(gpcf, x1, x2, varargin)
%GP_CONSTANT_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_CONSTANT_COV(GP, TX, X) takes in covariance function
%    of a Gaussian process GP and two matrixes TX and X that
%    contain input vectors to GP. Returns covariance matrix C.
%    Every element ij of C contains covariance between inputs i
%    in TX and j in X. This is a mandatory subfunction used for
%    example in prediction and energy computations.
%
%  See also
%    GPCF_CONSTANT_TRCOV, GPCF_CONSTANT_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C = ones(n1,n2)*gpcf.constSigma2;
end

function C = gpcf_constant_trcov(gpcf, x)
%GP_CONSTANT_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_CONSTANT_TRCOV(GP, TX) takes in covariance function
%    of a Gaussian process GP and matrix TX that contains
%    training input vectors. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i and j
%    in TX. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_CONSTANT_COV, GPCF_CONSTANT_TRVAR, GP_COV, GP_TRCOV

  n =size(x,1);
  C = ones(n,n)*gpcf.constSigma2;

end

function C = gpcf_constant_trvar(gpcf, x)
%GP_CONSTANT_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_CONSTANT_TRVAR(GPCF, TX) takes in covariance function
%    of a Gaussian process GPCF and matrix TX that contains
%    training inputs. Returns variance vector C. Every
%    element i of C contains variance of input i in TX. This is
%    a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_CONSTANT_COV, GP_COV, GP_TRCOV

  n =size(x,1);
  C = ones(n,1)*gpcf.constSigma2;

end

function reccf = gpcf_constant_recappend(reccf, ri, gpcf)
%RECAPPEND Record append
%
%  Description
%    RECCF = GPCF_CONSTANT_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains all
%    the old samples and the current samples from GPCF. This
%    subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_constant';

    % Initialize parameters
    reccf.constSigma2 = [];

    % Set the function handles
    reccf.fh.pak = @gpcf_constant_pak;
    reccf.fh.unpak = @gpcf_constant_unpak;
    reccf.fh.lp = @gpcf_constant_lp;
    reccf.fh.lpg = @gpcf_constant_lpg;
    reccf.fh.cfg = @gpcf_constant_cfg;
    reccf.fh.cov = @gpcf_constant_cov;
    reccf.fh.trcov  = @gpcf_constant_trcov;
    reccf.fh.trvar  = @gpcf_constant_trvar;
    reccf.fh.recappend = @gpcf_constant_recappend;
    reccf.p=[];
    reccf.p.constSigma2=[];
    if ~isempty(ri.p.constSigma2)
      reccf.p.constSigma2 = ri.p.constSigma2;
    end

  else
    % Append to the record
    gpp = gpcf.p;

    % record constSigma2
    reccf.constSigma2(ri,:)=gpcf.constSigma2;
    if isfield(gpp,'constSigma2') && ~isempty(gpp.constSigma2)
      reccf.p.constSigma2 = gpp.constSigma2.fh.recappend(reccf.p.constSigma2, ri, gpcf.p.constSigma2);
    end
  end
end

author: fieldtrip

repo: fieldtrip


text:
function hash = wsddnVOChash_init(strs)
% From the PASCAL VOC 2011 devkit

hsize=4999;
hash.key=cell(hsize,1);
hash.val=cell(hsize,1);

for i=1:numel(strs)
    s=strs{i};
    h=mod(str2double(s([4 6:end])),hsize)+1;
    j=numel(hash.key{h})+1;
    hash.key{h}{j}=strs{i};
    hash.val{h}(j)=i;
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function out = mapFeature(X1, X2)
% MAPFEATURE Feature mapping function to polynomial features
%
%   MAPFEATURE(X1, X2) maps the two input features
%   to quadratic features used in the regularization exercise.
%
%   Returns a new feature array with more features, comprising of
%   X1, X2, X1.^2, X2.^2, X1*X2, X1*X2.^2, etc..
%
%   Inputs X1, X2 must be the same size
%

degree = 6;
out = ones(size(X1(:,1)));
for i = 1:degree
    for j = 0:i
        out(:, end+1) = (X1.^(i-j)).*(X2.^j);
    end
end

end

repo: MachineLearningEx


text:
function test_failed=test_gdgt
%TEST_GDGT  Test GDGT
%
%  This script runs a throrough test of the COMP_GDGT routine, testing it on
%  a range of input parameters.
%

Lr=[24,144,108,144,24,135,35,77,20];
ar=[ 4,  9,  9, 12, 6,  9, 5, 7, 1];
Mr=[ 6, 16, 12, 24, 8,  9, 7,11,20];

R=1;

test_failed=0;

disp(' ===============  TEST_GDGT ================');

for ii=1:length(Lr);

  for ctii=0:1

    c_t=ctii*.5;

    for cfii=0:1

      c_f=cfii*.5;

      for W=1:3

	L=Lr(ii);

	M=Mr(ii);
	a=ar(ii);

	b=L/M;
	N=L/a;
	c=gcd(a,M);
	d=gcd(b,N);
	p=a/c;
	q=M/c;

	%g=(1:L)';
	%g=i*gabtight(a,b,L);
	f=tester_crand(L,W);
	g=tester_crand(L,R);

	gd=gabdual(g,a,M);
	gt=gabtight(g,a,M);

	cc=comp_gdgt(f,g,a,M,L,c_t,c_f,0,0);
	cc2=reshape(ref_gdgt(f,g,a,M,c_t,c_f,0),M,N,W);

	cdiff=cc-cc2;
	res=norm(cdiff(:));
        [test_failed,fail]=ltfatdiditfail(res,test_failed);
	s=sprintf(['REF L:%3i c_t: %0.5g c_f %0.5g ', ...
                   'W:%2i R:%2i a:%3i b:%3i c:%3i ', ...
                   'd:%3i p:%3i q:%3i %0.5g %s'],...
                  L,c_t,c_f,W,R,a,b,c,d,p,q,res,fail);
	disp(s)

	r=comp_igdgt(cc,gd,a,M,L,c_t,c_f,0,0);
	res=norm(f-r,'fro');
	[test_failed,fail]=ltfatdiditfail(res,test_failed);

	s=sprintf(['REC L:%3i c_t: %0.5g c_f %0.5g ',...
                   'W:%2i R:%2i a:%3i b:%3i c:%3i ',...
                   'd:%3i p:%3i q:%3i %0.5g %s'],...
                  L,c_t,c_f,W,R,a,b,c,d,p,q,res,fail);
	disp(s)

	res=norm(f-idgt(dgt(f,gt,a,M),gt,a),'fro');
        [test_failed,fail]=ltfatdiditfail(res,test_failed);
	s=sprintf(['TIG L:%3i c_t: %0.5g c_f %0.5g ',...
                   'W:%2i R:%2i a:%3i b:%3i c:%3i ',...
                   'd:%3i p:%3i q:%3i %0.5g %s'],...
                  L,c_t,c_f,W,R,a,b,c,d,p,q,res,fail);
	disp(s);
      end;

    end;
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function bc_start(mCatalog, hFigure, bMap, vCoastline, vFaults, rContainer, sName)
    % Starts the probabilistic forecast testing.
    %
    % bc_start(mCatalog, hFigure, bMap, rContainer)
    % ------------------------------------------------------
    % Starts the probabilistic forecast testing. The function invokes the parameter dialog, calls the
    %   calculation function and invokes the dialog for displaying the results
    %
    % Input parameters:
    %   mCatalog      Earthquake catalog to use for the testing
    %   hFigure       Handle of figure where the user should select the grid (i.e. the seismicity map)
    %   bMap          Map/cross-section switch. If the testing is carried out on a map set bMap = true,
    %                 on a cross-section set bMap = false
    %   rContainer    Just a container (structure proposed) to store any variables that should be available in the
    %                 probabilistic forecast testing code
    %
    % Danijel Schorlemmer
    % June 26, 2003

    report_this_filefun();

    % Launch GUI
    hMenuFig = bc_options(bMap);

    % Set up parameter struct and store input parameters
    params.mCatalog = mCatalog;
    params.sCatalogName = sName;
    params.bMap = bMap;
    params.vFaults = vFaults;
    params.vCoastline = vCoastline;
    params.rContainer = rContainer;

    % Analyze Output
    if ~ishandle(hMenuFig)
        answer = 0;
    else
        handles = guidata(hMenuFig);
        answer = handles.answer;
        % OK pressed
        if answer == 1
            % Get parameter values from the gui
            if get(handles.radNumber, 'Value') == 1
                rOptions.nGriddingMode = 0;
                rOptions.nNumberEvents = str2double(get(handles.txtNumber, 'String'));
                rOptions.fRadius = 0;
                rOptions.fSizeRectX = 0;
                rOptions.fSizeRectY = 0;
                rOptions.nMinimumNumber = 0;
                rOptions.fMaximumRadius = str2double(get(handles.txtMaximumRadius, 'String'));
            elseif get(handles.radRadius, 'Value') == 1
                rOptions.nGriddingMode = 1;
                rOptions.nNumberEvents = 0;
                rOptions.fRadius = str2double(get(handles.txtRadius, 'String'));
                rOptions.fSizeRectX = 0;
                rOptions.fSizeRectY = 0;
                rOptions.nMinimumNumber = str2double(get(handles.txtMinimumNumber, 'String'));
                rOptions.fMaximumRadius = 0;
            else % rectangle
                rOptions.nGriddingMode = 2;
                rOptions.nNumberEvents = 0;
                rOptions.fRadius = 0;
                rOptions.fSizeRectX = str2double(get(handles.txtSizeRectX, 'String'));
                rOptions.fSizeRectY = str2double(get(handles.txtSizeRectY, 'String'));
                rOptions.nMinimumNumber = str2double(get(handles.txtMinimumNumber, 'String'));
                rOptions.fMaximumRadius = 0;
            end
            rOptions.nCalculateMC = get(handles.cboCalculateMC, 'Value');
            rOptions.bConstrainMc = get(handles.chkConstrainMc, 'Value');
            rOptions.fMcMin = str2double(get(handles.txtMcMin, 'String'));
            rOptions.fMcMax = str2double(get(handles.txtMcMax, 'String'));
            if get(handles.cboBinning, 'Value') == 1
                rOptions.fBinning = 0.1;
            else
                rOptions.fBinning = 0.01;
            end

            params.rOptions = rOptions;

            % Options section
            params.bGridEntireArea = get(handles.chkGridEntireArea, 'Value');
            params.fSpacingX = str2double(get(handles.txtSpacingX, 'String'));
            params.fSpacingY = str2double(get(handles.txtSpacingY, 'String'));
            params.bUtsuTest = get(handles.chkUtsuTest, 'Value');
            params.fGamma = str2double(get(handles.txtGamma, 'String'));
            bSaveParameter = get(handles.chkSaveParameter, 'Value');
            if bSaveParameter
                sSaveParameter = get(handles.lblSaveParameter, 'String');
            end
            params.sComment = 'BeCubed';

            % Remove figure from memory
            delete(hMenuFig);

            % Select grid
            [params.mPolygon, params.vX, params.vY, params.vUsedNodes] = ex_selectgrid(hFigure, params.fSpacingX, params.fSpacingY, params.bGridEntireArea);

            % Validate polygonsize
            if length(params.vX) < 4  ||  length(params.vY) < 4
                errordlg('Selection is too small. Please select a larger polygon.');
                return;
            end

            % Create Indices to catalog
            params.rOptions.caNodeIndices = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, ...
                params.rOptions.nGriddingMode, params.rOptions.nNumberEvents, params.rOptions.fRadius, ...
                params.rOptions.fSizeRectX, params.rOptions.fSizeRectY);

            if bSaveParameter
                % Save parameters
                save(sSaveParameter, 'params');
            else
                % Perform the calculation
                [params] = bc_calc(params);
                % Display the results
                gui_result(params);
            end
        else
            delete(hMenuFig);
        end
    end
end

author: CelsoReyes


text:
%########################################################################
%
%	- PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	separate.m:
%
% Description:
%
%   Implements the Kalman based separation of frequencies based upon the diffusion process algorithm.
%
% References:
%
%	Christian S. Pilz, Jarek Krajewski, Vladimir Blazek.
%   On the Diffusion Process for Heart Rate Estimation from Face Videos under Realistic Conditions.
%   Pattern Recognition: 39th German Conference, GCPR 2017, Basel, Switzerland.
%   Proceedings (Lecture Notes in Computer Science), Springer, 2017
%

function [S,CS,QCS,SMM,SPP,FMM,FPP] = separate(Y,dt,FF,nharm,BF,BQ,BL,BH,R,qr)
% - Separate periodic signal from other signals
%
% Syntax:
%   [S,CS,QCS,SMM,SPP,FMM,FPP] = separate(Y,dt,FF,nharm,BF,BQ,BL,BH,R,qr)
%
% In:
%   Y  - Signal as 1xD vector
%   dt - Sampling period in seconds (e.g. 0.1)
%   FF - Fundamental frequencies as 1xD vector in BPM (from the IMM)
%   nharm - Number of harmonics to be estimated (e.g. 3)
%   BF - Feedback matrix for the bias model (e.g. [0 1; 0 0])
%   BQ - Process spectral density for the bias model (e.g. 0.01)
%   BL - Noise multiplier matrix for the bias model (e.g. [0;1])
%   BH - Measurement matrix for the bias model (e.g. [1 0])
%   R  - Measurement variance (e.g. 0.1^2)
%   qr - Resonator's process noise spectral density (e.g. 0.1)
%
% Out:
%   S   - Cleaned signal as 1xD vector
%   CS  - Fundamental signal and its harmonics as (nharm)xD vector
%   QCS - Quadrature periodic signals as (nharm)xD vector
%   SMM - The smoother means NxD
%   SPP - The smoother covariances NxNxD
%   FMM - The filter means NxD
%   FPP - The filter covariances NxNxD
%
% Description:
%   Separate given signal into true signal and periodic
%   signal by using the known frequency of the periodic signal.
%   See TRACK for the model used.

    % Number of harmonics (including fundamental)
    N = nharm;

    % Prior means and covariances
    m0 = zeros(2*N+size(BF,1),1);  % XXX: Could be function parameter
    P0 = 10*eye(size(m0,1));       % XXX: Could be function parameter

    % Mean to measurement at step 1
    m0(1) = Y(1);

    %
    % Form the constant matrices
    %
    F  = zeros(size(m0,1));
    H  = zeros(1,size(m0,1));
    L  = zeros(size(m0,1),N+size(BQ,1));
    Qc = zeros(N+size(BQ,1));

    for j=1:N
        i1 = 1+2*(j-1);
        i2 = 2+2*(j-1);
        L(i2,j) = 1;
        Qc(j,j) = qr;
        H(1,i1) = 1;
    end

    %
    % Form the bias model
    %
    i1 = 1+2*N;
    i2 = 2+2*N;
    j1 = N+2-size(BQ,1);
    j2 = N+1;
    F(i1:i2,i1:i2)  = BF;
    L(i1:i2,j1:j2)  = BL;
    H(1,i1:i2)      = BH;
    Qc(j1:j2,j1:j2) = BQ;

%% Run filter

    m = m0;
    P = P0;

    FMM = zeros(size(m,1),length(Y));
    FPP = zeros(size(P,1),size(P,2),length(Y));

    for k=1:length(Y)

        %
        % Update the resonators
        %
        w  = 2*pi*FF(k)/60;
        for j=1:N
            i1 = 1+2*(j-1);
            i2 = 2+2*(j-1);
            F(i1:i2,i1:i2) = [0 j*w; -j*w 0];
        end

        %
        % Discretize
        %
        [A,Q] = lti_disc(F,L,Qc,dt);

        %
        % Estimate with KF
        %
        [m,P] = kf_predict(m,P,A,Q);
        [m,P] = kf_update(m,P,Y(k),H,R);

        FMM(:,k) = m;
        FPP(:,:,k) = P;
    end

%% Run smoother

    % Allocate and set initial
    SMM = zeros(size(FMM));
    SPP = zeros(size(FPP));
    SMM(:,end)   = m;
    SPP(:,:,end) = P;

    for k=length(Y)-1:-1:1

        %
        % Update each resonator
        %
        w  = 2*pi*FF(k+1)/60;
        for j=1:N
            i1 = 1+2*(j-1);
            i2 = 2+2*(j-1);
            F(i1:i2,i1:i2) = [0 j*w; -j*w 0];
        end

        %
        % Discretize
        %
        [A,Q] = lti_disc(F,L,Qc,dt);

        % Smoother
        m_pred = A * FMM(:,k);
        P_pred = A * FPP(:,:,k) * A' + Q;
        D  = FPP(:,:,k) * A' / P_pred;
        m  = FMM(:,k) + D * (m - m_pred);
        P  = FPP(:,:,k) + D * (P - P_pred) * D';

        % Store results
        SMM(:,k)   = m;
        SPP(:,:,k) = P;

    end

    %
    % Extract the signals
    %
    S   = SMM(end-size(BF,1)+1,:);
    CS  = SMM(1:2:(2*nharm),:);
    QCS = SMM(2:2:(2*nharm+1),:);

function [A,Q] = lti_disc(F,L,Q,dt)
% LTI_DISC - Discretize LTI ODE with Gaussian Noise
%
% Syntax:
%   [A,Q] = lti_disc(F,L,Qc,dt)
%
% In:
%   F  - NxN Feedback matrix
%   L  - NxL Noise effect matrix        (optional, default identity)
%   Qc - LxL Diagonal Spectral Density  (optional, default zeros)
%   dt - Time Step                      (optional, default 1)
%
% Out:
%   A - Transition matrix
%   Q - Discrete Process Covariance
%
% Description:
%   Discretize LTI ODE with Gaussian Noise. The original
%   ODE model is in form
%
%     dx/dt = F x + L w,  w ~ N(0,Qc)
%
%   Result of discretization is the model
%
%     x[k] = A x[k-1] + q, q ~ N(0,Q)
%
%   Which can be used for integrating the model
%   exactly over time steps, which are multiples
%   of dt.
%

  %
  % Check number of arguments
  %
  if nargin < 1
    error('Too few arguments');
  end
  if nargin < 2
    L = [];
  end
  if nargin < 3
    Q = [];
  end
  if nargin < 4
    dt = [];
  end

  if isempty(L)
    L = eye(size(F,1));
  end
  if isempty(Q)
    Q = zeros(size(F,1),size(F,1));
  end
  if isempty(dt)
    dt = 1;
  end

  %
  % Closed form integration of transition matrix
  %
  A = expm(F*dt);

  %
  % Closed form integration of covariance
  % by matrix fraction decomposition
  %
  n   = size(F,1);
  Phi = [F L*Q*L'; zeros(n,n) -F'];
  AB  = expm(Phi*dt)*[zeros(n,n);eye(n)];
  Q   = AB(1:n,:)/AB((n+1):(2*n),:);

 function [x,P] = kf_predict(x,P,A,Q,B,u)
% KF_PREDICT - Perform Kalman Filter prediction step
%
% Syntax:
%   [X,P] = KF_PREDICT(X,P,A,Q,B,U)
%
% In:
%   X - Nx1 mean state estimate of previous step
%   P - NxN state covariance of previous step
%   A - Transition matrix of discrete model (optional, default identity)
%   Q - Process noise of discrete model     (optional, default zero)
%   B - Input effect matrix                 (optional, default identity)
%   U - Constant input                      (optional, default empty)
%
% Out:
%   X - Predicted state mean
%   P - Predicted state covariance
%
% Description:
%   Perform Kalman Filter prediction step. The model is
%
%     x[k] = A*x[k-1] + B*u[k-1] + q,  q ~ N(0,Q).
%
%   The predicted state is distributed as follows:
%
%     p(x[k] | x[k-1]) = N(x[k] | A*x[k-1] + B*u[k-1], Q[k-1])
%
%   The predicted mean x-[k] and covariance P-[k] are calculated
%   with the following equations:
%
%     m-[k] = A*x[k-1] + B*u[k-1]
%     P-[k] = A*P[k-1]*A' + Q.
%
%   If there is no input u present then the first equation reduces to
%     m-[k] = A*x[k-1]
%

  %
  % Check arguments
  %
  if nargin < 3
    A = [];
  end
  if nargin < 4
    Q = [];
  end
  if nargin < 5
    B = [];
  end
  if nargin < 6
    u = [];
  end

  %
  % Apply defaults
  %
  if isempty(A)
    A = eye(size(x,1));
  end
  if isempty(Q)
    Q = zeros(size(x,1));
  end
  if isempty(B) & ~isempty(u)
    B = eye(size(x,1),size(u,1));
  end

  %
  % Perform prediction
  %
  if isempty(u)
    x = A * x;
    P = A * P * A' + Q;
  else
    x = A * x + B * u;
    P = A * P * A' + Q;
  end

 function [X,P,K,IM,IS,LH] = kf_update(X,P,y,H,R)
% KF_UPDATE - Kalman Filter update step
%
% Syntax:
%   [X,P,K,IM,IS,LH] = KF_UPDATE(X,P,Y,H,R)
%
% In:
%   X - Nx1 mean state estimate after prediction step
%   P - NxN state covariance after prediction step
%   Y - Dx1 measurement vector.
%   H - Measurement matrix.
%   R - Measurement noise covariance.
%
% Out:
%   X  - Updated state mean
%   P  - Updated state covariance
%   K  - Computed Kalman gain
%   IM - Mean of predictive distribution of Y
%   IS - Covariance or predictive mean of Y
%   LH - Predictive probability (likelihood) of measurement.
%
% Description:
%   Kalman filter measurement update step. Kalman Filter
%   model is
%
%     x[k] = A*x[k-1] + B*u[k-1] + q,  q ~ N(0,Q)
%     y[k] = H*x[k]   + r,             r ~ N(0,R)
%
%   Prediction step of Kalman filter computes predicted
%   mean m-[k] and covariance P-[k] of state:
%
%     p(x[k] | y[1:k-1]) = N(x[k] | m-[k], P-[k])
%
%   See for instance KF_PREDICT how m-[k] and P-[k] are
%   calculated.
%
%   Update step computes the posterior mean m[k] and
%   covariance P[k]  of state given new measurement:
%
%     p(x[k] |�y[1:k]) = N(x[k] |�m[k], P[k])
%
%   Innovation distribution is defined as
%
%     p(y[k] |�y[1:k-1]) = N(y[k] | IM[k], IS[k])
%
%   Updated mean x[k] and covarience P[k] are given by
%   the following equations (not the only possible ones):
%
%     v[k] = y[k] - H[k]*m-[k]
%     S[k] = H[k]*P-[k]*H[k]' + R[k]
%     K[k] = P-[k]*H[k]'*[S[k]]^(-1)
%     m[k] = m-[k] + K[k]*v[k]
%     P[k] = P-[k] - K[k]*S[k]*K[k]'
%
% Example:
%   m = m0;
%   P = P0;
%   M = m0;
%   for i=1:size(Y,2)
%     [m,P] = kf_predict(m,P,A,Q);
%     [m,P] = kf_update(m,P,Y(:,i),H,R);
%     M = [M m];
%   end
%

  %
  % Check which arguments are there
  %
  if nargin < 5
    error('Too few arguments');
  end

  %
  % update step
  %
  IM = H*X;
  IS = (R + H*P*H');
  K = P*H'/IS;
  X = X + K * (y-IM);
  P = P - K*IS*K';
  if nargout > 5
    LH = gauss_pdf(y,IM,IS);
  end

  function [P,E] = gauss_pdf(X,M,S)
%GAUSS_PDF  Multivariate Gaussian PDF
%
% Syntax:
%   [P,E] = GAUSS_PDF(X,M,S)
%
% In:
%   X - Dx1 value or N values as DxN matrix
%   M - Dx1 mean of distibution or N values as DxN matrix.
%   S - DxD covariance matrix
%
% Out:
%   P - Probability of X.
%   E - Negative logarithm of P
%
% Description:
%   Calculate values of PDF (Probability Density
%   Function) of multivariate Gaussian distribution
%
%    N(X | M, S)
%
%   Function returns probability of X in PDF. If multiple
%   X's or M's are given (as multiple columns), function
%   returns probabilities for each of them. X's and M's are
%   repeated to match each other, S must be the same for all.
%

  if size(M,2) == 1
    DX = X-repmat(M,1,size(X,2));
    E = 0.5*sum(DX.*(S\DX),1);
    d = size(M,1);
    E = E + 0.5 * d * log(2*pi) + 0.5 * log(det(S));
    P = exp(-E);
  elseif size(X,2) == 1
    DX = repmat(X,1,size(M,2))-M;
    E = 0.5*sum(DX.*(S\DX),1);
    d = size(M,1);
    E = E + 0.5 * d * log(2*pi) + 0.5 * log(det(S));
    P = exp(-E);
  else
    DX = X-M;
    E = 0.5*DX'*(S\DX);
    d = size(M,1);
    E = E + 0.5 * d * log(2*pi) + 0.5 * log(det(S));
    P = exp(-E);
  end

author: partofthestars

repo: PPGI-Toolbox


text:
function f = note2freq(N)
% MUSIC.NOTE2FREQ Converts a scientific pitch note to a frequency.
%   F = MUSIC.NOTE2FREQ(N) returns the frequency of note N. N is a note string
%   in scientific pitch notation, or a cell array of strings.
%
%   Example
%      f = music.note2freq('F#6')  % returns 1480.0
%
%   See also music.note2tone, music.note2cent, music.freq2note.

%    Author: E. Johnson
%    Copyright 2010 The MathWorks, Inc.

T = music.note2tone(N);
f = music.tone2freq(T);

author: Sable

repo: mcbench-benchmarks


text:
function [net, info] = cnn_mnist_tt(varargin)
% CNN_MNIST_TT  Demonstrated TensorNet on MNIST.

% Fix the random seed.
rng(0);

opts.expDir = fullfile('data','mnist-baseline') ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.dataDir = fullfile('data','mnist') ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.train.batchSize = 100 ;
opts.train.numEpochs = 100 ;
opts.train.continue = true ;
opts.train.gpus = [] ;
opts.train.learningRate = logspace(-2, -5, 45) ;
opts.train.expDir = opts.expDir ;
opts = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = getMnistImdb(opts) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

net = cnn_mnist_init() ;

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------

[net, info] = cnn_train(net, imdb, @getBatch, ...
    opts.train, ...
    'val', find(imdb.images.set == 3)) ;

% --------------------------------------------------------------------
function [im, labels] = getBatch(imdb, batch)
% --------------------------------------------------------------------
im = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;

% --------------------------------------------------------------------
function imdb = getMnistImdb(opts)
% --------------------------------------------------------------------
% Prepare the imdb structure, returns image data with images of size 32 x 32
% and the mean image subtracted.
files = {'train-images-idx3-ubyte', ...
         'train-labels-idx1-ubyte', ...
         't10k-images-idx3-ubyte', ...
         't10k-labels-idx1-ubyte'} ;

if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
end

for i=1:4
  if ~exist(fullfile(opts.dataDir, files{i}), 'file')
    url = sprintf('http://yann.lecun.com/exdb/mnist/%s.gz',files{i}) ;
    fprintf('downloading %s\n', url) ;
    gunzip(url, opts.dataDir) ;
  end
end

f=fopen(fullfile(opts.dataDir, 'train-images-idx3-ubyte'),'r') ;
x1=fread(f,inf,'uint8');
fclose(f) ;
x1=permute(reshape(x1(17:end),28,28,60e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 't10k-images-idx3-ubyte'),'r') ;
x2=fread(f,inf,'uint8');
fclose(f) ;
x2=permute(reshape(x2(17:end),28,28,10e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 'train-labels-idx1-ubyte'),'r') ;
y1=fread(f,inf,'uint8');
fclose(f) ;
y1=double(y1(9:end)')+1 ;

f=fopen(fullfile(opts.dataDir, 't10k-labels-idx1-ubyte'),'r') ;
y2=fread(f,inf,'uint8');
fclose(f) ;
y2=double(y2(9:end)')+1 ;

set = [ones(1,numel(y1)) 3*ones(1,numel(y2))] ;
dataSmall = single(reshape(cat(3, x1, x2),28,28,1,[])) ;
% Fill the image with zeros on the border to resize it to 32 x 32.
data = zeros(32, 32, 1, size(dataSmall, 4)) ;
data(3:30, 3:30, :, :) = dataSmall ;
dataMean = mean(data(:,:,:,set == 1), 4) ;
data = bsxfun(@minus, data, dataMean) ;

imdb.images.data = data ;
imdb.images.data_mean = dataMean ;
imdb.images.labels = cat(2, y1, y2) ;
imdb.images.set = set ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = arrayfun(@(x)sprintf('%d',x),0:9,'uniformoutput',false) ;

author: Bihaqo

repo: TensorNet

save_path: github-repos/MATLAB/Bihaqo-TensorNet


text:
function b = r8gb_vxm ( m, n, ml, mu, a, x )

%*****************************************************************************80
%
%% R8GB_VXM multiplies a vector by a R8GB matrix.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "R8GB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%    LINPACK and LAPACK storage of general band matrices requires
%    an extra ML upper diagonals for possible fill in entries during
%    Gauss elimination.  This routine does not access any entries
%    in the fill in diagonals, because it assumes that the matrix
%    has NOT had Gauss elimination applied to it.  If the matrix
%    has been Gauss eliminated, then the routine R8GB_MU must be
%    used instead.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than min(M,N)-1.
%
%    Input, real A(2*ML+MU+1,N), the R8GB matrix.
%
%    Input, real X(M), the vector to be multiplied by A.
%
%    Output, real B(N), the product X*A.
%
  b(1:n) = 0.0;

  for j = 1 : n
    ilo = max ( 1, j - mu );
    ihi = min ( m, j + ml );
    for i = ilo : ihi
      b(j) = b(j) + x(i) * a(i-j+ml+mu+1,j);
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ev,ee,ebound,xyp] = q1q1grid(x,y,xy,mv,bound,mbound);
%q1q1grid   Q1-Q1 element grid generator
%   [ev,ee,ebound,xyp] = q1q1grid(x,y,xy,mv,bound,mbound);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          xy         nodal coordinate vector
%          mv         Q2 macroelement mapping matrix
%          bound      boundary vertex vector
%          mbound     macroelement boundary vertex vector
%   output
%          ev         element vertex matrix
%          ee         element edge connection matrix
%          ebound     element boundary edge matrix
%          xyp        vertex coordinate vector
%
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
xx=xy(:,1); yy=xy(:,2); nvtx=length(xx);
adj=sparse(nvtx,nvtx);
mel=length(mv(:,1)); nel=4*mel;
ev=zeros(nel,4);
%
%% loop over macroelements
k=1:mel;
% first element
ke=4*k-3;
ev(ke,1)=mv(k,1);
ev(ke,2)=mv(k,5);
ev(ke,3)=mv(k,9);
ev(ke,4)=mv(k,8);
% second element
ke=4*k-2;
ev(ke,1)=mv(k,5);
ev(ke,2)=mv(k,2);
ev(ke,3)=mv(k,6);
ev(ke,4)=mv(k,9);
% third element
ke=4*k-1;
ev(ke,1)=mv(k,9);
ev(ke,2)=mv(k,6);
ev(ke,3)=mv(k,3);
ev(ke,4)=mv(k,7);
% fourth element
ke=4*k;
ev(ke,1)=mv(k,8);
ev(ke,2)=mv(k,9);
ev(ke,3)=mv(k,7);
ev(ke,4)=mv(k,4);
%
%% define element edges
ect=1;
% bottom boundary edges
k1=find(mbound(:,2)==1)';
for k=mbound(k1)
   ebound(ect,1)=4*k-3; ebound(ect+1,1)=4*k-2;
   ebound(ect,2)=1    ; ebound(ect+1,2)=1;
   ect=ect+2;
end
% right boundary edges
k2=find(mbound(:,2)==2)';
for k=mbound(k2)
   ebound(ect,1)=4*k-2; ebound(ect+1,1)=4*k-1;
   ebound(ect,2)=2    ; ebound(ect+1,2)=2;
   ect=ect+2;
end
% top boundary edges
k3=find(mbound(:,2)==3)';
for k=mbound(k3)
   ebound(ect,1)=4*k-1; ebound(ect+1,1)=4*k;
   ebound(ect,2)=3    ; ebound(ect+1,2)=3;
   ect=ect+2;
end
% left boundary edges
k4=find(mbound(:,2)==4)';
for k=mbound(k4)
   ebound(ect,1)=4*k; ebound(ect+1,1)=4*k-3;
   ebound(ect,2)=4    ; ebound(ect+1,2)=4;
   ect=ect+2;
end
%%
xyp=xy;
%
%% compute edge to edge connection array ee
      np=nel;
% initialise global matrices
      adj = sparse(nvtx,nvtx);
      ee = zeros(nel,4);
%
% evaluate element number on each edge in turn
% and assemble into adjacency matrix
%% nx= 0, ny=-1
		 adj=adj + sparse(ev(:,1),ev(:,2),1:np,nvtx,nvtx);
%% nx= 1, ny= 0
		 adj=adj + sparse(ev(:,2),ev(:,3),1:np,nvtx,nvtx);
%% nx= 0, ny= 1
		 adj=adj + sparse(ev(:,3),ev(:,4),1:np,nvtx,nvtx);
%% nx=-1, ny= 0
		 adj=adj + sparse(ev(:,4),ev(:,1),1:np,nvtx,nvtx);
%
       for el=1:nel
		   [ii,jj]=find(adj==el);
           ee(el,:)=diag(adj(jj,ii))';
		   end
           ee=ee(:,[2,4,3,1]);
%
% plotting of the grid
%
%if mel <=64,
	adj=sparse(nvtx,nvtx);
    for i=1:nel
	adj(ev(i,1),ev(i,2)) =1;
	adj(ev(i,2),ev(i,3)) =1;
	adj(ev(i,3),ev(i,4)) =1;
	adj(ev(i,4),ev(i,1)) =1;
    end
    figure(30)
    gplot(adj,xy,'b')
    axis('square')
    hold on
    adj=sparse(nvtx,nvtx);
    k1=find(ebound(:,2)==1);
    for k=1:length(k1)
    kk=ebound(k1(k));
    adj(ev(kk,1),ev(kk,2))=1;
    end
    k2=find(ebound(:,2)==2);
    for k=1:length(k2)
    kk=ebound(k2(k));
    adj(ev(kk,2),ev(kk,3))=1;
    end
    k3=find(ebound(:,2)==3);
    for k=1:length(k3)
    kk=ebound(k3(k));
    adj(ev(kk,3),ev(kk,4))=1;
    end
    k4=find(ebound(:,2)==4);
    for k=1:length(k4)
    kk=ebound(k4(k));
    adj(ev(kk,4),ev(kk,1))=1;
    end
%   gplot(adj,xy,'r')
%   axis('off')
plot(xy(:,1),xy(:,2),'ro')
xybd=xy(bound,:);
plot(xybd(:,1),xybd(:,2),'ko')
plot(xyp(:,1),xyp(:,2),'k*')
hold off
drawnow
title('Q1-Q1 finite element subdivision')
%end
return


text:
%

function [f, df] = sampleEijOpt_6(x)
% optimization code for six wall case

% param
vert = x(6:11);
w_res = x(12);

x0 = x(1); y0 = x(2); % camera center
yc = x(3); % corner
xd = x(4);
yf = x(5);

% assume vertical lines are fixed
a_aob = (vert(2) - vert (1))/w_res*2*pi;
a_boc = (vert(3) - vert (2))/w_res*2*pi;
a_cod = (vert(4) - vert (3))/w_res*2*pi;
a_doe = (vert(5) - vert (4))/w_res*2*pi;
a_eof = (vert(6) - vert (5))/w_res*2*pi;
a_foa = (vert(1) + w_res - vert (6))/w_res*2*pi;

% energy
v_ao  = [x0 y0]; v_bo = [x0-1 y0]; v_co = [x0-1 y0 - yc];
v_do = [x0-xd y0-yc];  v_eo = [x0-xd y0-yf]; v_fo = [x0 y0-yf];

n_v_ao = norm(v_ao); n_v_bo = norm(v_bo); n_v_co = norm(v_co);
n_v_do = norm(v_do); n_v_eo = norm(v_eo); n_v_fo = norm(v_fo);

b_aob = acos(dot(v_ao, v_bo)/n_v_ao/n_v_bo);
if det([v_ao;v_bo]) < 0
    b_aob = 2*pi - b_aob;
end
b_boc = acos(dot(v_bo, v_co)/n_v_bo/n_v_co);
if det([v_bo;v_co]) < 0
    b_boc = 2*pi - b_boc;
end
b_cod = acos(dot(v_co, v_do)/n_v_co/n_v_do);
if det([v_co;v_do]) < 0
    b_cod = 2*pi - b_cod;
end
b_doe = acos(dot(v_do, v_eo)/n_v_do/n_v_eo);
if det([v_do;v_eo]) < 0
    b_doe = 2*pi - b_doe;
end
b_eof = acos(dot(v_eo, v_fo)/n_v_eo/n_v_fo);
if det([v_eo;v_fo]) < 0
    b_eof = 2*pi - b_eof;
end
b_foa = acos(dot(v_fo, v_ao)/n_v_fo/n_v_ao);
if det([v_fo;v_ao]) < 0
    b_foa = 2*pi - b_foa;
end

f = (b_aob - a_aob)^2 + (b_boc - a_boc)^2 + (b_cod - a_cod)^2 + ...
    (b_doe - a_doe)^2 + (b_eof - a_eof)^2 + (b_foa - a_foa)^2;

% gradient
% x0
d_aob_x0 = (2*x0-1)*n_v_ao*n_v_bo + dot(v_ao, v_bo) * (x0*n_v_bo/n_v_ao + (x0-1)*n_v_ao/n_v_bo);
d_aob_x0 = d_aob_x0 * (-1/(sqrt(1-cos(b_aob)*cos(b_aob))+eps))/n_v_ao/n_v_ao/n_v_bo/n_v_bo;
if det([v_ao;v_bo]) < 0
    d_aob_x0 = -d_aob_x0;
end
d_aob_x0 = 2*(b_aob - a_aob) * d_aob_x0;

d_boc_x0 = 2*(x0-1)*n_v_bo*n_v_co + dot(v_bo, v_co) * ((x0-1)*n_v_co/n_v_bo + (x0-1)*n_v_bo/n_v_co);
d_boc_x0 = d_boc_x0 * (-1/(sqrt(1-cos(b_boc)*cos(b_boc))+eps))/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_x0 = -d_boc_x0;
end
d_boc_x0 = 2*(b_boc - a_boc) * d_boc_x0;

d_cod_x0 = (2*x0-xd-1)*n_v_co*n_v_do + dot(v_co, v_do) * ((x0-1)*n_v_do/n_v_co + (x0-xd)*n_v_co/n_v_do);
d_cod_x0 = d_cod_x0 * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_x0 = -d_cod_x0;
end
d_cod_x0 = 2*(b_cod - a_cod) * d_cod_x0;

d_doe_x0 = 2*(x0-xd)*n_v_do*n_v_eo + dot(v_do, v_eo) * ((x0-xd)*n_v_eo/n_v_do + (x0-xd)*n_v_do/n_v_eo);
d_doe_x0 = d_doe_x0 * (-1/(sqrt(1-cos(b_doe)*cos(b_doe))+eps))/n_v_do/n_v_do/n_v_eo/n_v_eo;
if det([v_do;v_eo]) < 0
    d_doe_x0 = -d_doe_x0;
end
d_doe_x0 = 2*(b_doe - a_doe) * d_doe_x0;

d_eof_x0 = (2*x0-xd)*n_v_eo*n_v_fo + dot(v_eo, v_fo) * ((x0-xd)*n_v_fo/n_v_eo + x0*n_v_eo/n_v_fo);
d_eof_x0 = d_eof_x0 * (-1/(sqrt(1-cos(b_eof)*cos(b_eof))+eps))/n_v_eo/n_v_eo/n_v_fo/n_v_fo;
if det([v_eo;v_fo]) < 0
    d_eof_x0 = -d_eof_x0;
end
d_eof_x0 = 2*(b_eof - a_eof) * d_eof_x0;

d_foa_x0 = 2*x0*n_v_fo*n_v_ao + dot(v_fo, v_ao) * (x0*n_v_ao/n_v_fo + x0*n_v_fo/n_v_ao);
d_foa_x0 = d_foa_x0 * (-1/(sqrt(1-cos(b_foa)*cos(b_foa))+eps))/n_v_fo/n_v_fo/n_v_ao/n_v_ao;
if det([v_fo;v_ao]) < 0
    d_foa_x0 = -d_foa_x0;
end
d_foa_x0 = 2*(b_foa - a_foa) * d_foa_x0;

% y0
d_aob_y0 = 2*y0*n_v_ao*n_v_bo + dot(v_ao, v_bo) * (y0*n_v_bo/n_v_ao + y0*n_v_ao/n_v_bo);
d_aob_y0 = d_aob_y0 * (-1/(sqrt(1-cos(b_aob)*cos(b_aob))+eps))/n_v_ao/n_v_ao/n_v_bo/n_v_bo;
if det([v_ao;v_bo]) < 0
    d_aob_y0 = -d_aob_y0;
end
d_aob_y0 = 2*(b_aob - a_aob) * d_aob_y0;

d_boc_y0 = (2*y0-yc)*n_v_bo*n_v_co + dot(v_bo, v_co) * (y0*n_v_co/n_v_bo + (y0-yc)*n_v_bo/n_v_co );
d_boc_y0 = d_boc_y0 * (-1/(sqrt(1-cos(b_boc)*cos(b_boc))+eps))/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_y0 = -d_boc_y0;
end
d_boc_y0 = 2*(b_boc - a_boc) * d_boc_y0;

d_cod_y0 = 2*(y0-yc)*n_v_co*n_v_do + dot(v_co, v_do) * ((y0-yc)/n_v_co*n_v_do + (y0-yc)/n_v_do*n_v_co);
d_cod_y0 = d_cod_y0 * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_y0 = -d_cod_y0;
end
d_cod_y0 = 2*(b_cod - a_cod) * d_cod_y0;

d_doe_y0 = (2*y0-yf-yc)*n_v_do*n_v_eo + dot(v_do, v_eo)*((y0-yc)/n_v_do*n_v_eo + (y0-yf)/n_v_eo*n_v_do);
d_doe_y0 = d_doe_y0 * (-1/(sqrt(1-cos(b_doe)*cos(b_doe))+eps))/n_v_do/n_v_do/n_v_eo/n_v_eo;
if det([v_do;v_eo]) < 0
    d_doe_y0 = -d_doe_y0;
end
d_doe_y0 = 2*(b_doe - a_doe) * d_doe_y0;

d_eof_y0 = 2*(y0-yf)*n_v_eo*n_v_fo + dot(v_eo, v_fo)*((y0-yf)/n_v_eo*n_v_fo + (y0-yf)/n_v_fo*n_v_eo);
d_eof_y0 = d_eof_y0 * (-1/(sqrt(1-cos(b_eof)*cos(b_eof))+eps))/n_v_eo/n_v_eo/n_v_fo/n_v_fo;
if det([v_eo;v_fo]) < 0
    d_eof_y0 = -d_eof_y0;
end
d_eof_y0 = 2*(b_eof - a_eof) * d_eof_y0;

d_foa_y0 = (2*y0-yf)*n_v_fo*n_v_ao + dot(v_fo, v_ao)*((y0-yf)/n_v_fo*n_v_ao + y0/n_v_ao*n_v_fo);
d_foa_y0 = d_foa_y0 * (-1/(sqrt(1-cos(b_foa)*cos(b_foa))+eps))/n_v_fo/n_v_fo/n_v_ao/n_v_ao;
if det([v_fo;v_ao]) < 0
    d_foa_y0 = -d_foa_y0;
end
d_foa_y0 = 2*(b_foa - a_foa) * d_foa_y0;

% yc
d_boc_yc = (-y0)*n_v_bo*n_v_co + dot(v_bo, v_co) * n_v_bo* (yc-y0)/n_v_co;
d_boc_yc = d_boc_yc/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_yc = -d_boc_yc;
end
d_boc_yc = 2*(b_boc - a_boc) * d_boc_yc;

d_cod_yc = 2*(yc-y0)*n_v_co*n_v_do + dot(v_co, v_do)*((yc-y0)/n_v_co*n_v_do + (yc-y0)/n_v_do*n_v_co);
d_cod_yc = d_cod_yc * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_yc = -d_cod_yc;
end
d_cod_yc = 2*(b_cod - a_cod) * d_cod_yc;

d_doe_yc = (yf-y0)*n_v_do*n_v_eo + dot(v_do, v_eo)*(yc-y0)/n_v_do*n_v_eo;
d_doe_yc = d_doe_yc * (-1/(sqrt(1-cos(b_doe)*cos(b_doe))+eps))/n_v_do/n_v_do/n_v_eo/n_v_eo;
if det([v_do;v_eo]) < 0
    d_doe_yc = -d_doe_yc;
end
d_doe_yc = 2*(b_doe - a_doe) * d_doe_yc;

% xd
d_cod_xd = (1-x0)*n_v_co*n_v_do + dot(v_co, v_do)*(xd-x0)*n_v_co/n_v_do;
d_cod_xd = d_cod_xd * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_xd = -d_cod_xd;
end
d_cod_xd = 2*(b_cod - a_cod) * d_cod_xd;

d_doe_xd = 2*(xd-x0)*n_v_do*n_v_eo + dot(v_do, v_eo)*((xd-x0)/n_v_do*n_v_eo + (xd-x0)/n_v_eo*n_v_do);
d_doe_xd = d_doe_xd * (-1/(sqrt(1-cos(b_doe)*cos(b_doe))+eps))/n_v_do/n_v_do/n_v_eo/n_v_eo;
if det([v_do;v_eo]) < 0
    d_doe_xd = -d_doe_xd;
end
d_doe_xd = 2*(b_doe - a_doe) * d_doe_xd;

d_eof_xd = (-x0)*n_v_eo*n_v_fo + dot(v_eo, v_fo)*(xd-x0)/n_v_eo*n_v_fo;
d_eof_xd = d_eof_xd * (-1/(sqrt(1-cos(b_eof)*cos(b_eof))+eps))/n_v_eo/n_v_eo/n_v_fo/n_v_fo;
if det([v_eo;v_fo]) < 0
    d_eof_xd = -d_eof_xd;
end
d_eof_xd = 2*(b_eof - a_eof) * d_eof_xd;

% yf
d_doe_yf = (yc-y0)*n_v_do*n_v_eo + dot(v_do, v_eo)*(yf-y0)/n_v_eo*n_v_do;
d_doe_yf = d_doe_yf * (-1/(sqrt(1-cos(b_doe)*cos(b_doe))+eps))/n_v_do/n_v_do/n_v_eo/n_v_eo;
if det([v_do;v_eo]) < 0
    d_doe_yf = -d_doe_yf;
end
d_doe_yf = 2*(b_doe - a_doe) * d_doe_yf;

d_eof_yf = 2*(yf-y0)*n_v_eo*n_v_fo + dot(v_eo, v_fo)*((yf-y0)/n_v_eo*n_v_fo + (yf-y0)/n_v_fo*n_v_eo);
d_eof_yf = d_eof_yf * (-1/(sqrt(1-cos(b_eof)*cos(b_eof))+eps))/n_v_eo/n_v_eo/n_v_fo/n_v_fo;
if det([v_eo;v_fo]) < 0
    d_eof_yf = -d_eof_yf;
end
d_eof_yf = 2*(b_eof - a_eof) * d_eof_yf;

d_foa_yf = (-y0)*n_v_fo*n_v_ao + dot(v_fo, v_ao)*(yf-y0)/n_v_fo*n_v_ao;
d_foa_yf = d_foa_yf * (-1/(sqrt(1-cos(b_foa)*cos(b_foa))+eps))/n_v_fo/n_v_fo/n_v_ao/n_v_ao;
if det([v_fo;v_ao]) < 0
    d_foa_yf = -d_foa_yf;
end
d_foa_yf = 2*(b_foa - a_foa) * d_foa_yf;

d_x0 = d_aob_x0 + d_boc_x0 + d_cod_x0 + d_doe_x0 + d_eof_x0 + d_foa_x0;
d_y0 = d_aob_y0 + d_boc_y0 + d_cod_y0 + d_doe_y0 + d_eof_y0 + d_foa_y0;
d_yc = d_boc_yc + d_cod_yc + d_doe_yc;
d_xd = d_cod_xd + d_doe_xd + d_eof_xd;
d_yf = d_doe_yf + d_eof_yf + d_foa_yf;

df = [d_x0; d_y0; d_yc; d_xd; d_yf; 0; 0; 0; 0; 0;0;0];

author: zouchuhang

repo: LayoutNet


text:
function [x, mx, sx] = standardise(x, dim, lim)

% STANDARDISE computes the zscore of a matrix along dimension dim
% has similar functionality as the stats-toolbox's zscore function
%
% Use as
%   x = standardise(x, dim)
%
% See also ZSCORE

% Copyright (C) 2009, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin == 1,
  dim = find(size(x)>1,1,'first');
end

if nargin == 3,
  ft_error('third input argument is not used');
end

switch dim
case 1
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(ones(1,n),:,:,:,:,:,:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(ones(1,n),:,:,:,:,:,:,:);
case 2
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,ones(1,n),:,:,:,:,:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,ones(1,n),:,:,:,:,:,:);
case 3
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,ones(1,n),:,:,:,:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,ones(1,n),:,:,:,:,:);
case 4
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,:,ones(1,n),:,:,:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,:,ones(1,n),:,:,:,:);
case 5
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,:,:,ones(1,n),:,:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,:,:,ones(1,n),:,:,:);
case 6
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,:,:,:,ones(1,n),:,:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,:,:,:,ones(1,n),:,:);
case 7
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,:,:,:,:,ones(1,n),:);
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,:,:,:,:,ones(1,n),:);
case 8
  n  = size(x,dim);
  mx = mean(x,dim);
  x  = x-mx(:,:,:,:,:,:,:,ones(1,n));
  sx = sqrt(sum(abs(x).^2,dim)./n);
  x  = x./sx(:,:,:,:,:,:,:,ones(1,n));
otherwise
  ft_error('dim too large, standardise currently supports dimensionality up to 8');
end

author: spm


text:
function a2 = r8vec_copy ( n, a1 )

%*****************************************************************************80
%
%% R8VEC_COPY copies an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the length of the vectors.
%
%    Input, real A1(N), the vector to be copied.
%
%    Output, real A2(N), a copy of A1.
%
  a2(1:n) = a1(1:n);

  return
end


text:
function [pc,r]=circumcenter(p,t)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

nt=size(t,1);
pc=zeros(nt,2);
r=zeros(nt,1);

for it=1:nt
  ct=t(it,:);
  dp1=p(ct(2),:)-p(ct(1),:);
  dp2=p(ct(3),:)-p(ct(1),:);

  mid1=(p(ct(2),:)+p(ct(1),:))/2;
  mid2=(p(ct(3),:)+p(ct(1),:))/2;

  s=[-dp1(2),dp2(2);dp1(1),-dp2(1)]\[-mid1+mid2]';

  cpc=mid1+s(1)*[-dp1(2),dp1(1)];
  cr=norm(p(ct(1),:)-cpc);

  pc(it,:)=cpc;
  r(it,1)=cr;
end

repo: ifem


text:
%fill_mm_bundle Proj. reconstruction from MM [with bundle adjustment].
%
%  Call fill_mm [and bundle adjustment].
%
%  [ P,X, u1,u2, info ] = fill_mm_bundle(M, imsize, nl_params_all_cams [,opt])
%
%  Parameters:
%    M .. measurement matrix (MM) with homogeneous image points with NaNs
%         standing for the unknown elements in all three coordinates
%    imsize .. double(2,m), image sizes: imsize(:,i) is size of image i
%                           m .. No. of cameras
%    opt .. options with default values in ():
%           .no_BA(0) .. whether refine using bundle adjustment
%           .verbose(1) .. whether display info
%           .verbose_short .. see opt in bundle_PX_proj
%           ... other options see in fill_mm
%
%  Return parameters:
%    info.R_lin .. linear estimation of filled M
%    ... other parameters see in fill_mm

function [ P,X, u1,u2, info ] = fill_mm_bundle(M, imsize, nl_params_all_cams, opt)

if nargin < 4, opt = []; end
if ~isfield(opt, 'no_BA')
  opt.no_BA = 0; end
if ~isfield(opt, 'verbose'),
  opt.verbose = 1; end

[P,X, u1,u2, info] = fill_mm(M, opt);

info.R_lin = P*X;

if ~opt.no_BA && length(u1) < size(M,1)/3 && length(u2) < size(M,2)
  if opt.verbose, fprintf(1, 'Bundle adjustment...\n'); tic; end

  [m,n] = size(M); m = m/3; r1 = setdiff(1:m,u1); r2 = setdiff(1:n,u2);
  [P,X] = bundle_PX_proj(P,X, normalize_cut(M(k2i(r1),r2)), imsize, nl_params_all_cams, opt);
  % old bundler:
  %[P,X] = qPXbundle_cmp(P,X, normalize_cut(M(k2i(r1),r2)));

  if opt.verbose, disp(['(' num2str(toc) ' sec)']); end
  info.err.BA = dist(M(k2i(r1),r2), P*X, info.opt.metric);
  if opt.verbose, fprintf('Error (after BA): %f\n', info.err.BA);
  else, fprintf(' %f\n', info.err.BA); end
else, if ~opt.verbose, fprintf('\n'); end; end

author: strawlab

repo: MultiCamSelfCal


text:
function [ v, more ] = triangle_lattice_point_next ( c, v, more )

%*****************************************************************************80
%
%% TRIANGLE_LATTICE_POINT_NEXT returns the next triangle lattice point.
%
%  Discussion:
%
%    The lattice triangle is defined by the vertices:
%
%      (0,0), (C(3)/C(1), 0) and (0,C(3)/C(2))
%
%    The lattice triangle is bounded by the lines
%
%      0 <= X,
%      0 <= Y
%      X / C(1) + Y / C(2) <= C(3)
%
%    Lattice points are listed one at a time, starting at the origin,
%    with X increasing first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer C(3), coefficients defining the
%    lattice triangle.  These should be positive.
%
%    Input/output, integer V(2).  On first call, the input
%    value is not important.  On a repeated call, the input value should
%    be the output value from the previous call.  On output, V contains
%    the next lattice point.
%
%    Input/output, logical MORE.  On input, set MORE to FALSE to indicate
%    that this is the first call for a given triangle.  Thereafter, the input
%    value should be the output value from the previous call.  On output,
%    MORE is TRUE if the returned value V is a new lattice point.
%    If the output value is FALSE, then no more lattice points were found,
%    and V was reset to 0, and the routine should not be called further
%    for this triangle.
%
  n = 2;

  if ( ~more )

    v(1:2) = 0;
    more = 1;

  else

    c1n = i4vec_lcm ( n, c );
    rhs = c1n * c(n+1);

    if ( c(2) * ( v(1) + 1 ) + c(1) * v(2) <= rhs )
      v(1) = v(1) + 1;
    else
      v(1) = 0;
      if ( c(2) * v(1) + c(1) * ( v(2) + 1 ) <= rhs )
        v(2) = v(2) + 1;
      else
        v(2) = 0;
        more = 0;
      end
    end

  end

  return
end


text:
function [varargout] = ft_plot_patch(hdat, vdat, varargin)

% FT_PLOT_PATCH plot a colored shape, similar to the MATLAB patch() function. It is
% similar in usage as ft_plot_vector, and they can be combined, for example,
% to plot an area equivalent to a SEM or STD-DEV around a line.
%
% Use as
%   ft_plot_patch(X, Y, ...)
% where X and Y are similar as the input to the MATLAB patch() function.
%
% Optional arguments should come in key-value pairs and can include
%   'axis'            = draw the local axis,  can be 'yes', 'no', 'xy', 'x' or 'y'
%   'box'             = draw a box around the local axes, can be 'yes' or 'no'
%   'tag'             = string, the name assigned to the object. All tags with the same name can be deleted in a figure, without deleting other parts of the figure.
%   'facecolor'       = see MATLAB standard patch properties
%   'facealpha'       = see MATLAB standard patch properties (note, approx. transparency can be achieved using 'facecolor')
%   'edgecolor'       = see MATLAB standard patch properties (default is 'none') (equivalent to 'linecolor' in PLOT)
%   'linestyle'       = see MATLAB standard patch properties
%   'linewidth'       = see MATLAB standard patch properties
%
% The color of the patchand the edges (i.e. border lines) can be specified in a variety of ways
%   - as a string with one character per line that you want to plot. Supported colors are the same as in PATCH, i.e. 'bgrcmykw'.
%   - as an 'RGB triplet', a 1x3 vector with values between 0 and 1
%   - as 'none' if you do not want the face of the patch to be filled (useful when you want to plot an empty box).
%
% It is possible to plot the object in a local pseudo-axis (c.f. subplot), which is specfied as follows
%   'hpos'            = horizontal position of the center of the local axes
%   'vpos'            = vertical position of the center of the local axes
%   'width'           = width of the local axes
%   'height'          = height of the local axes
%   'hlim'            = horizontal scaling limits within the local axes
%   'vlim'            = vertical scaling limits within the local axes
%
% Example
%   hdat = [1:10 10:-1:1];
%   vdat = rand(1,10);
%   vdat = [vdat vdat(end:-1:1)+1];
%   ft_plot_patch(hdat, vdat)
%
% See also FT_PLOT_VECTOR, PATCH, PLOT

% Copyrights (C) 2015, Roemer van der Meij
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

ws = warning('on', 'MATLAB:divideByZero');

% get the optional input arguments
hpos            = ft_getopt(varargin, 'hpos');
vpos            = ft_getopt(varargin, 'vpos');
width           = ft_getopt(varargin, 'width');
height          = ft_getopt(varargin, 'height');
hlim            = ft_getopt(varargin, 'hlim', 'maxmin');
vlim            = ft_getopt(varargin, 'vlim', 'maxmin');
axis            = ft_getopt(varargin, 'axis', false);
box             = ft_getopt(varargin, 'box', false);
tag             = ft_getopt(varargin, 'tag', '');
parent          = ft_getopt(varargin, 'parent', []);
facecolor       = ft_getopt(varargin, 'facecolor', 'b');
facealpha       = ft_getopt(varargin, 'facealpha', 1);
edgecolor       = ft_getopt(varargin, 'edgecolor', 'none');
linestyle       = ft_getopt(varargin, 'linestyle', 'none');
linewidth       = ft_getopt(varargin, 'linewidth', .5);

% convert the yes/no strings into boolean values
box = istrue(box);

% color management
if ischar(facecolor) && exist([facecolor '.m'], 'file')
  facecolor = eval(facecolor);
end
if ischar(edgecolor) && exist([edgecolor '.m'], 'file')
  edgecolor = eval(edgecolor);
end

% this should be a string, because valid options include yes, no, xy, x, y
if isequal(axis, true)
  axis = 'yes';
elseif isequal(axis, false)
  axis = 'no';
end

% everything is added to the current figure
holdflag = ishold;
if ~holdflag
  hold on
end

if ischar(hlim)
  switch hlim
    case 'maxmin'
      hlim = [min(hdat) max(hdat)];
    case 'maxabs'
      hlim = max(abs(hdat));
      hlim = [-hlim hlim];
    otherwise
      ft_error('unsupported option for hlim')
  end % switch
end % if ischar

if ischar(vlim)
  switch vlim
    case 'maxmin'
      vlim = [min(vdat(:)) max(vdat(:))];
    case 'maxabs'
      vlim = max(abs(vdat(:)));
      vlim = [-vlim vlim];
    otherwise
      ft_error('unsupported option for vlim')
  end % switch
end % if ischar

if vlim(1)==vlim(2)
  % vertical scaling cannot be determined, behave consistent to the plot() function
  vlim = [-1 1];
end

% these must be floating point values and not integers, otherwise the scaling fails
hdat = double(hdat);
vdat = double(vdat);
hlim = double(hlim);
vlim = double(vlim);

if isempty(hpos) && ~isempty(hlim)
  hpos = (hlim(1)+hlim(2))/2;
end
if isempty(vpos) && ~isempty(vlim)
  vpos = (vlim(1)+vlim(2))/2;
end

if isempty(width) && ~isempty(hlim)
  width = hlim(2)-hlim(1);
end

if isempty(height) && ~isempty(vlim)
  height = vlim(2)-vlim(1);
end

% first shift the horizontal axis to zero
if any(hlim) ~= 0
  hdat = hdat - (hlim(1)+hlim(2))/2;
  % then scale to length 1
  if hlim(2)-hlim(1)~=0
    hdat = hdat ./ (hlim(2)-hlim(1));
  else
    hdat = hdat /hlim(1);
  end
  % then scale to the new width
  hdat = hdat .* width;
end
% then shift to the new horizontal position
hdat = hdat + hpos;

if any(vlim) ~= 0
  % first shift the vertical axis to zero
  vdat = vdat - (vlim(1)+vlim(2))/2;
  % then scale to length 1
  vdat = vdat / (vlim(2)-vlim(1));
  % then scale to the new width
  vdat = vdat .* height;
end
% then shift to the new vertical position
vdat = vdat + vpos;

% plot the patch
h = patch(hdat, vdat, facecolor, 'facealpha', facealpha, 'edgecolor', edgecolor, 'linestyle', linestyle, 'linewidth', linewidth);

if box
  % this plots a box around the original hpos/vpos with appropriate width/height
  x1 = hpos - width/2;
  x2 = hpos + width/2;
  y1 = vpos - height/2;
  y2 = vpos + height/2;

  X = [x1 x2 x2 x1 x1];
  Y = [y1 y1 y2 y2 y1];
  h = line(X, Y);
  set(h, 'color', 'k');

  % this plots a box around the original hpos/vpos with appropriate width/height
  % boxposition = zeros(1,4);
  % boxposition(1) = hpos - width/2;
  % boxposition(2) = hpos + width/2;
  % boxposition(3) = vpos - height/2;
  % boxposition(4) = vpos + height/2;
  % ft_plot_box(boxposition, 'facecolor', 'none', 'edgecolor', 'k');

  % this plots a box around the complete data
  % boxposition = zeros(1,4);
  % boxposition(1) = hlim(1);
  % boxposition(2) = hlim(2);
  % boxposition(3) = vlim(1);
  % boxposition(4) = vlim(2);
  % ft_plot_box(boxposition, 'hpos', hpos, 'vpos', vpos, 'width', width, 'height', height, 'hlim', hlim, 'vlim', vlim);
end

if ~isempty(axis) && ~strcmp(axis, 'no')
  switch axis
    case {'yes' 'xy'}
      xaxis = true;
      yaxis = true;
    case {'x'}
      xaxis = true;
      yaxis = false;
    case {'y'}
      xaxis = false;
      yaxis = true;
    otherwise
      ft_error('invalid specification of the "axis" option')
  end

  if xaxis
    % x-axis should touch 0,0
    xrange = hlim;
    if sign(xrange(1))==sign(xrange(2))
      [dum minind] = min(abs(hlim));
      xrange(minind) = 0;
    end
    ft_plot_line(xrange, [0 0], 'hpos', hpos, 'vpos', vpos, 'hlim', hlim, 'vlim', vlim, 'width', width, 'height', height);
  end
  if yaxis
    % y-axis should touch 0,0
    yrange = vlim;
    if sign(yrange(1))==sign(yrange(2))
      [dum minind] = min(abs(vlim));
      yrange(minind) = 0;
    end
    ft_plot_line([0 0], yrange, 'hpos', hpos, 'vpos', vpos, 'hlim', hlim, 'vlim', vlim, 'width', width, 'height', height);
  end
end

set(h, 'tag', tag);

if ~isempty(parent)
  set(h, 'Parent', parent);
end

% the (optional) output is the handle
if nargout == 1;
  varargout{1} = h;
end

if ~holdflag
  hold off
end

warning(ws); % revert to original state


text:
% Color Similarity Non-local Pixel Affinities
% This function implements the affinity based on color differences
% first used for image matting in the paper
% Qifeng Chen, Dingzeyu Li, Chi-Keung Tang, "KNN Matting", IEEE
% TPAMI, 2013.
% All parameters other than image are optional. The output is a sparse
% matrix which has non-zero element for the non-local neighbors of
% the pixels given by binary map inMap.
% - K defines the number of neighbors from which LLE weights are
%   computed.
% - outMap is a binary map that defines where the nearest neighbor
%   search is done.
% - xyWeight determines how much importance is given to the spatial
%   coordinates in the nearest neighbor selection.
% - When useHSV is false (default), the search is done i [r g b x y] space,
%   otherwise the feature space is [cos(h) sin(h), s, v, x, y].

function Wcs = colorSimilarityAffinities(image, K, inMap, outMap, xyWeight, useHSV)

    [h, w, ~] = size(image);
    N = h * w;

    if ~exist('K', 'var') || isempty(K)
        K = 5;
    end
    if ~exist('inMap', 'var') || isempty(inMap)
        inMap = true(h, w);
    end
    if ~exist('outMap', 'var') || isempty(outMap)
        outMap = true(h, w);
    end
    if ~exist('xyWeight', 'var') || isempty(xyWeight)
        xyWeight = 0.05;
    end
    if ~exist('useHSV', 'var') || isempty(useHSV)
        useHSV = false;
    end

    if useHSV
        image = rgb2hsv(image);
        image = cat(3, cos(image(:, :, 1)) * 2 * pi, sin(image(:, :, 1)) * 2 * pi, image(:,:,2:3));
    end

    [~, neighInd, ~] = findNonlocalNeighbors(image, K, xyWeight, inMap, outMap);

    % This behaviour below, decreasing the xy-weight and finding a new set of neighbors, is taken
    % from the public implementation of KNN matting by Chen et al.
    [inInd, neighInd2, features] = findNonlocalNeighbors(image, ceil(K / 5), xyWeight / 100, inMap, outMap);
    neighInd = [neighInd, neighInd2];
    features(:, end-1 : end) = features(:, end-1 : end) / 100;

    inInd = repmat(inInd, [1, size(neighInd, 2)]);
    flows = max(1 - sum(abs(features(inInd(:), :) - features(neighInd(:), :)), 2) / size(features, 2), 0);

    Wcs = sparse(inInd(:), neighInd(:), flows, N, N);
    Wcs = (Wcs + Wcs') / 2; % If p is a neighbor of q, make q a neighbor of p
end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
function ASK_FSK_PSK(msglen)
%msglen= number of bits to be transmitted
%take msglen=10000, or 20000 for more accuracy
%If you have any problem or feedback please contact me @
%%===============================================
% NIKESH BAJAJ
% Asst. Prof., Lovely Professional University, India
% Almameter: Aligarh Muslim University, India
% +919915522564, bajaj.nikkey@gmail.com
%%===============================================
n=msglen;
b=randint(1,n);
f1=1;f2=2;
t=0:1/30:1-1/30;
%ASK
sa1=sin(2*pi*f1*t);
E1=sum(sa1.^2);
sa1=sa1/sqrt(E1); %unit energy
sa0=0*sin(2*pi*f1*t);
%FSK
sf0=sin(2*pi*f1*t);
E=sum(sf0.^2);
sf0=sf0/sqrt(E);
sf1=sin(2*pi*f2*t);
E=sum(sf1.^2);
sf1=sf1/sqrt(E);
%PSK
sp0=-sin(2*pi*f1*t)/sqrt(E1);
sp1=sin(2*pi*f1*t)/sqrt(E1);

%MODULATION
ask=[];psk=[];fsk=[];
for i=1:n
    if b(i)==1
        ask=[ask sa1];
        psk=[psk sp1];
        fsk=[fsk sf1];
    else
        ask=[ask sa0];
        psk=[psk sp0];
        fsk=[fsk sf0];
    end
end
figure(1)
subplot(411)
stairs(0:10,[b(1:10) b(10)],'linewidth',1.5)
axis([0 10 -0.5 1.5])
title('Message Bits');grid on
subplot(412)
tb=0:1/30:10-1/30;
plot(tb, ask(1:10*30),'b','linewidth',1.5)
title('ASK Modulation');grid on
subplot(413)
plot(tb, fsk(1:10*30),'r','linewidth',1.5)
title('FSK Modulation');grid on
subplot(414)
plot(tb, psk(1:10*30),'k','linewidth',1.5)
title('PSK Modulation');grid on
xlabel('Time');ylabel('Amplitude')
%AWGN
for snr=0:20
    askn=awgn(ask,snr);
    pskn=awgn(psk,snr);
    fskn=awgn(fsk,snr);

    %DETECTION
    A=[];F=[];P=[];
    for i=1:n
        %ASK Detection
        if sum(sa1.*askn(1+30*(i-1):30*i))>0.5
            A=[A 1];
        else
            A=[A 0];
        end
        %FSK Detection
        if sum(sf1.*fskn(1+30*(i-1):30*i))>0.5
            F=[F 1];
        else
            F=[F 0];
        end
        %PSK Detection
        if sum(sp1.*pskn(1+30*(i-1):30*i))>0
            P=[P 1];
        else
            P=[P 0];
        end
    end

    %BER
    errA=0;errF=0; errP=0;
    for i=1:n
        if A(i)==b(i)
            errA=errA;
        else
            errA=errA+1;
        end
        if F(i)==b(i)
            errF=errF;
        else
            errF=errF+1;
        end
        if P(i)==b(i)
            errP=errP;
        else
            errP=errP+1;
        end
    end
    BER_A(snr+1)=errA/n;
    BER_F(snr+1)=errF/n;
    BER_P(snr+1)=errP/n;
end

figure(2)
subplot(411)
stairs(0:10,[b(1:10) b(10)],'linewidth',1.5)
axis([0 10 -0.5 1.5]);grid on
title('Received signal after AWGN Channel')
subplot(412)
tb=0:1/30:10-1/30;
plot(tb, askn(1:10*30),'b','linewidth',1.5)
title('Received ASK signal');grid on
subplot(413)
plot(tb, fskn(1:10*30),'r','linewidth',1.5)
title('Received FSK signal');grid on
subplot(414)
plot(tb, pskn(1:10*30),'k','linewidth',1.5)
title('Received PSK signal');grid on
figure(3)
semilogy(0:20,BER_A, 'b','linewidth',2)
title('BER Vs SNR')
grid on;
hold on
semilogy(0:20,BER_F,'r','linewidth',2)
semilogy(0:20,BER_P, 'k','linewidth',2)
xlabel('Eo/No(dB)')
ylabel('BER')
hold off
legend('ASK','FSK','PSK');


text:
classdef MimBrushColourEditBox < MimGuiPluginEditBox
    % MimBrushColourEditBox. Gui Plugin for changing the size of the edit brush
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Custom'
        SelectedText = 'Custom'

        ToolTip = 'Change the colour of the editing paint brush'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Location = 40

        MinValue = 0
        MaxValue = 255
        DefaultValue = 1

        EditBoxPosition = 90
        EditBoxWidth = 25

        StackVertically = false
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = true;
        end

        function [value_instance_handle, value_property_name, limits_instance_handle, limits_property_name] = GetHandleAndProperty(gui_app)
            value_instance_handle = gui_app.ImagePanel;
            value_property_name = 'PaintBrushColour';
            limits_instance_handle = [];
            limits_property_name = [];
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function res = utl_aggregate_results(varargin)
% Internal. Aggregate the given results (in any format allowed for ml_predict) into a single array.
% In:
%   Results... : results as produced by ml_predict
%
% Out:
%   Aggregate : all results concatenated into a single array
%
% See also:
%   ml_predict
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-07

if isnumeric(varargin{1})
    try
        res = vertcat(varargin{:});
    catch e
        sizes = cellfun(@size,varargin,'UniformOutput',false);
        error('Failed to concatenate numeric prediction arrays -- make sure that the dimensionality of your predictions (%s) matches in each array; error message: %s.',hlp_tostring(sizes),e.message);
    end
elseif iscell(varargin{1}) && length(varargin{1}) > 1
    % make sure that the result formats are all equal
    v1 = varargin{1}([1 3:end]);
    for i=2:length(varargin)
        if ~isequal(v1,varargin{i}([1 3:end]))
            error('The formats of all prediction arrays to concatenate need to be identical.'); end
    end
    % grab the actual parameters
    params = cellfun(@(x)x{2},varargin,'UniformOutput',false);
    if isnumeric(varargin{1}{2})
        % numeric parameters are vertically concatenated
        try
            res = {varargin{1}{1} vertcat(params{:}) varargin{1}{3:end}};
        catch e
            sizes = cellfun(@size,params,'UniformOutput',false);
            error('Failed to concatenate prediction arrays -- make sure that the dimensionality of your distributions (%s) matches in each array; error message: %s.',hlp_tostring(sizes),e.message);
        end
    elseif iscell(varargin{1}{2})
        % cell parameters are horizontally concatenated
        try
            res = {varargin{1}{1} horzcat(params{:}) varargin{1}{3:end}};
        catch e
            sizes = cellfun(@size,params,'UniformOutput',false);
            error('Failed to concatenate prediction arrays -- make sure that the dimensionality of your distributions (%s) matches in each array; error message: %s.',hlp_tostring(sizes),e.message);
        end
    else
        error('This function can only concatenate prediction arrays of the form {X,Y, ...} where Y is either numeric or a cell array.');
    end
else
    error('This function can only concatenate arrays of predictions that are either all numeric or in the cell array form described in ml_predict.');
end

repo: REST


text:
%% DEMO 16: TIGRE 2D tomography
%
%
%  In demo we show how to do 2D tomography in TIGRE. It is exactly the same
%  as 3D.
%
%
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
%% Initialize
clear;
close all;
%% PARALLEL BEAM 2D

% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
geo.DSD = 1536;                             % Distance Source Detector      (mm)
geo.DSO = 1000;                             % Distance Source Origin        (mm)
% Detector parameters

%
% 1- Make the detector Nx1
% 2- Make the image NxMx1

% Image parameters
geo.nVoxel=[256;256;1];                   % number of voxels              (vx)
geo.sVoxel=[256;256;1];                   % total size of the image       (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)

geo.nDetector=[512;  1];					% number of pixels              (px)
geo.dDetector=[0.8; geo.dVoxel(3)]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)

% MAKE SURE THAT THE DETECTOR PIXELS SIZE IN V IS THE SAME AS THE IMAGE!

% Offsets
geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)
geo.offDetector=[0; 0];

% Auxiliary
geo.accuracy=0.5;                           % Accuracy of FWD proj          (vx/sample)
geo.mode='parallel';
%% Define angles of projection and load phatom image

angles=linspace(0,2*pi,100);
phatom=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
projections=Ax(phatom,geo,angles);
%% recosntruct

imgOSSART=OS_SART(projections,geo,angles,40);
imgASDPOCS=ASD_POCS(projections,geo,angles,40);

%% Plot
plotImg([imgOSSART;imgASDPOCS ],'Dim',3,'Slice',1); %top is ASD_POCS, bottom OSSART

%%
%%
%% And now Fan Beam
%%
% The same thing!

%% FAN BEAM 2D

% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
geo.DSD = 1536;                             % Distance Source Detector      (mm)
geo.DSO = 1000;                             % Distance Source Origin        (mm)
% Detector parameters

% Image parameters
geo.nVoxel=[256;256;1];                   % number of voxels              (vx)
geo.sVoxel=[256;256;1];                   % total size of the image       (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)

geo.nDetector=[512;  1];					% number of pixels              (px)
geo.dDetector=[0.8; 1]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)

% Offsets
geo.offOrigin =[0;0;0];                        % Offset of image from origin   (mm)              ;
geo.offDetector=[0; 0];

% Auxiliary
geo.accuracy=0.5;                           % Accuracy of FWD proj          (vx/sample)
geo.mode='cone';
%% Define angles of projection and load phatom image

angles=linspace(0,2*pi,100);
phatom=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
% phatom=cat(3,phatom,phatom);
projections=Ax(phatom,geo,angles,'interpolated');
%% recosntruct

imgOSSART=OS_SART(projections,geo,angles,40);
% in 2D, the TV norm minimization happens in CPU, this causes the
% algorithms to be slower
imgASDPOCS=ASD_POCS(projections,geo,angles,40);

%% Plot
plotImg([imgOSSART;imgASDPOCS ],'Dim',3,'Slice',1); %top is ASD_POCS, nottom OSSART

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (var_elim)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)

% we could pre-process the evidence here, to prevent repeated work, but we don't.
engine.evidence = evidence;

if nargout == 2
  [m, loglik] = marginal_nodes(engine, [1]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% RANGE = showSpyr (PYR, INDICES, RANGE, GAP, LEVEL_SCALE_FACTOR)
%
% Display a steerable pyramid, specified by PYR and INDICES
% (see buildSpyr), in the current figure.  The highpass band is not shown.
%
% RANGE is a 2-vector specifying the values that map to black and
% white, respectively.  These values are scaled by
% LEVEL_SCALE_FACTOR^(lev-1) for bands at each level.  Passing a value
% of 'auto1' sets RANGE to the min and max values of MATRIX.  'auto2'
% sets RANGE to 3 standard deviations below and above 0.0.  In both of
% these cases, the lowpass band is independently scaled.  A value of
% 'indep1' sets the range of each subband independently, as in a call
% to showIm(subband,'auto1').  Similarly, 'indep2' causes each subband
% to be scaled independently as if by showIm(subband,'indep2').
% The default value for RANGE is 'auto2'.
%
% GAP (optional, default=1) specifies the gap in pixels to leave
% between subbands.
%
% LEVEL_SCALE_FACTOR indicates the relative scaling between pyramid
% levels.  This should be set to the sum of the kernel taps of the
% lowpass filter used to construct the pyramid (default is 2, which is
% correct for L2-normalized filters.

% Eero Simoncelli, 2/97.

function [range] = showSpyr(pyr, pind, range, gap, scale);

nbands = spyrNumBands(pind);

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('range') ~= 1)
  range = 'auto2';
end

if (exist('gap') ~= 1)
  gap = 1;
end

if (exist('scale') ~= 1)
  scale = 2;
end

%------------------------------------------------------------

ht = spyrHt(pind);
nind = size(pind,1);

%% Auto range calculations:
if strcmp(range,'auto1')
  range = ones(nind,1);
  band = spyrHigh(pyr,pind);
  [mn,mx] = range2(band);
  for lnum = 1:ht
    for bnum = 1:nbands
      band = spyrBand(pyr,pind,lnum,bnum)/(scale^(lnum-1));
      range((lnum-1)*nbands+bnum+1) = scale^(lnum-1);
      [bmn,bmx] = range2(band);
      mn = min(mn, bmn);
      mx = max(mx, bmx);
    end
  end
  range = range * [mn mx]; 		% outer product
  band = pyrLow(pyr,pind);
  [mn,mx] = range2(band);
  range(nind,:) = [mn, mx];

elseif strcmp(range,'indep1')
  range = zeros(nind,2);
  for bnum = 1:nind
    band = pyrBand(pyr,pind,bnum);
    [mn,mx] = range2(band);
    range(bnum,:) =  [mn mx];
  end

elseif strcmp(range,'auto2')
  range = ones(nind,1);
  band = spyrHigh(pyr,pind);
  sqsum = sum(sum(band.^2));  numpixels = prod(size(band));
  for lnum = 1:ht
    for bnum = 1:nbands
      band = spyrBand(pyr,pind,lnum,bnum)/(scale^(lnum-1));
      sqsum = sqsum + sum(sum(band.^2));
      numpixels = numpixels + prod(size(band));
      range((lnum-1)*nbands+bnum+1) = scale^(lnum-1);
    end
  end
  stdev = sqrt(sqsum/(numpixels-1));
  range = range * [ -3*stdev 3*stdev ]; % outer product
  band = pyrLow(pyr,pind);
  av = mean2(band);   stdev = sqrt(var2(band));
  range(nind,:) = [av-2*stdev,av+2*stdev];

elseif strcmp(range,'indep2')
  range = zeros(nind,2);
  for bnum = 1:(nind-1)
    band = pyrBand(pyr,pind,bnum);
    stdev = sqrt(var2(band));
    range(bnum,:) =  [ -3*stdev 3*stdev ];
  end
  band = pyrLow(pyr,pind);
  av = mean2(band);   stdev = sqrt(var2(band));
  range(nind,:) = [av-2*stdev,av+2*stdev];

elseif isstr(range)
  error(sprintf('Bad RANGE argument: %s',range))

elseif ((size(range,1) == 1) & (size(range,2) == 2))
  scales = scale.^[0:(ht-1)];
  scales = ones(nbands,1) * scales;   %outer product
  scales = [1; scales(:); scale^ht];  %tack on highpass and lowpass
  range = scales * range;		% outer product
  band = pyrLow(pyr,pind);
  range(nind,:) = range(nind,:) + mean2(band) - mean(range(nind,:));

end

% CLEAR FIGURE:
clf;

colormap(gray);
cmap = get(gcf,'Colormap');
nshades = size(cmap,1);

%  Find background color index:
clr = get(gcf,'Color');
bg = 1;
dist = norm(cmap(bg,:)-clr);
for n = 1:nshades
  ndist = norm(cmap(n,:)-clr);
  if (ndist < dist)
    dist = ndist;
    bg = n;
  end
end

%% Compute positions of subbands:
llpos = ones(nind,2);

if (nbands == 2)
  ncols = 1;  nrows = 2;
else
  ncols = ceil((nbands+1)/2);   nrows = ceil(nbands/2);
end
relpos = [ (1-nrows):0, zeros(1,(ncols-1)); ...
           zeros(1,nrows), -1:-1:(1-ncols) ]';
if (nbands > 1)
  mvpos = [-1 -1];
else
  mvpos = [0 -1];
end
basepos = [0 0];

for lnum = 1:ht
  ind1 = (lnum-1)*nbands + 2;
  sz = pind(ind1,:)+gap;
  basepos = basepos + mvpos .* sz;
  if (nbands < 5)			% to align edges...
    sz = sz + gap*(ht-lnum+1);
  end
  llpos(ind1:ind1+nbands-1,:) = relpos * diag(sz) + ones(nbands,1)*basepos;
end

% lowpass band
sz = pind(nind-1,:)+gap;
basepos = basepos + mvpos .* sz;
llpos(nind,:) = basepos;

%% Make position list positive, and allocate appropriate image:
llpos = llpos - ones(nind,1)*min(llpos) + 1;
llpos(1,:) = [1 1];
urpos = llpos + pind - 1;
d_im = bg + zeros(max(urpos));

%% Paste bands into image, (im-r1)*(nshades-1)/(r2-r1) + 1.5
for bnum=2:nind
  mult = (nshades-1) / (range(bnum,2)-range(bnum,1));
  d_im(llpos(bnum,1):urpos(bnum,1), llpos(bnum,2):urpos(bnum,2)) = ...
      mult*pyrBand(pyr,pind,bnum) + (1.5-mult*range(bnum,1));
end

hh = image(d_im);
axis('off');
pixelAxes(size(d_im),'full');
set(hh,'UserData',range);

repo: SelfExSR


text:
% funSimLogNormProbCovFade returns SINR-based k-coverage probability
% under Rayleigh (mean one) fading model and log-normal shadowing based on
% repeated simulations of of model outlined in [1]
%
% simPCovFade=funSimLogNormProbCovFade(tValues,betaConst,K,lambda,sigma,W,diskRadius,simNumb)
% simPCovFade is the 1-coverage probability
% tValues are the SINR threshold values. tValues can be a vector
% betaConst is the pathloss exponent
% K = path-loss constant
% lambda = density of base station/nodes of cellular network
% sigma = standard deviation of log-normal shadowing (not in dB)
% W = noise constant
% diskRadius = radius of simulation region (Warning: if too small, edge
% effects will not be sufficiently included and disagreement with analytic
% results may occurr)
% simNumb = number of simulations
% All input values (except tValues) are scalars
%
% Author: H.P. Keeler, Inria Paris/ENS, 2013
%
% References
% [1] H.P. Keeler, B. Błaszczyszyn and M. Karray,
% 'SINR-based k-coverage probability in cellular networks with arbitrary
% shadowing', accepted at ISIT, 2013

function simPCovFade=funSimLogNormProbCovFade(tValues,betaConst,K,lambda,sigma,W,diskRadius,simNumb)

tNumb=length(tValues);

%%% Simulation Section %%%
%(uniformly) randomly places nodes on a disk of radius diskRadius
diskArea=pi*diskRadius^2;
coveredNumb=zeros(size(tValues));
ESTwoBeta=exp(sigma^2*(2-betaConst)/betaConst^2);
%rescale lambda - see foot note 5 in [1]
lambdaSim=lambda*ESTwoBeta;
for i=1:simNumb

    randNumb=poissrnd(lambdaSim*diskArea);
    %shadowing distribution can be constant if lambda is rescaled - see [1], [2] and [3]
    shadowRand=ones(randNumb,1);
    fadeRand=exprnd(1,randNumb,1); %Rayleigh fading corresponds to exponential random variables
    %random distances from the typical node
    rRand=diskRadius*sqrt(rand(randNumb,1)); %uniform in cartesion, not polar coordinates
    signalRand=(K*rRand).^(betaConst)./shadowRand;
    [Y1 indexY1]=min(signalRand);    %find Y_1 (first order statistics)

    interferTotal=sum(fadeRand./signalRand); %total inteference in network
    newExp=fadeRand(indexY1); %use corresponding exponential variable
    SINR=(newExp/Y1)./((interferTotal-(newExp/Y1))+W);
    for j=1:tNumb
        T=tValues(j);
        %counts how many nodes are connected/covered
        if sum(SINR>=T)>=1
            coveredNumb(j)=coveredNumb(j)+1;
        end
    end
end

simPCovFade=coveredNumb/simNumb;


text:
function Rb2s = flu2rdf

% FLU2RDF Camera body to camera sensor rotation matrix.
%   FLU2RDF computes the rotation matrix for a camera whos body
%   is in the FLU frame (x-front, y-left, z-up) and its sensor in
%   the RDF frame (x-right, y-down, z-front).

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

Rb2s = e2R([-pi/2 0 -pi/2]');

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb