text:
function [estimate] = ft_inverse_dics(sourcemodel, sens, headmodel, dat, C, varargin)

% FT_INVERSE_DICS scans on pre-defined dipole locations with a single dipole
% and returns the beamformer spatial filter output for a dipole on every
% location.
%
% Use as
%   [estimate] = ft_inverse_dics(sourcemodel, sens, headmodel, dat, cov, ...)
% where
%   sourcemodel is the input source model, see FT_PREPARE_SOURCEMODEL
%   sens        is the gradiometer or electrode definition, see FT_DATATYPE_SENS
%   headmodel   is the volume conductor definition, see FT_PREPARE_HEADMODEL
%   dat         is the data matrix with the ERP or ERF
%   cov         is the data covariance or cross-spectral density matrix
% and
%   estimate    contains the estimated source parameters
%
% Additional input arguments should be specified as key-value pairs and can include
%   'Pr'               = power of the external reference channel
%   'Cr'               = cross spectral density between all data channels and the external reference channel
%   'refdip'           = location of dipole with which coherence is computed
%   'powmethod'        = can be 'trace' or 'lambda1'
%   'feedback'         = can be 'none', 'gui', 'dial', 'textbar', 'text', 'textcr', 'textnl' (default = 'text')
%   'fixedori'         = use fixed or free orientation,                 can be 'yes' or 'no'
%   'projectnoise'     = project noise estimate through filter,         can be 'yes' or 'no'
%   'realfilter'       = construct a real-valued filter,                can be 'yes' or 'no'
%   'keepfilter'       = remember the beamformer filter,                can be 'yes' or 'no'
%   'keepleadfield'    = remember the forward computation,              can be 'yes' or 'no'
%   'keepcsd'          = remember the estimated cross-spectral density, can be 'yes' or 'no'
%   'weightnorm'       = normalize the beamformer weights,              can be 'no', 'unitnoisegain' or 'nai'
%
% These options influence the forward computation of the leadfield
%   'reducerank'      = 'no' or number (default = 3 for EEG, 2 for MEG)
%   'backproject'     = 'yes' or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace (default = 'yes')
%   'normalize'       = 'no', 'yes' or 'column' (default = 'no')
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or Nx1 vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%
% These options influence the mathematical inversion of the cross-spectral density matrix
%  'lambda'           = regularisation parameter
%  'kappa'            = parameter for covariance matrix inversion
%  'tol'              = parameter for covariance matrix inversion
%
% If the dipole definition only specifies the dipole location, a rotating
% dipole (regional source) is assumed on each location. If a dipole moment
% is specified, its orientation will be used and only the strength will
% be fitted to the data.
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Copyright (C) 2003-2020, Robert Oostenveld and 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 mod(nargin-5,2)
  % the first 5 arguments are fixed, the other arguments should come in pairs
  ft_error('invalid number of optional arguments');
end

% get the optional input arguments, or use defaults
Pr             = ft_getopt(varargin, 'Pr');
Cr             = ft_getopt(varargin, 'Cr');
refdip         = ft_getopt(varargin, 'refdip');
powmethod      = ft_getopt(varargin, 'powmethod');  % the default for this is set below
realfilter     = ft_getopt(varargin, 'realfilter'); % the default for this is set below
subspace       = ft_getopt(varargin, 'subspace');
feedback       = ft_getopt(varargin, 'feedback',      'text');
keepcsd        = ft_getopt(varargin, 'keepcsd',       'no');
keepfilter     = ft_getopt(varargin, 'keepfilter',    'no');
keepleadfield  = ft_getopt(varargin, 'keepleadfield', 'no');
projectnoise   = ft_getopt(varargin, 'projectnoise',  'yes');
fixedori       = ft_getopt(varargin, 'fixedori',      'no');
weightnorm     = ft_getopt(varargin, 'weightnorm',    'no');

% construct the low-level options for the covariance matrix inversion as key-value pairs, these are passed to FT_INV
invopt = {};
invopt = ft_setopt(invopt, 'lambda',    ft_getopt(varargin, 'lambda', 0));
invopt = ft_setopt(invopt, 'kappa',     ft_getopt(varargin, 'kappa'));
invopt = ft_setopt(invopt, 'tolerance', ft_getopt(varargin, 'tol'));
invopt = ft_setopt(invopt, 'method',    ft_getopt(varargin, 'invmethod'));

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(varargin, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(varargin, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(varargin, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(varargin, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(varargin, 'weight'));

% convert the yes/no arguments to the corresponding logical values
keepfilter     = istrue(keepfilter);
keepleadfield  = istrue(keepleadfield);
keepcsd        = istrue(keepcsd);
projectnoise   = istrue(projectnoise);
fixedori       = istrue(fixedori);

% FIXME besides regular/complex lambda1, also implement a real version

% default is to use the largest singular value of the csd matrix, see Gross 2001
if isempty(powmethod)
  powmethod = 'lambda1';
end

% default is to be consistent with the original description of DICS in Gross 2001
if isempty(realfilter)
  realfilter = 'no';
end

% use these two logical flags instead of doing the string comparisons each time again
powtrace   = strcmp(powmethod, 'trace');
powlambda1 = strcmp(powmethod, 'lambda1');

if ~isempty(Cr)
  % ensure that the cross-spectral density with the reference signal is a column matrix
  Cr = Cr(:);
end

if isfield(sourcemodel, 'mom') && fixedori
  ft_error('you cannot specify a dipole orientation and fixedori simultaneously');
end

% flags to avoid calling isfield repeatedly in the loop over grid positions (saves a lot of time)
hasmom        = isfield(sourcemodel, 'mom');
hasleadfield  = isfield(sourcemodel, 'leadfield');
hasfilter     = isfield(sourcemodel, 'filter');
hassubspace   = isfield(sourcemodel, 'subspace');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the dipole positions that are inside/outside the brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(sourcemodel, 'inside')
  if hasfilter
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.filter);
  elseif hasleadfield
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.leadfield);
  else
    sourcemodel.inside = ft_inside_headmodel(sourcemodel.pos, headmodel);
  end
end

% convert to logical representation
sourcemodel = fixinside(sourcemodel);

if hasfilter && (fixedori || ~isequal(weightnorm, 'no'))
  ft_warning('with precomputed spatial filters a fixed orientation constraint or weight normalisation options are not applied');
end

% keep the original details on inside and outside positions
originside = sourcemodel.inside;
origpos    = sourcemodel.pos;

% select only the dipole positions inside the brain for scanning
sourcemodel.pos    = sourcemodel.pos(originside,:);
sourcemodel.inside = true(size(sourcemodel.pos,1),1);

if hasmom
  sourcemodel.mom = sourcemodel.mom(:,originside);
end

if hasfilter
  ft_info('using precomputed filters\n');
  sourcemodel.filter = sourcemodel.filter(originside);
elseif hasleadfield
  ft_info('using precomputed leadfields\n');
  sourcemodel.leadfield = sourcemodel.leadfield(originside);
else
  ft_info('computing forward model on the fly\n');
end

if hassubspace
  ft_info('using subspace projection\n');
  sourcemodel.subspace = sourcemodel.subspace(originside);
end

% dics has the following sub-methods, which depend on the function input arguments
% power only, cortico-muscular coherence and cortico-cortical coherence
if ~isempty(Cr) && ~isempty(Pr) && isempty(refdip)
  % compute cortico-muscular coherence, using reference cross spectral density
  submethod = 'dics_refchan';
elseif isempty(Cr) && isempty(Pr) && ~isempty(refdip)
  % compute cortico-cortical coherence with a dipole at the reference position
  submethod = 'dics_refdip';
  % check for forbidden options in combination with this submethod
  if hassubspace || ~isempty(subspace)
    ft_error('subspace projections are not supported for beaming cortico-cortical coherence');
  end
  if fixedori
    ft_error('fixed orientations are not supported for beaming cortico-cortical coherence');
  end
elseif isempty(Cr) && isempty(Pr) && isempty(refdip)
  % only compute power of a dipole at the grid positions
  submethod = 'dics_power';
else
  ft_error('invalid combination of input arguments for dics');
end

% it is difficult to give a quantitative estimate of lambda, therefore also
% support relative (percentage) measure that can be specified as string (e.g. '10%')
% the converted value needs to be passed on to ft_inv
lambda = ft_getopt(invopt, 'lambda');
if ~isempty(lambda) && ischar(lambda) && lambda(end)=='%'
  ratio  = sscanf(lambda, '%f%%');
  ratio  = ratio/100;
  lambda = ratio * trace(C)/size(C,1);
  if ~hassubspace
    invopt = ft_setopt(invopt, 'lambda', lambda);
  end
end

if projectnoise || strcmp(weightnorm, 'nai')
  % estimate the noise level in the covariance matrix by the smallest (non-zero) singular value
  % always needed for the NAI weight normalization case
  noise = svd(C);
  noise = noise(rank(C));
  % estimated noise floor is equal to or higher than lambda
  noise = max(noise, lambda);
end

% the inverse of the cross-spectral density matrix only has to be computed once for all dipoles
if hassubspace
  ft_notice('using source-specific subspace projection\n');
  % remember the original data prior to the voxel dependent subspace projection
  dat_pre_subspace = dat;
  C_pre_subspace = C;
  if strcmp(submethod, 'dics_refchan')
    Cr_pre_subspace = Cr;
    Pr_pre_subspace = Pr;
  end
elseif ~isempty(subspace)
  ft_notice('using data-specific subspace projection\n');
  % TODO implement an "eigenspace beamformer" as described in Sekihara et al. 2002 in HBM
  if numel(subspace)==1
    % interpret this as a truncation of the eigenvalue-spectrum
    % if <1 it is a fraction of the largest eigenvalue
    % if >=1 it is the number of largest eigenvalues
    dat_pre_subspace = dat;
    C_pre_subspace  = C;
    [u, s, v] = svd(real(C));
    if subspace<1
      subspace  = find(diag(s)./s(1,1) > subspace, 1, 'last');
    end

    C       = s(1:subspace,1:subspace);
    % this is equivalent to subspace*C*subspace' but behaves well numerically
    % by construction.
    invC     = diag(1./diag(C + lambda * eye(size(C))));
    subspace = u(:,1:subspace)';
    if ~isempty(dat), dat = subspace*dat; end

    if strcmp(submethod, 'dics_refchan')
      Cr = subspace*Cr;
    end

  else
    C_pre_subspace  = C;
    C               = subspace*C*subspace';
    % here the subspace can be different from the singular vectors of C, so we
    % have to do the sandwiching as opposed to line 256
    if strcmp(realfilter, 'yes')
      invC = ft_inv(real(C), invopt{:});
    else
      invC = ft_inv(C, invopt{:});
    end

    if strcmp(submethod, 'dics_refchan')
      Cr = subspace*Cr;
    end
  end
else
  if strcmp(realfilter, 'yes')
    % the filter is computed using only the leadfield and the inverse covariance or CSD matrix
    % therefore using the real-valued part of the CSD matrix here ensures a real-valued filter
    invC = ft_inv(real(C), invopt{:});
  else
    invC = ft_inv(C, invopt{:});
  end
end

if ~hassubspace
  % compute the square of invC, which might be needed for unitnoisegain or NAI constraint
  invC_squared = invC^2;
end

if strcmp(submethod, 'dics_refdip')
  % handle the reference dipole, this needs to be done only once, since in this implementation
  % no fancy pairwise dipole source model is assumed
  if isstruct(refdip) && isfield(refdip, 'filter') % check if precomputed filter is present
    assert(iscell(refdip.filter) && numel(refdip.filter)==1);
    filt1 = refdip.filter{1};
  else
    if ~isstruct(refdip)
      % this is very old behavior: refdip is assumed to be a 1x3 position vector
      refpos     = refdip; clear refdip;
      refdip.pos = refpos;
    end
    % now refdip is always a struct
    if isfield(refdip, 'leadfield') % check if precomputed leadfield is present
      assert(iscell(refdip.leadfield) && numel(refdip.leadfield)==1);
      lf1 = refdip.leadfield{1};
    elseif isfield(refdip, 'pos')
      assert(isnumeric(refdip.pos) && numel(refdip.pos)==3);
      lf1 = ft_compute_leadfield(refdip.pos, sens, headmodel, leadfieldopt{:});
    end
    if isfield(refdip, 'mom') % check for fixed orientation
      lf1 = lf1 * refdip.mom(:);
    end

    % construct the spatial filter
    switch weightnorm
      case 'nai'
        % Van Veen's Neural Activity Index
        ft_error('vector version of nai weight normalization is not implemented');
      case 'unitnoisegain'
        % filt*filt' = I
        % Unit-noise gain minimum variance (aka Borgiotti-Kaplan) beamformer
        denom = pinv(lf1' * invC * lf1);
        gamma = denom * (lf1' * invC_squared * lf1) * denom;

        % compute the spatial filter, as per eqn. 4.85
        filt1 = diag(1./sqrt(diag(gamma))) * denom * lf1' * invC;
      case 'arraygain'
        % filt*lf = ||lf||, applies to scalar leadfield, and to one of the possibilities of the vector version, eqn. 4.75
        lfn   = lf1./norm(lf1);
        filt1 = pinv(lfn' * invC * lfn) * lfn' * invC; % S&N eqn. 4.09 (scalar version), and eqn. 4.75 (vector version)
      case {'unitgain' 'no'}
        % this is the 'standard' unit gain constraint spatial filter: filt*lf=I, applies both to vector and scalar leadfields
        filt1 = pinv(lf1' * invC * lf1) * lf1' * invC; % Gross eqn. 3 & van Veen eqn. 23, use PINV/SVD to cover rank deficient leadfield
      otherwise
        ft_error('unsupported option for weightnorm');
    end
  end

  % compute the power of the reference dipole location
  if powlambda1
    Pr = lambda1(filt1 * C * ctranspose(filt1));      % compute the power at the first dipole location, Gross eqn. 8
  elseif powtrace
    Pr = real(trace(filt1 * C * ctranspose(filt1)));  % compute the power at the first dipole location
  end
end % if dics_refdip

% start the scanning with the proper metric
ft_progress('init', feedback, 'scanning grid');
for i=1:size(sourcemodel.pos,1)
  ft_progress(i/size(sourcemodel.pos,1), 'scanning grid %d/%d\n', i, size(sourcemodel.pos,1));

  if hasfilter
    % precomputed filter is provided, the leadfield is not needed
    filt = sourcemodel.filter{i};

  else
    if hasleadfield && hasmom && size(sourcemodel.mom, 1)==size(sourcemodel.leadfield{i}, 2)
      % reuse the leadfield that was previously computed and project
      lf = sourcemodel.leadfield{i} * sourcemodel.mom(:,i);
    elseif  hasleadfield &&  hasmom
      % reuse the leadfield that was previously computed but don't project
      lf = sourcemodel.leadfield{i};
    elseif  hasleadfield && ~hasmom
      % reuse the leadfield that was previously computed
      lf = sourcemodel.leadfield{i};
    elseif ~hasleadfield &&  hasmom
      % compute the leadfield for a fixed dipole orientation
      lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:}) * sourcemodel.mom(:,i);
    else
      % compute the leadfield
      lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:});
    end

    if hassubspace
      % do subspace projection of the forward model
      lforig = lf;
      lf = sourcemodel.subspace{i} * lf;
      % the cross-spectral density becomes voxel dependent due to the projection
      C    = sourcemodel.subspace{i} * C_pre_subspace * sourcemodel.subspace{i}';
      if strcmp(realfilter, 'yes')
        invC = ft_inv(sourcemodel.subspace{i} * real(C_pre_subspace) * sourcemodel.subspace{i}', invopt{:});
      else
        invC = ft_inv(sourcemodel.subspace{i} * C_pre_subspace * sourcemodel.subspace{i}', invopt{:});
      end
    elseif ~isempty(subspace)
      % do subspace projection of the forward model only
      lforig = lf;
      lf     = subspace * lf;

      % according to Kensuke's paper, the eigenspace bf boils down to projecting
      % the 'traditional' filter onto the subspace spanned by the first k eigenvectors
      % [u,s,v] = svd(C); filt = ESES*filt; ESES = u(:,1:k)*u(:,1:k)';
      % however, even though it seems that the shape of the filter is identical to
      % the shape it is obtained with the following code, the w*lf=I does not hold.
    end

    if fixedori
      switch(weightnorm)
        case {'unitnoisegain','nai'}
          % optimal orientation calculation for unit-noise gain beamformer,
          % (also applies nai weightnorm constraint), based on equation 4.47 from Sekihara & Nagarajan (2008)
          % the following is a reformulation of the generalized eigenproblem
          [v, d]        = eig(pinv(lf' * invC_squared *lf)*(lf' * invC *lf));
          [d, iv]       = sort(diag(d), 'descend');
          unitnoiseori  = v(:,iv(1));
          lf            = lf * unitnoiseori;
          estimate.ori{i} = unitnoiseori;
          estimate.eta(i) = d(1)./d(2); % ratio between largest and second largest eigenvalues
          if hassubspace, lforig = lforig * unitnoiseori; end

        case 'arraygain'
          % optimal orientation calculation for array-gain beamformer, Sekihara & Nagarajan eqn. 4.44
          [v, d]        = eig(pinv(lf' * invC *lf)*(lf' * lf));
          [d, iv]       = sort(diag(d), 'descend');
          arraygainori  = v(:,iv(1));
          lf            = lf * arraygainori;
          estimate.ori{i} = arraygainori;
          estimate.eta(i) = d(1)./d(2); % ratio between largest and second largest eigenvalues
          if hassubspace, lforig = lforig * arraygainori; end

        otherwise
          % compute the leadfield for the optimal dipole orientation that maximizes spatial filter output
          % subsequently the leadfield for only that dipole orientation will be used for the final filter computation
          % filt = pinv(lf' * invC * lf) * lf' * invC;
          % [u, s, v] = svd(real(filt * C * ctranspose(filt)));
          % in this step the filter computation is not necessary, use the quick way to compute the voxel level covariance (cf. van Veen 1997)
          [u, s, v]     = svd(real(pinv(lf' * invC *lf)));
          maxpowori     = u(:,1);
          lf            = lf * maxpowori;
          estimate.ori{i} = maxpowori;
          estimate.eta(i) = s(1,1)./s(2,2); % ratio between the first and second singular values
          if hassubspace, lforig = lforig * maxpowori; end
      end
    end

    % construct the spatial filter
    switch weightnorm
      case 'nai'
        % Van Veen's Neural Activity Index
        % below equation is equivalent to following:
        % filt = pinv(lf' * invC * lf) * lf' * invC;
        % filt = filt/sqrt(noise*filt*filt');
        % the scaling term in the denominator is sqrt of projected noise, as per eqn. 2.67 of Sekihara & Nagarajan 2008 (S&N)
        if fixedori
          filt = pinv(sqrt(noise * lf' * invC_squared * lf)) * lf' *invC; % based on S&N eqn. 4.08
        else
          ft_error('vector version of nai weight normalization is not implemented');
        end
      case 'unitnoisegain'
        % filt*filt' = I
        % Unit-noise gain minimum variance (aka Borgiotti-Kaplan) beamformer
        % below equation is equivalent to following:
        % filt = pinv(lf' * invC * lf) * lf' * invC;
        % filt = filt/sqrt(filt*filt');
        if fixedori
          filt = pinv(sqrt(lf' * invC_squared * lf)) * lf' *invC; % S&N eqn. 4.15
        else
          % compute the matrix that is used for scaling of the filter's rows, as per eqn. 4.83
          denom = pinv(lf' * invC * lf);
          gamma = denom * (lf' * invC_squared * lf) * denom;

          % compute the spatial filter, as per eqn. 4.85
          filt = diag(1./sqrt(diag(gamma))) * denom * lf' * invC;
        end
      case 'arraygain'
        % filt*lf = ||lf||, applies to scalar leadfield, and to one of the possibilities of the vector version, eqn. 4.75
        lfn  = lf./norm(lf);
        filt = pinv(lfn' * invC * lfn) * lfn' * invC; % S&N eqn. 4.09 (scalar version), and eqn. 4.75 (vector version)

      case {'unitgain' 'no'}
        % this is the 'standard' unit gain constraint spatial filter: filt*lf=I, applies both to vector and scalar leadfields
        filt = pinv(lf' * invC * lf) * lf' * invC; % Gross eqn. 3 & van Veen eqn. 23, use PINV/SVD to cover rank deficient leadfield

      otherwise
    end
  end

  cfilt = filt * C * ctranspose(filt); % Gross eqn. 4 and 5
  if powlambda1
    if size(cfilt,1) == 1
      % only 1 orientation, no need to do svd
      estimate.pow(i,1) = real(cfilt);
    else
      estimate.pow(i,1) = lambda1(cfilt);   % compute the power at the dipole location, Gross eqn. 8
    end
  elseif powtrace
    estimate.pow(i,1) = real(trace(cfilt)); % compute the power at the dipole location
  end

  if strcmp(submethod, 'dics_refchan')
    % compute the cross-spectrum between the dipole and the reference channel
    cfilt = filt*Cr; % Gross eqn. 6, replaces cfilt computed above
    if powlambda1
      coh = lambda1(cfilt)^2 / (estimate.pow(i,1) * Pr); % Gross eqn. 9
    elseif powtrace
      coh = norm(cfilt)^2 / (estimate.pow(i,1) * Pr);
    end
    estimate.coh(i,1) = coh;
  elseif strcmp(submethod, 'dics_refdip')
    % compute the cross-spectrum between the reference dipole and the dipole
    cfilt = filt1 * C * ctranspose(filt); % Gross eqn. 6, replaces cfilt computed above
    if powlambda1
      coh = lambda1(cfilt)^2 / (estimate.pow(i,1) * Pr); % Gross eqn. 9
    elseif powtrace
      coh = norm(cfilt)^2 / (estimate.pow(i,1) * Pr);
    end
    estimate.coh(i,1) = coh;
  end

  if keepcsd
    % keep the dipole's csd or the dipole-ref cross-spectrum in the output
    estimate.csd{i,1} = cfilt;
  end

  if projectnoise
    if powlambda1
      estimate.noise(i,1) = noise * lambda1(filt * ctranspose(filt));
    elseif powtrace
      estimate.noise(i,1) = noise * real(trace(filt * ctranspose(filt)));
    end
    if keepcsd
      estimate.noisecsd{i,1} = noise * filt * ctranspose(filt);
    end
  end

  if keepfilter
    if hassubspace
      estimate.filter{i,1} = filt*sourcemodel.subspace{i};
    elseif ~isempty(subspace)
      estimate.filter{i,1} = filt*subspace; % FIXME should this be subspace, or pinv(subspace)?
    else
      estimate.filter{i,1} = filt;
    end
  end

  if keepleadfield
    if ~isempty(subspace)
      estimate.leadfield{i,1} = lforig;
    else
      estimate.leadfield{i,1} = lf;
    end
  end

end % for each dipole position
ft_progress('close');

% reassign the estimated values over the inside and outside grid positions
estimate.inside  = originside;
estimate.pos     = origpos;

fnames_cell   = {'leadfield' 'filter' 'ori' 'csd' 'noisecsd' 'subspace'};
for k = 1:numel(fnames_cell)
  if isfield(estimate, fnames_cell{k})
    estimate.(fnames_cell{k})( originside) = estimate.(fnames_cell{k});
    estimate.(fnames_cell{k})(~originside) = {[]};
  end
end

fnames_scalar = {'pow' 'noise' 'eta' 'coh'};
for k = 1:numel(fnames_scalar)
  if isfield(estimate, fnames_scalar{k})
    estimate.(fnames_scalar{k})( originside) = estimate.(fnames_scalar{k});
    estimate.(fnames_scalar{k})(~originside) = nan;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to obtain the largest singular value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s, ori] = lambda1(x)
% determine the largest singular value, which corresponds to the power along the dominant direction
[u, s, v] = svd(x);
s   = s(1);
ori = u(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute the pseudo inverse. This is the same as the
% standard MATLAB function, except that the default tolerance is twice as
% high.
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision$  $Date: 2009/06/17 13:40:37 $
%   default tolerance increased by factor 2 (Robert Oostenveld, 7 Feb 2004)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = pinv(A,varargin)
[m,n] = size(A);
if n > m
  X = pinv(A',varargin{:})';
else
  [U,S,V] = svd(A,0);
  if m > 1, s = diag(S);
  elseif m == 1, s = S(1);
  else s = 0;
  end
  if nargin == 2
    tol = varargin{1};
  else
    tol = 10 * max(m,n) * max(s) * eps;
  end
  r = sum(s > tol);
  if (r == 0)
    X = zeros(size(A'),class(A));
  else
    s = diag(ones(r,1)./s(1:r));
    X = V(:,1:r)*s*U(:,1:r)';
  end
end


text:
function fwseg_dist= comp_fwseg(cleanFile, enhancedFile);

% ----------------------------------------------------------------------
%      Frequency weighted SNRseg Objective Speech Quality Measure
%
%   This function implements the frequency-weighted SNRseg measure [1]
%   using a different weighting function, the clean spectrum.
%
%   Usage:  fwSNRseg=comp_fwseg(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         fwSNRseg      - computed frequency weighted SNRseg in dB
%
%         Note that large numbers of fwSNRseg are better.
%
%  Example call:  fwSNRseg =comp_fwseg('sp04.wav','enhanced.wav')
%
%
%  References:
%   [1]  Tribolet, J., Noll, P., McDermott, B., and Crochiere, R. E. (1978).
%        A study of complexity and quality of speech waveform coders. Proc.
%        IEEE Int. Conf. Acoust. , Speech, Signal Processing, 586-590.
%
%   Author: Philipos C. Loizou
%  (critical-band filtering routines were written by Bryan Pellom & John Hansen)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
% ----------------------------------------------------------------------

if nargin~=2
    fprintf('USAGE: fwSNRseg=comp_fwseg(cleanFile.wav, enhancedFile.wav)\n');
    fprintf('For more help, type: help comp_fwseg\n\n');
    return;
end

[data1, Srate1, Nbits1]= wavread(cleanFile);
[data2, Srate2, Nbits2]= wavread(enhancedFile);
if ( Srate1~= Srate2) | ( Nbits1~= Nbits2)
    error( 'The two files do not match!\n');
end

len= min( length( data1), length( data2));
data1= data1( 1: len)+eps;
data2= data2( 1: len)+eps;

wss_dist_vec= fwseg( data1, data2,Srate1);

fwseg_dist=mean(wss_dist_vec);

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

function distortion = fwseg(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Files  must have same length.');
  return
end

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); 	   % window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
max_freq    = sample_rate/2;	   % maximum bandwidth
num_crit    = 25;		   % number of critical bands
USE_25=1;
n_fft       = 2^nextpow2(2*winlength);
n_fftby2    = n_fft/2;		   % FFT size/2
gamma=0.2;  % power exponent

% ----------------------------------------------------------------------
% Critical Band Filter Definitions (Center Frequency and Bandwidths in Hz)
% ----------------------------------------------------------------------

cent_freq(1)  = 50.0000;   bandwidth(1)  = 70.0000;
cent_freq(2)  = 120.000;   bandwidth(2)  = 70.0000;
cent_freq(3)  = 190.000;   bandwidth(3)  = 70.0000;
cent_freq(4)  = 260.000;   bandwidth(4)  = 70.0000;
cent_freq(5)  = 330.000;   bandwidth(5)  = 70.0000;
cent_freq(6)  = 400.000;   bandwidth(6)  = 70.0000;
cent_freq(7)  = 470.000;   bandwidth(7)  = 70.0000;
cent_freq(8)  = 540.000;   bandwidth(8)  = 77.3724;
cent_freq(9)  = 617.372;   bandwidth(9)  = 86.0056;
cent_freq(10) = 703.378;   bandwidth(10) = 95.3398;
cent_freq(11) = 798.717;   bandwidth(11) = 105.411;
cent_freq(12) = 904.128;   bandwidth(12) = 116.256;
cent_freq(13) = 1020.38;   bandwidth(13) = 127.914;
cent_freq(14) = 1148.30;   bandwidth(14) = 140.423;
cent_freq(15) = 1288.72;   bandwidth(15) = 153.823;
cent_freq(16) = 1442.54;   bandwidth(16) = 168.154;
cent_freq(17) = 1610.70;   bandwidth(17) = 183.457;
cent_freq(18) = 1794.16;   bandwidth(18) = 199.776;
cent_freq(19) = 1993.93;   bandwidth(19) = 217.153;
cent_freq(20) = 2211.08;   bandwidth(20) = 235.631;
cent_freq(21) = 2446.71;   bandwidth(21) = 255.255;
cent_freq(22) = 2701.97;   bandwidth(22) = 276.072;
cent_freq(23) = 2978.04;   bandwidth(23) = 298.126;
cent_freq(24) = 3276.17;   bandwidth(24) = 321.465;
cent_freq(25) = 3597.63;   bandwidth(25) = 346.136;

W=[  % articulation index weights
0.003
0.003
0.003
0.007
0.010
0.016
0.016
0.017
0.017
0.022
0.027
0.028
0.030
0.032
0.034
0.035
0.037
0.036
0.036
0.033
0.030
0.029
0.027
0.026
0.026];

W=W';

if USE_25==0  % use 13 bands
    % ----- lump adjacent filters together ----------------
    k=2;
    cent_freq2(1)=cent_freq(1);
    bandwidth2(1)=bandwidth(1)+bandwidth(2);
    W2(1)=W(1);
    for i=2:13
        cent_freq2(i)=cent_freq2(i-1)+bandwidth2(i-1);
        bandwidth2(i)=bandwidth(k)+bandwidth(k+1);
        W2(i)=0.5*(W(k)+W(k+1));
        k=k+2;
    end

    sumW=sum(W2);
    bw_min      = bandwidth2 (1);	   % minimum critical bandwidth
else
    sumW=sum(W);
    bw_min=bandwidth(1);
end

% ----------------------------------------------------------------------
% Set up the critical band filters.  Note here that Gaussianly shaped
% filters are used.  Also, the sum of the filter weights are equivalent
% for each critical band filter.  Filter less than -30 dB and set to
% zero.
% ----------------------------------------------------------------------

min_factor = exp (-30.0 / (2.0 * 2.303));       % -30 dB point of filter
if USE_25==0

    num_crit=length(cent_freq2);

    for i = 1:num_crit
        f0 = (cent_freq2 (i) / max_freq) * (n_fftby2);
        all_f0(i) = floor(f0);
        bw = (bandwidth2 (i) / max_freq) * (n_fftby2);
        norm_factor = log(bw_min) - log(bandwidth2(i));
        j = 0:1:n_fftby2-1;
        crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
        crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
    end

else
    for i = 1:num_crit
        f0 = (cent_freq (i) / max_freq) * (n_fftby2);
        all_f0(i) = floor(f0);
        bw = (bandwidth (i) / max_freq) * (n_fftby2);
        norm_factor = log(bw_min) - log(bandwidth(i));
        j = 0:1:n_fftby2-1;
        crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
        crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
    end
end

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1:num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the magnitude Spectrum of Clean and Processed
   % ----------------------------------------------------------

       clean_spec     = abs(fft(clean_frame,n_fft));
       processed_spec = abs(fft(processed_frame,n_fft));

    % normalize spectra to have area of one
    %
    clean_spec=clean_spec/sum(clean_spec(1:n_fftby2));
    processed_spec=processed_spec/sum(processed_spec(1:n_fftby2));

   % ----------------------------------------------------------
   % (3) Compute Filterbank Output Energies
   % ----------------------------------------------------------

   clean_energy=zeros(1,num_crit);
   processed_energy=zeros(1,num_crit);
   error_energy=zeros(1,num_crit);
   W_freq=zeros(1,num_crit);

   for i = 1:num_crit
      clean_energy(i) = sum(clean_spec(1:n_fftby2) ...
                         	.*crit_filter(i,:)');
      processed_energy(i) = sum(processed_spec(1:n_fftby2) ...
          .*crit_filter(i,:)');

        error_energy(i)=max((clean_energy(i)-processed_energy(i))^2,eps);
        W_freq(i)=(clean_energy(i))^gamma;

   end
   SNRlog=10*log10((clean_energy.^2)./error_energy);

   fwSNR=sum(W_freq.*SNRlog)/sum(W_freq);

   distortion(frame_count)=min(max(fwSNR,-10),35);

   start = start + skiprate;

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [x,y,z,nx,ny,nz,red,green,blue] = import_3D_model(FileName,prefer_mat_file)

if(~exist(FileName,'file'))
    error('Specified file does not exist:%s',FileName);
end

if(~exist('prefer_mat_file','var'))
    prefer_mat_file = 0;
end

[pathstr, name, ext] = fileparts(FileName);

if(strcmpi(ext,'.ply'))

    % check to see if there's a .mat file corresponding to the .ply file
    ply_info = dir(FileName);
    ply_date = datenum(ply_info.date);

    [pathstr, name, ext] = fileparts(FileName);
    mat_file = [pathstr, '/' ,name, '.mat'];
    mat_info = dir(mat_file);

    mat_file_is_there = 0;
    if(length(mat_info) ~= 0) %if there is a mat file (and it is newer than the ply file)
        mat_date = datenum(mat_info.date);
        if  mat_date>ply_date
            mat_file_is_there = 1;
        end
    end

    %now try to open the right file
    if(prefer_mat_file)
        if(mat_file_is_there)
            load(mat_file);
        else
            data = nplyRead(FileName);
            save([pathstr, '/' ,name, '.mat'],'data');
            msgbox('A .mat version of your .ply file is now saved in the same location. You can use that instead!','replace');
        end
    else
        data = nplyRead(FileName);
    end

elseif(strcmpi(ext,'.mat'))
    load(FileName);
else
    error('Invalid file extension');
end

x = data.vertex.x;
y = data.vertex.y;
z = data.vertex.z;

try
    nx = data.vertex.nx;
    ny = data.vertex.ny;
    nz = data.vertex.nz;
catch
    try
        nx = data.vertex.vsfm_cnx;
        ny = data.vertex.vsfm_cny;
        nz = data.vertex.vsfm_cnz;
    catch
        nx = [];
        ny = [];
        nz = [];
    end
end

try
    red = data.vertex.r/255;
    green = data.vertex.g/255;
    blue = data.vertex.b/255;
catch
    try
        red = data.vertex.red/255;
        green = data.vertex.green/255;
        blue = data.vertex.blue/255;
    catch
        try
            red = data.vertex.diffuse_red/255;
            green = data.vertex.diffuse_green/255;
            blue = data.vertex.diffuse_blue/255;
        catch
            red = .2;
            green = .2;
            blue = .2;
        end
    end
end

if(nargout == 1)
    model.x = x;
    model.y = y;
    model.z = z;
    model.nx = nx;
    model.ny = ny;
    model.nz = nz;
    model.r = red;
    model.g = green;
    model.b = blue;

    if(isfield(data,'face'))
        model.face = data.face;
    end

    x = model;
end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function index = randIndex(maxIndex,len)
%INDEX = RANDINDEX(MAXINDEX,LEN)
%   randomly, non-repeatedly select LEN integers from 1:MAXINDEX

if len > maxIndex
	index = [];
	return
end

index = zeros(1,len);
available = 1:maxIndex;
rs = ceil(rand(1,len).*(maxIndex:-1:maxIndex-len+1));
for p = 1:len
	while rs(p) == 0
		rs(p) = ceil(rand(1)*(maxIndex-p+1));
	end
	index(p) = available(rs(p));
	available(rs(p)) = [];
end

author: Sable

repo: mcbench-benchmarks


text:
function [bs] = A2bs(A)
% Convert length from angstroms to beard-seconds.
% Chad A. Greene 2012
bs = A/100;


text:
function [ x, know ] = p02_sol ( m, know )

%*****************************************************************************80
%
%% P02_SOL returns known solutions for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution.
%
  if ( know == 0 )
    know = 1;
    x = [ ...
      0.390500591228663, ...
      0.392051909813608, ...
      0.393601661544812, ...
      0.395149843840982 ]';
  else
    know = 0;
    x = zeros ( m, 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Iphi,SigmaPhi,deltaMuPhi,suffStat] = VBA_Iphi_binomial(phi,y,posterior,suffStat,dim,u,options)
% Gauss-Newton update of the observation parameters, for binomial data

if options.DisplayWin % Display progress
    if isequal(options.g_fname,@VBA_odeLim)
        STR = 'VB Gauss-Newton on observation/evolution parameters... ';
    else
        STR = 'VB Gauss-Newton on observation parameters... ';
    end
    set(options.display.hm(1),'string',STR);
    set(options.display.hm(2),'string','0%');
    drawnow
end

%  Look-up which evolution parameter to update
indIn = options.params2update.phi;

% Preallocate intermediate variables
Q = options.priors.SigmaPhi(indIn,indIn);
iQ = VBA_inv(Q);
muPhi0 = options.priors.muPhi;
Phi = muPhi0;
Phi(indIn) = phi;
dphi0 = muPhi0-Phi;
ddydphi = 0;
dy = zeros(dim.p,dim.n_t);
vy = zeros(dim.p,dim.n_t);
gx = zeros(dim.p,dim.n_t);
d2gdx2 = zeros(dim.n_phi,dim.n_phi);
logL = 0;
if isequal(options.g_fname,@VBA_odeLim)
    clear VBA_odeLim
    muX = zeros(options.inG.old.dim.n,dim.n_t);
    SigmaX = cell(dim.n_t,1);
end
div = 0;

%--- Loop over time series ---%
for t=1:dim.n_t

    % evaluate observation function at current mode
    [gx(:,t),dG_dX,dG_dPhi] = VBA_evalFun('g',posterior.muX(:,t),Phi,u(:,t),options,dim,t);

    % fix numerical instabilities
    gx(:,t) = VBA_finiteBinomial (gx(:,t));

    % store states dynamics if ODE mode
    if isequal(options.g_fname,@VBA_odeLim)
        % get sufficient statistics of the hidden states from unused i/o in
        % VBA_evalFun.
        muX(:,t) = dG_dX.xt;
        SigmaX{t} = dG_dX.dx'*posterior.SigmaPhi*dG_dX.dx;
    end

    % predicted variance over binomial data
    vy(:,t) = gx(:,t).*(1-gx(:,t));

    % remove irregular trials
    yin = find(~options.isYout(:,t));

    % accumulate log-likelihood
    logL = logL + y(yin,t)'*log(gx(yin,t)) + (1-y(yin,t))'*log(1-gx(yin,t));

    % prediction error
    dy(yin,t) = y(yin,t) - gx(yin,t);

    % gradient and Hessian
    ddydphi = ddydphi + dG_dPhi(:,yin)*(dy(yin,t)./vy(yin,t));
    tmp = y(yin,t)./gx(yin,t).^2 - (y(yin,t)-1)./(1-gx(yin,t)).^2;
    d2gdx2 = d2gdx2 + dG_dPhi(:,yin)*diag(tmp)*dG_dPhi(:,yin)';

    % Display progress
    if mod(t,dim.n_t./10) < 1
        if options.DisplayWin
            set(options.display.hm(2),'string',[num2str(floor(100*t/dim.n_t)),'%']);
            drawnow
        end
    end

    % Accelerate divergent update
    if VBA_isWeird ({dy, dG_dPhi, dG_dX})
        div = 1;
        break
    end

end

% Display progress
if options.DisplayWin
    set(options.display.hm(2),'string','OK');
    drawnow
end

% posterior covariance matrix
iSigmaPhi = iQ + d2gdx2(indIn,indIn);
SigmaPhi = VBA_inv(iSigmaPhi);

% mode
tmp = iQ*dphi0(indIn) + ddydphi(indIn);
deltaMuPhi = SigmaPhi*tmp;

% variational energy
Iphi = -0.5.*dphi0(indIn)'*iQ*dphi0(indIn) + logL;
if VBA_isWeird ({Iphi, SigmaPhi}) || div
    Iphi = -Inf;
end

% update sufficient statistics
suffStat.gx = gx;
suffStat.dy = dy;
suffStat.logL = logL;
suffStat.vy = vy;
suffStat.dphi = dphi0;
if isequal(options.g_fname,@VBA_odeLim)
    suffStat.muX = muX;
    suffStat.SigmaX = SigmaX;
end
suffStat.div = div;

author: MBB-team

repo: VBA-toolbox


text:
% this code is revised based on ILSVRC 2013 (http://www.image-net.org/challenges/LSVRC/2013/)

function  top_recall = top_recall_Phrase(Nre, tuple_confs_cell, tuple_labels_cell, sub_bboxes_cell, obj_bboxes_cell)

load('gt.mat','gt_tuple_label','gt_obj_bboxes','gt_sub_bboxes');

%num_imgs = length(gt_tuple_label);
num_imgs = 1000;
for i=1:num_imgs
    [tuple_confs_cell{i}, ind] = sort(tuple_confs_cell{i},'descend');
    if length(ind) >= Nre
        tuple_confs_cell{i} = tuple_confs_cell{i}(1:Nre);
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind(1:Nre),:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind(1:Nre),:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind(1:Nre),:);
    else
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind,:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind,:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind,:);
    end
end

num_pos_tuple = 0;
for ii = 1 : num_imgs
    num_pos_tuple = num_pos_tuple + size(gt_tuple_label{ii},1);
end

tp_cell = cell(1,num_imgs);
fp_cell = cell(1,num_imgs);

gt_thr = 0.5;
% iterate over images
for i=1:num_imgs

    gt_tupLabel = gt_tuple_label{i};

    if ~isempty(gt_obj_bboxes{i})
        gt_box_entity = [min(gt_obj_bboxes{i}(:,1:2),gt_sub_bboxes{i}(:,1:2)),max(gt_obj_bboxes{i}(:,3:4),gt_sub_bboxes{i}(:,3:4))];
    else
        gt_box_entity = [];
    end

    num_gt_tuple = size(gt_tupLabel,1);
    gt_detected = zeros(1,num_gt_tuple);

    labels = tuple_labels_cell{i};

    boxObj = obj_bboxes_cell{i};
    boxSub = sub_bboxes_cell{i};
    if ~isempty(boxObj)
        box_entity_our  = [min(boxObj(:,1:2), boxSub(:,1:2)), max(boxObj(:,3:4), boxSub(:,3:4))];
    else
        box_entity_our  = [];
    end

    num_obj = size(labels,1);
    tp = zeros(1,num_obj);
    fp = zeros(1,num_obj);
    for j=1:num_obj

        bbO = box_entity_our(j,:);
        ovmax = -inf;
        kmax = -1;

        for k=1:num_gt_tuple
            if norm(labels(j,:) - gt_tupLabel(k,:),2) ~= 0
                continue;
            end
            if gt_detected(k) > 0
                continue;
            end

            bbgtO = gt_box_entity(k,:);

            biO=[max(bbO(1),bbgtO(1)) ; max(bbO(2),bbgtO(2)) ; min(bbO(3),bbgtO(3)) ; min(bbO(4),bbgtO(4))];
            iwO=biO(3)-biO(1)+1;
            ihO=biO(4)-biO(2)+1;

            if iwO>0 & ihO>0
                % compute overlap as area of intersection / area of union
                uaO=(bbO(3)-bbO(1)+1)*(bbO(4)-bbO(2)+1)+...
                   (bbgtO(3)-bbgtO(1)+1)*(bbgtO(4)-bbgtO(2)+1)-...
                   iwO*ihO;
                ov =iwO*ihO/uaO;

                % makes sure that this object is detected according
                % to its individual threshold
                if ov >= gt_thr && ov > ovmax
                    ovmax=ov;
                    kmax=k;
                end
            end
        end

        if kmax > 0
            tp(j) = 1;
            gt_detected(kmax) = 1;
        else
            fp(j) = 1;
        end
    end

    % put back into global vector
    tp_cell{i} = tp;
    fp_cell{i} = fp;

end

t = tic;
tp_all = [];
fp_all = [];
confs = [];
for ii = 1 : num_imgs
tp_all = [tp_all; tp_cell{ii}(:) ];
fp_all = [fp_all; fp_cell{ii}(:) ];
confs = [confs; tuple_confs_cell{ii}(:)];
end

[confs, ind] = sort(confs,'descend');
tp_all = tp_all(ind);
fp_all = fp_all(ind);

tp = cumsum(tp_all );
fp = cumsum(fp_all );
recall =(tp/num_pos_tuple);
top_recall = recall(end);

end

author: Prof-Lu-Cewu

repo: Visual-Relationship-Detection


text:
function [output_dir] = run_bu_experiment(bu_dir, verbose, varargin)

    if(isunix)
        executable = '"../../build/bin/FeatureExtraction"';
    else
        executable = '"../../x64/Release/FeatureExtraction.exe"';
    end

    output_dir = 'experiments/bu_out/';

    buFiles = dir([bu_dir '*.avi']);

    % Only outputing the pose (-pose)
    command = sprintf('%s -inroot "%s" -out_dir "%s" -fx 500 -fy 500 -cx 160 -cy 120 -pose -vis-track ', executable, bu_dir, output_dir);

    for i=1:numel(buFiles)
        inputFile = [buFiles(i).name];
        command = cat(2, command, sprintf(' -f "%s" ', inputFile));
    end

    if(verbose)
        command = cat(2, command, [' -tracked ' outputVideo]);
    end

    if(any(strcmp('model', varargin)))
        command = cat(2, command, [' -mloc "', varargin{find(strcmp('model', varargin))+1}, '"']);
    end

    if(isunix)
        unix(command, '-echo')
    else
        dos(command);
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function cnn = init_cnn_model(varargin)
% cnn = init_cnn_model
% Initialize a CNN with caffe
%
% Optional arguments
% net_file   network binary file
% def_file   network prototxt file
% use_gpu    set to false to use CPU (default: true)
% use_caffe  set to false to avoid using caffe (default: true)
%            useful for running on the cluster (must use cached pyramids!)

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the DeepPyramid code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

% ------------------------------------------------------------------------
% Options
ip = inputParser;

% network binary file
ip.addParamValue('net_file', ...
    './data/caffe_nets/ilsvrc_2012_train_iter_310k', ...
    @isstr);

% network prototxt file
ip.addParamValue('def_file', ...
    './model-defs/pyramid_cnn_output_conv5_scales_7_plane_1713.prototxt', ...
    @isstr);

% Set use_gpu to false to use the CPU
ip.addParamValue('use_gpu', true, @islogical);

% Set use_caffe to false to avoid using caffe
% (must be used in conjunction with cached features!)
ip.addParamValue('use_caffe', true, @islogical);

ip.parse(varargin{:});
opts = ip.Results;
% ------------------------------------------------------------------------

cnn.binary_file = opts.net_file;
cnn.definition_file = opts.def_file;
cnn.init_key = -1;
cnn.mu = get_channelwise_mean;
cnn.pyra.dimx = 1713;
cnn.pyra.dimy = 1713;
cnn.pyra.stride = 16;
cnn.pyra.num_levels = 7;
cnn.pyra.num_channels = 256;
cnn.pyra.scale_factor = sqrt(2);

if opts.use_caffe
  cnn.init_key = ...
      caffe('init', cnn.definition_file, cnn.binary_file);
  caffe('set_phase_test');
  if opts.use_gpu
    caffe('set_mode_gpu');
  else
    caffe('set_mode_cpu');
  end
end

% ------------------------------------------------------------------------
function mu = get_channelwise_mean()
% ------------------------------------------------------------------------
% load the ilsvrc image mean
data_mean_file = 'ilsvrc_2012_mean.mat';
assert(exist(data_mean_file, 'file') ~= 0);
% input size business isn't likley necessary, but we're doing it
% to be consistent with previous experiments
ld = load(data_mean_file);
mu = ld.image_mean; clear ld;
input_size = 227;
off = floor((size(mu,1) - input_size)/2)+1;
mu = mu(off:off+input_size-1, off:off+input_size-1, :);
mu = sum(sum(mu, 1), 2) / size(mu, 1) / size(mu, 2);

author: rbgirshick

repo: DeepPyramid


text:
function varargout = process_average_freq( varargin )
% PROCESS_AVERAGE_FREQ: For each file in input, compute the average of the different frequency bands.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Average frequency bands';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Average';
    sProcess.Index       = 303;
    sProcess.Description = '';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'timefreq'};
    sProcess.OutputTypes = {'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % === OVERWRITE
    sProcess.options.overwrite.Comment = 'Overwrite input files';
    sProcess.options.overwrite.Type    = 'checkbox';
    sProcess.options.overwrite.Value   = 0;
    sProcess.options.overwrite.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFile = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFile = [];
    isOverwrite = sProcess.options.overwrite.Value;
    % Load TF file
    TimefreqMat = in_bst_timefreq(sInput.FileName, 0);
    % Check for measure
    if strcmpi(TimefreqMat.Measure, 'none')
        bst_report('Error', sProcess, sInput, 'Cannot average complex values. Please apply a measure to the values before calling this function.');
        return;
    end
    % Average time-freq values
    TimefreqMat.TF = mean(TimefreqMat.TF, 3);
    % Remove the frequency information
    if iscell(TimefreqMat.Freqs)
        TimefreqMat.Freqs = {'Avg freq', TimefreqMat.Freqs{1,2}, TimefreqMat.Freqs{end,3}};
    else
        TimefreqMat.Freqs = TimefreqMat.Freqs(end);
    end
    % Do not keep the Std field in the output
    if isfield(TimefreqMat, 'Std') && ~isempty(TimefreqMat.Std)
        TimefreqMat.Std = [];
    end
    % TFmask: Keep only the values good everywhere
    if isfield(TimefreqMat, 'TFmask') && ~isempty(TimefreqMat.TFmask)
        TimefreqMat.TFmask = all(TimefreqMat.TFmask, 1);
    end

    % Add file tag
    TimefreqMat.Comment = [TimefreqMat.Comment, ' | avgfreq'];
    % Add history entry
    TimefreqMat = bst_history('add', TimefreqMat, 'avgfreq', 'Averaged all the frequency bands.');

    % Overwrite the input file
    if isOverwrite
        OutputFile = file_fullpath(sInput.FileName);
        bst_save(OutputFile, TimefreqMat, 'v6');
        % Reload study
        db_reload_studies(sInput.iStudy);
    % Save new file
    else
        % Output filename: add file tag
        OutputFile = strrep(file_fullpath(sInput.FileName), '.mat', '_avgfreq.mat');
        OutputFile = file_unique(OutputFile);
        % Save file
        bst_save(OutputFile, TimefreqMat, 'v6');
        % Add file to database structure
        db_add_data(sInput.iStudy, OutputFile, TimefreqMat);
    end
end

author: brainstorm-tools


text:
% StackExchange Mathematics Q3674241
% https://math.stackexchange.com/questions/3674241
% Solve argmina∥ax−y∥1 - Minimizer of the L1 Norm of the Difference of a Vectors
% References:
%   1.
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     14/05/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 42; %<! Set to 0 for Random

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numElements = 5; %<! Symmetric Matrix
numPts      = numElements * 1000;

%% Load / Generate Data

vX = randn(numElements, 1);
vY = randn(numElements, 1);

hObjFun = @(valA) sum(abs(valA * vX - vY)); %<! Scaled L1 Norm
hS = @(valA) sign(valA * vX - vY); %<! The 'vS' vector

%% Solution by CVX

solverString = 'CVX';

cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
% cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
% cvx_begin()
    % cvx_precision('best');
    variable valA(1, 1);
    minimize( norm(valA * vX - vY, 1) );
cvx_end

runTime = toc(hRunTime);

% vX = mX(:);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The ', solverString, ' Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(valA))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(valA), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

%% Solution by Minimum Search

solverString = 'Solution by 1D Search of the Minimum';

hRunTime = tic();

valA = fminsearch(hObjFun, mean(vY ./ vX));

runTime = toc(hRunTime);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(valA))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(valA), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

%% Solution by Root Search

solverString = 'Solution by 1D Search of the Roto of the Gradient';

hF = @(valA) sign(valA * vX - vY).' * vX;

hRunTime = tic();

valA = fzero(hF, mean(vY ./ vX));

runTime = toc(hRunTime);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(valA))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(valA), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

%% Solution by Analytic Form

solverString = 'Solution by Analytic Form';

hRunTime = tic();

valA = SolveScaledL1(vX, vY);

runTime = toc(hRunTime);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(valA))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(valA), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

%% Display Results

vA = sort(vY ./ vX, 'ascend');
vB = zeros(numElements, 1);
vC = zeros(numElements, 1);

for ii = 1:numElements
    vB(ii) = hS(vA(ii)).' * vX;
    vC(ii) = hObjFun(vA(ii));
end

vD = linspace(vA(1) - 1, vA(numElements) + 1, numPts);
vD = sort([vD(:); vA], 'ascend');
vE = zeros(numPts + numElements, 1);
vF = zeros(numPts + numElements, 1);

for ii = 1:(numPts + numElements)
    vE(ii) = hS(vD(ii)).' * vX;
    vF(ii) = hObjFun(vD(ii));
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);

hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vD, vE);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries = plot(vD, zeros(numPts + numElements, 1));
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries = plot(vA, vB);
set(hLineSeries, 'LineStyle', 'none', 'Marker', '*', 'MarkerSize', markerSizeLarge);
set(get(hAxes, 'Title'), 'String', ['Value of f''(a)'], 'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'a', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ {f}^{''} \left( a \right) $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
hLegend = ClickableLegend({['Values of All Points'], ['Zero Range'], ['Junction Points']});

if(generateFigures == ON)
    % saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);

hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vD, vF);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries = plot(vA, vC);
set(hLineSeries, 'LineStyle', 'none', 'Marker', '*', 'MarkerSize', markerSizeLarge);
set(get(hAxes, 'Title'), 'String', ['Value of f(a)'], 'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'a', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ f \left( a \right) $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
hLegend = ClickableLegend({['Values of All Points'], ['Junction Points']});

if(generateFigures == ON)
    % saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function dat = read_spike6mat_data(filename, varargin)

% read_spike6mat_data() - read Matlab files exported from Spike 6
%
% Usage:
%   >> header = read_spike6mat_data(filename, varargin);
%
% Inputs:
%   filename - [string] file name
%
% Optional inputs:
%   'begsample'      first sample to read
%   'endsample'      last sample to read
%   'chanindx'  -    list with channel indices to read
%   'header'    - FILEIO structure header
%
% Outputs:
%   dat    - data over the specified range
% _______________________________________________________________________
% Copyright (C) 2008 Institute of Neurology, UCL
% Vladimir Litvak

if nargin < 1
    help read_spike6mat_data;
    return;
end;

header    = ft_getopt(varargin, 'header');
begsample = ft_getopt(varargin, 'begsample');
endsample = ft_getopt(varargin, 'endsample');
chanindx  = ft_getopt(varargin, 'chanindx');

if isempty(header)
    header = read_spike6mat_header(filename);
end

if isempty(begsample), begsample = 1; end;
if isempty(endsample), endsample = header.nSamples; end;

try
    vars = struct2cell(load(filename));
catch
    ft_error('File not found or wrong format.');
end

if isempty(chanindx)
    chanindx = 1:numel(vars);
end

dat = zeros(length(chanindx), endsample-begsample+1);

for i = 1:length(chanindx)
    dat(i, :) = vars{chanindx(i)}.values(begsample:endsample);
end

author: spm


text:
function [ err_outline, err_no_outline ] = Run_OF_on_images(output_loc, database_root, varargin)

dataset_dirs = {};

if(any(strcmp(varargin, 'use_afw')))
    afw_loc = [database_root '/AFW/'];
    dataset_dirs = cat(1, dataset_dirs, afw_loc);
end
if(any(strcmp(varargin, 'use_lfpw')))
    lfpw_loc = [database_root 'lfpw/testset/'];
    dataset_dirs = cat(1, dataset_dirs, lfpw_loc);
end
if(any(strcmp(varargin, 'use_ibug')))
    ibug_loc = [database_root 'ibug/'];
    dataset_dirs = cat(1, dataset_dirs, ibug_loc);
end
if(any(strcmp(varargin, 'use_helen')))
    helen_loc = [database_root '/helen/testset/'];
    dataset_dirs = cat(1, dataset_dirs, helen_loc);
end

if(any(strcmp(varargin, 'verbose')))
    verbose = true;
else
    verbose = false;
end

if(isunix)
    executable = '"../../build/bin/FaceLandmarkImg"';
else
    executable = '"../../x64/Release/FaceLandmarkImg.exe"';
end

if(any(strcmp(varargin, 'model')))
    model = varargin{find(strcmp(varargin, 'model')) + 1};
else
    % the default model is the 68 point in the wild one
    model = '"model/main_wild.txt"';
end

if(any(strcmp(varargin, 'multi_view')))
    multi_view = varargin{find(strcmp(varargin, 'multi_view')) + 1};
else
    multi_view = 0;
end

command = sprintf('%s -mloc %s -multi_view %s -2Dfp ', executable, model, num2str(multi_view));

tic
for i=1:numel(dataset_dirs)

    command_c = sprintf('%s -fdir "%s" -bboxdir "%s" -out_dir "%s" -wild ',...
        command, dataset_dirs{i},  dataset_dirs{i}, output_loc);

    if(isunix)
        unix(command_c, '-echo');
    else
        dos(command_c);
    end

end
toc

%%

% Extract the error sizes
dirs = {[database_root '/AFW/'];
    [database_root '/ibug/'];
    [database_root '/helen/testset/'];
    [database_root 'lfpw/testset/'];};

landmark_dets = dir([output_loc '/*.csv']);

landmark_det_dir = [output_loc '/'];

num_imgs = size(landmark_dets,1);

labels = zeros(68,2,num_imgs);

shapes = zeros(68,2,num_imgs);

curr = 0;

% work out which columns in the csv file are relevant
tab = readtable([landmark_det_dir, landmark_dets(1).name]);
column_names = tab.Properties.VariableNames;
landmark_inds_x = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'x_'));
landmark_inds_y = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'y_'));

for i=1:numel(dirs)

    gt_labels = dir([dirs{i}, '*.pts']);

    for g=1:numel(gt_labels)
        curr = curr+1;

        gt_landmarks = dlmread([dirs{i}, gt_labels(g).name], ' ', 'A4..B71');
        [~, name, ~] = fileparts(gt_labels(g).name);
        % find the corresponding detection
        all_params  = dlmread([landmark_det_dir, name, '.csv'], ',', 1, 0);

        landmark_det = [all_params(landmark_inds_x); all_params(landmark_inds_y)]';

        labels(:,:,curr) = gt_landmarks;

        if(size(landmark_det,1) == 66)
            inds_66 = [[1:60],[62:64],[66:68]];
            shapes(inds_66,:,curr) = landmark_det;
        else
            shapes(:,:,curr) = landmark_det;
        end
    end

end

% Convert to correct format, so as to have same feature points in ground
% truth and detections
if(size(shapes,2) == 66 && size(labels,2) == 68)
    inds_66 = [[1:60],[62:64],[66:68]];

    labels = labels(inds_66,:,:);
    shapes = shapes(inds_66,:,:);
end

% Center the pixel, and convert to OCV format
labels = labels - 2.0;

err_outline = compute_error(labels, shapes);

labels_no_out = labels(18:end,:,:);
shapes_no_out = shapes(18:end,:,:);

err_no_outline = compute_error(labels_no_out, shapes_no_out);

%%

save([output_loc, 'res.mat'], 'labels', 'shapes', 'err_outline', 'err_no_outline');

end


text:
% -------------------------------------------------------------------------------------------------------------------
function save_crops(imdb_video,v_1,v_end, root_original, root_crops)
    % Extract and save crops from video v_1 (start from 1) to v_end (check num video in imdb)
	% e.g. save_crops(imdb_video, 1, 1000, '/path/to/original/ILSVRC15/', '/path/to/new/curated/ILSVRC15/')
% -------------------------------------------------------------------------------------------------------------------
    rootDataDir_src = [root_original '/Data/VID/train/'];
    rootDataDir_dest = [root_crops '/Data/VID/train/'];
    % sides of the crops for z and x saved on disk
    exemplar_size = 127;
    instance_size = 255;
    context_amount = 0.5;

    saved_crops = 0;
    for v=v_1:v_end
        valid_trackids = find(imdb_video.valid_trackids(:,v));
        for ti=1:numel(valid_trackids)
            valid_objects = imdb_video.valid_per_trackid{valid_trackids(ti),v};
            for o = 1:numel(valid_objects)
                obj = imdb_video.objects{v}{valid_objects(o)};
                assert(obj.valid)
                fprintf('%d %d: %s\n', v, obj.track_id, obj.frame_path);
                im = imread([rootDataDir_src obj.frame_path]);

                [im_crop_z, bbox_z, pad_z, im_crop_x, bbox_x, pad_x] = get_crops(im, obj, exemplar_size, instance_size, context_amount);

                root = [rootDataDir_dest strrep(obj.frame_path,'.JPEG','') '.' num2str(obj.track_id,'%02d')];
                pz = fopen([root '.pad.z.txt'],'w');
                px = fopen([root '.pad.x.txt'],'w');
                fprintf(pz,'%.2f,%.2f,%.2f,%.2f\n', pad_z(1),pad_z(2),pad_z(3),pad_z(4));
                fprintf(px,'%.2f,%.2f,%.2f,%.2f\n', pad_x(1),pad_x(2),pad_x(3),pad_x(4));
                fclose(pz);
                fclose(px);
                imwrite(im_crop_z, [root '.crop.z.jpg'], 'Quality', 90);
                imwrite(im_crop_x, [root '.crop.x.jpg'], 'Quality', 90);

                saved_crops = saved_crops+1;
            end
        end
    end

    fprintf('\n:: SAVED %d crops ::\n', saved_crops);

end

% -------------------------------------------------------------------------------------------------------------------
function [im_crop_z, bbox_z, pad_z, im_crop_x, bbox_x, pad_x] = get_crops(im, object, size_z, size_x, context_amount)
% -------------------------------------------------------------------------------------------------------------------
    %% Get exemplar sample
    % take bbox with context for the exemplar

    bbox = double(object.extent);
    [cx, cy, w, h] = deal(bbox(1)+bbox(3)/2, bbox(2)+bbox(4)/2, bbox(3), bbox(4));
    wc_z = w + context_amount*(w+h);
    hc_z = h + context_amount*(w+h);
    s_z = sqrt(single(wc_z*hc_z));
    scale_z = size_z / s_z;
    [im_crop_z, left_pad_z, top_pad_z, right_pad_z, bottom_pad_z] = get_subwindow_avg(im, [cy cx], [size_z size_z], [round(s_z) round(s_z)]);
    pad_z = ceil([scale_z*(left_pad_z+1) scale_z*(top_pad_z+1) size_z-scale_z*(right_pad_z+left_pad_z) size_z-scale_z*(top_pad_z+bottom_pad_z+1)]);
    %% Get instance sample
    d_search = (size_x - size_z)/2;
    pad = d_search/scale_z;
    s_x = s_z + 2*pad;
    scale_x = size_x / s_x;
    [im_crop_x, left_pad_x, top_pad_x, right_pad_x, bottom_pad_x] = get_subwindow_avg(im, [cy cx], [size_x size_x], [round(s_x) round(s_x)]);
    pad_x = ceil([scale_x*(left_pad_x+1) scale_x*(top_pad_x+1) size_x-scale_x*(right_pad_x+left_pad_x) size_x-scale_x*(top_pad_x+bottom_pad_x+1)]);
    % Size of object within the crops
    ws_z = w * scale_z;
    hs_z = h * scale_z;
    ws_x = w * scale_x;
    hs_x = h * scale_x;
    bbox_z = [(size_z-ws_z)/2, (size_z-hs_z)/2, ws_z, hs_z];
    bbox_x = [(size_x-ws_x)/2, (size_x-hs_x)/2, ws_x, hs_x];
end

% ---------------------------------------------------------------------------------------------------------------
function [im_patch, left_pad, top_pad, right_pad, bottom_pad] = get_subwindow_avg(im, pos, model_sz, original_sz)
%GET_SUBWINDOW_AVG Obtain image sub-window, padding with avg channel if area goes outside of border
% ---------------------------------------------------------------------------------------------------------------

    avg_chans = [mean(mean(im(:,:,1))) mean(mean(im(:,:,2))) mean(mean(im(:,:,3)))];

    if isempty(original_sz)
        original_sz = model_sz;
    end
    sz = original_sz;
    im_sz = size(im);
    %make sure the size is not too small
    assert(all(im_sz(1:2) > 2));
    c = (sz+1) / 2;

    %check out-of-bounds coordinates, and set them to avg_chans
    context_xmin = round(pos(2) - c(2));
    context_xmax = context_xmin + sz(2) - 1;
    context_ymin = round(pos(1) - c(1));
    context_ymax = context_ymin + sz(1) - 1;
    left_pad = double(max(0, 1-context_xmin));
    top_pad = double(max(0, 1-context_ymin));
    right_pad = double(max(0, context_xmax - im_sz(2)));
    bottom_pad = double(max(0, context_ymax - im_sz(1)));

    context_xmin = context_xmin + left_pad;
    context_xmax = context_xmax + left_pad;
    context_ymin = context_ymin + top_pad;
    context_ymax = context_ymax + top_pad;

    if top_pad || left_pad
        R = padarray(im(:,:,1), [top_pad left_pad], avg_chans(1), 'pre');
        G = padarray(im(:,:,2), [top_pad left_pad], avg_chans(2), 'pre');
        B = padarray(im(:,:,3), [top_pad left_pad], avg_chans(3), 'pre');
        im = cat(3, R, G, B);
    end

    if bottom_pad || right_pad
        R = padarray(im(:,:,1), [bottom_pad right_pad], avg_chans(1), 'post');
        G = padarray(im(:,:,2), [bottom_pad right_pad], avg_chans(2), 'post');
        B = padarray(im(:,:,3), [bottom_pad right_pad], avg_chans(3), 'post');
        im = cat(3, R, G, B);
    end

    xs = context_xmin : context_xmax;
    ys = context_ymin : context_ymax;

    im_patch_original = im(ys, xs, :);
    if ~isequal(model_sz, original_sz)
        im_patch = imresize(im_patch_original, model_sz);
    else
        im_patch = im_patch_original;
    end
end

author: bertinetto

repo: siamese-fc


text:
classdef SMD10 < PROBLEM
% <multi> <real> <constrained> <bilevel>
% Bilevel optimization problems proposed by Sinha, Malo, and Deb
% maxFElower --- 500 --- Maximum number of lower level function evaluations for each solution

%------------------------------- Reference --------------------------------
% A. Sinha, P. Malo, K. Deb, Test problem construction for single-objective
% bilevel optimization, Evolutionary Computation, 2014, 22(3): 439-477.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(SetAccess = private)
        maxFElower; % Maximum number of lower level function evaluations for each solution
        DU;         % Number of decision variables of the upper level
        DL;         % Number of decision variables of the lower level
        C;       	% Number of upper constraints
        p;          % The length of xu1
        q;          % The length of xl1
        r;          % The length of xu2 and xl2
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.maxFElower = obj.ParameterSet(500);
            obj.M = 2;
            if isempty(obj.D); obj.D = 5; end
            obj.DU = floor(obj.D/2);
            obj.DL = obj.D - obj.DU;
            obj.C  = obj.DU;
            obj.r  = floor(obj.DU/2);
            obj.p  = obj.DU - obj.r;
            obj.q  = obj.DL - obj.r;
            obj.lower    = [-5*ones(1,obj.DU), -5*ones(1,obj.q), -pi/2*ones(1,obj.r)+1e-6];
            obj.upper    = [10*ones(1,obj.DU), 10*ones(1,obj.q),  pi/2*ones(1,obj.r)-1e-6];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate upper level and lower level objective values
        function PopObj = CalObj(obj,PopDec)
            xu1 = PopDec(:,1:obj.p);
            xu2 = PopDec(:,obj.p+1:obj.p+obj.r);
            xl1 = PopDec(:,obj.p+obj.r+1:obj.p+obj.r+obj.q);
            xl2 = PopDec(:,obj.p+obj.r+obj.q+1:end);
            % Calculate  UL objective values
            PopObj(:,1) = sum((xu1-2).^2,2) + sum(xl1.^2,2) + sum((xu2-2).^2,2) - sum((xu2-tan(xl2)).^2,2);
            % Calculate  LL objective values
            PopObj(:,2) = sum(xu1.^2,2) + sum((xl1-2).^2,2) + sum((xu2-tan(xl2)).^2,2);
        end
        %% Calculate upper level and lower level constraint violations
        function PopCon = CalCon(obj,PopDec)
            xu1 = PopDec(:,1:obj.p);
            xu2 = PopDec(:,obj.p+1:obj.p+obj.r);
            xl1 = PopDec(:,obj.p+obj.r+1:obj.p+obj.r+obj.q);
            xl2 = PopDec(:,obj.p+obj.r+obj.q+1:end);
            % Upper level constraint violation
            PopCon(:,1:obj.p) = -xu1 - xu1.^3 + sum(xu1.^3,2) + sum(xu2.^3,2);
            PopCon(:,obj.p+1:obj.p+obj.r) = -xu2 - xu2.^3 + sum(xu2.^3,2) + sum(xu1.^3,2);
            % Lower level constraint violation
            PopCon(:,obj.p+obj.r+1:obj.p+obj.r+obj.q) = -xl1 - xl1.^3 + sum(xl1.^3,2);
        end
        %% Calculate lower level objective values
        function llPopulation = EvaluationLower(obj,varargin)
            PopDec            = obj.CalDec(varargin{1});
            PopObj            = obj.CalObj(PopDec);
            PopObj(:,1)       = nan;
            PopCon            = obj.CalCon(PopDec);
            PopCon(:,1:obj.C) = nan;
            llPopulation      = SOLUTION(PopDec,PopObj,PopCon,varargin{2:end});
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [edgeResponse, commonSize] = calculateEdgeCost(edgeImg, labels,labelIndices, labOne, labTwo)
edgeResponse = 0;
commonSize = 0;
% Find common border..

[rC,cC] = size(labels);
temp = labelIndices{1,labOne};
rowArr = temp(:,1);
colArr = temp(:,2);
[realSize,~] = size(rowArr);

totalResp  = 0;
commonSize = 0;
for i=1:realSize
    xCor = rowArr(i,1);
    yCor = colArr(i,1);

    if xCor + 1 <= rC
        if labels(xCor+1,yCor) == labTwo % Common border!
            commonSize = commonSize + 1;
            if edgeImg(xCor,yCor) ~= edgeImg(xCor+1,yCor)
                totalResp = totalResp + 1;
            end
        end
    end

    if xCor - 1 > 0
         if labels(xCor-1,yCor) == labTwo % Common border!
            commonSize = commonSize + 1;
            if edgeImg(xCor,yCor) ~= edgeImg(xCor-1,yCor)
                totalResp = totalResp + 1;
            end
        end
    end

    if yCor + 1 <= cC
         if labels(xCor,yCor+1) == labTwo % Common border!
            commonSize = commonSize + 1;
            if edgeImg(xCor,yCor) ~= edgeImg(xCor,yCor+1)
                totalResp = totalResp + 1;
            end
        end
    end

    if yCor - 1 > 0
        if labels(xCor,yCor-1) == labTwo % Common border!
            commonSize = commonSize + 1;
            if edgeImg(xCor,yCor) ~= edgeImg(xCor,yCor-1)
                totalResp = totalResp + 1;
            end
        end
    end
end

if commonSize ~=0
    edgeResponse = totalResp / commonSize;
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/calculateEdgeCost.m


text:
function output=fcn_LeaderAcceleration(u)
global leaderNum dim samplePos_all nodenum formationType simTime
global rotateFlag tempp2 tempp3

u_all = reshape(u(1:dim*nodenum*3),dim,nodenum*3);
v_all=u_all(:,1:nodenum);
p_all=u_all(:,nodenum+1:2*nodenum);
control_all=u_all(:,2*nodenum+1:3*nodenum);
currentTime=u(end);
if currentTime-floor(currentTime)==0
    disp(strcat('Completed %', num2str(currentTime/simTime*100,'%.1f')))
end
%output
a_leader_all=zeros(dim,leaderNum);

% speed up to speed=1
time_start_speedup=0;
time_span_speedup=15;
if currentTime>time_start_speedup && currentTime<time_start_speedup+time_span_speedup
    ax=fcn_sin(time_start_speedup,time_span_speedup,currentTime)/pi;
    a_leader_all(:,1)=[ax,0]';
    a_leader_all(:,2)=[ax,0]';
    a_leader_all(:,3)=[ax,0]';
end

% scale down to pass through obstacles
time_start_scaledown=time_start_speedup+time_span_speedup;
time_span_scaledown=10;
if currentTime>=time_start_scaledown && currentTime<=time_start_scaledown+time_span_scaledown
    % agent 3
    a_leader_all(2,3)=-fcn_SlowDownSpeedUp(time_start_scaledown, time_span_scaledown, currentTime, 1);
    % agent 3
    a_leader_all(2,2)=fcn_SlowDownSpeedUp(time_start_scaledown, time_span_scaledown, currentTime, 1);
end
% scale up after passing through obstacles
time_start_scaleup=time_start_scaledown+time_span_scaledown+5;
time_span_scaleup=10;
if currentTime>=time_start_scaleup && currentTime<=time_start_scaleup+time_span_scaleup
    % agent 3
    a_leader_all(2,3)=fcn_SlowDownSpeedUp(time_start_scaleup, time_span_scaleup, currentTime, 1);
    % agent 3
    a_leader_all(2,2)=-fcn_SlowDownSpeedUp(time_start_scaleup, time_span_scaleup, currentTime, 1);
end

% slow down to speed=0
time_start_slowdown=time_start_scaleup+time_span_scaleup+5;
time_span_slowdown=5;
if currentTime>time_start_slowdown && currentTime<time_start_slowdown+time_span_slowdown
    %ax=fcn_normalDistribution(time_start_slowdown,time_span_slowdown, 2, currentTime);
    ax=-fcn_sin(time_start_slowdown,time_span_slowdown,currentTime)/pi;
    a_leader_all(:,1)=[ax,0]';
    a_leader_all(:,2)=[ax,0]';
    a_leader_all(:,3)=[ax,0]';
end

% rotate
time_start_rotate=time_start_slowdown+time_span_slowdown;
time_span_rotate=10;
if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate
    % agent 2
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate, time_span_rotate, currentTime, 4);
    acce_vec=[1 0]';
    a_leader_all(:,2)=acce_mag*acce_vec;
    % agent 3
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate, time_span_rotate, currentTime, 4);
    acce_vec=[0 1]';
    a_leader_all(:,3)=acce_mag*acce_vec;
end

% move downward
time_start_downward=time_start_rotate;
time_span_downward=time_span_rotate;
if currentTime>=time_start_downward && currentTime<=time_start_downward+time_span_downward
    % agent 1
    ay=fcn_sin(time_start_downward,time_span_downward,currentTime)/pi;
    a_leader_all(:,1)=[0,-ay]';
    % agent 2
    a_leader_all(:,2)=a_leader_all(:,2)+[0,-ay]';
    % agent 3
    a_leader_all(:,3)=a_leader_all(:,3)+[0,-ay]';
end

% % one leader move to avoid obstacle
% time_start_agent2rightward=time_start_downward+time_span_rotate+5;
% time_span_agent2rightward=10;
% if currentTime>=time_start_agent2rightward && currentTime<=time_start_agent2rightward+time_span_agent2rightward
%     % agent 2
%     a_leader_all(1,2)=-fcn_SlowDownSpeedUp(time_start_agent2rightward, time_span_agent2rightward, currentTime, 4);
%     a_leader_all(1,3)=fcn_SlowDownSpeedUp(time_start_agent2rightward, time_span_agent2rightward, currentTime, 4);
% end
% %
% time_start_agent2leftward=time_start_agent2rightward+time_span_agent2rightward+5;
% time_span_agent2leftward=10;
% if currentTime>=time_start_agent2leftward && currentTime<=time_start_agent2leftward+time_span_agent2leftward
%     % agent 2
%     a_leader_all(1,2)=fcn_SlowDownSpeedUp(time_start_agent2leftward, time_span_agent2leftward, currentTime, 4);
%     a_leader_all(1,3)=-fcn_SlowDownSpeedUp(time_start_agent2leftward, time_span_agent2leftward, currentTime, 4);
% end

% slow down to speed=0
time_start_slowdown2=time_start_rotate+time_span_downward+10;
time_span_slowdown2=5;
if currentTime>time_start_slowdown2 && currentTime<time_start_slowdown2+time_span_slowdown2
    %ax=fcn_normalDistribution(time_start_slowdown,time_span_slowdown, 2, currentTime);
    ay=-fcn_sin(time_start_slowdown2,time_span_slowdown2,currentTime)/pi;
    a_leader_all(:,1)=[0,-ay]';
    a_leader_all(:,2)=[0,-ay]';
    a_leader_all(:,3)=[0,-ay]';
end

% rotate
time_start_rotate2=time_start_slowdown2+time_span_slowdown2;
time_span_rotate2=10;
if currentTime>=time_start_rotate2 && currentTime<=time_start_rotate2+time_span_rotate2
    % agent 2
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate2, time_span_rotate2, currentTime, 4);
    acce_vec=[0 -1]';
    a_leader_all(:,2)=acce_mag*acce_vec;
    % agent 3
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate2, time_span_rotate2, currentTime, 4);
    acce_vec=[1 0]';
    a_leader_all(:,3)=acce_mag*acce_vec;
end

% speed up to move leftward
time_start_leftward=time_start_rotate2;
time_span_leftward=time_span_rotate2;
if currentTime>=time_start_leftward && currentTime<=time_start_leftward+time_span_leftward
    % agent 1
    ax=-fcn_sin(time_start_leftward,time_span_leftward,currentTime)/pi;
    a_leader_all(:,1)=[ax,0]';
    % agent 2
    a_leader_all(:,2)=a_leader_all(:,2)+[ax,0]';
    % agent 3
    a_leader_all(:,3)=a_leader_all(:,3)+[ax,0]';
end

% affine transformation: agent 3 move forward+scale down in the vertical direction
time_start_affineVertical=time_start_leftward+time_span_leftward+5;
time_span_affineVertical=20;
if currentTime>=time_start_affineVertical && currentTime<=time_start_affineVertical+time_span_affineVertical
    % agent 3
    a_leader_all(1,3)=-fcn_SlowDownSpeedUp(time_start_affineVertical, time_span_affineVertical, currentTime, 1.5);
    %
    a_leader_all(2,3)=fcn_SlowDownSpeedUp(time_start_affineVertical, time_span_affineVertical, currentTime, 2);
    % agent 3
    a_leader_all(2,2)=-fcn_SlowDownSpeedUp(time_start_affineVertical, time_span_affineVertical, currentTime, 2);
end

% affine transformation: back to normal formation
time_start_backtonormal=time_start_affineVertical+time_span_affineVertical+10;
time_span_backtonormal=5;
if currentTime>=time_start_backtonormal && currentTime<=time_start_backtonormal+time_span_backtonormal
    % x-direction: agent 3
    a_leader_all(1,3)=fcn_SlowDownSpeedUp(time_start_backtonormal, time_span_backtonormal, currentTime, 1.5);
    % y-direction
    % agent 3
    a_leader_all(2,3)=-fcn_SlowDownSpeedUp(time_start_backtonormal, time_span_backtonormal, currentTime, 2);
    % agent 2
    a_leader_all(2,2)=fcn_SlowDownSpeedUp(time_start_backtonormal, time_span_backtonormal, currentTime, 2);
end

% output: replace the leaders' velocity with the one prescribed above
control_all(:,1:leaderNum) = a_leader_all;
output=reshape(control_all, dim*nodenum, 1);

repo: paper-simulation


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Created: May 27th, 2015
% Institution Created: UNC-CH
% Date Modified: July 13, 2021
% Institution Modified: TCNJ
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in (x,y) positions of the immersed boundary from .vtk
%           format
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = read_Lagrangian_Data_From_vtk(path,simNums)

cd(path);

filename = ['lagsPts.' num2str(simNums) '.vtk'];  % desired lagsPts.xxxx.vtk file

fileID = fopen(filename);
if ( fileID== -1 )
    error('\nCANNOT OPEN THE FILE!');
end

str = fgets(fileID); %-1 if eof
if ~strcmp( str(3:5),'vtk');
    error('\nNot in proper VTK format');
end

% read in the header info %
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID); % Up to white space in previous lagsPts.X files (pre June 2021)

% Check whether VTK file has time info. This is a VTK file with time,
%       need to read the next 3 lines to have read in appropriate
%       number of header lines.
if ~strcmp( str(1), 'P')
	str = fgets(fileID);
    str = fgets(fileID);
    str = fgets(fileID);
end

% stores # of Lagrangian Pts. as stated in .vtk file
numLagPts = sscanf(str,'%*s %f %*s',1);

% read in the vertices %
[mat,count] = fscanf(fileID,'%f %f %f',3*numLagPts);
if count ~= 3*numLagPts
   error('\nProblem reading in Lagrangian Pts.');
end

mat = reshape(mat, 3, count/3); % Reshape vector -> matrix (every 3 entries in vector make into matrix row)
vertices = mat';              % Store vertices in new matrix

fclose(fileID);               % Closes the data file.

xLag = vertices(:,1);         % x-Lagrangian Pts.
yLag = vertices(:,2);         % y-Lagrangian Pts.

cd ..;                        % Change directory back to ../viz_IB2d/ directory

clear mat str filename fileID count;

author: nickabattista


text:
classdef GeometricRefFrameSet < matlab.mixin.SetGet
    %GeometricRefFrameSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        refFrames(1,:) AbstractGeometricRefFrame = AbstractGeometricRefFrame.empty(1,0);

        lvdData LvdData
    end

    methods
        function obj = GeometricRefFrameSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addRefFrame(obj, refFrame)
            obj.refFrames(end+1) = refFrame;
        end

        function removeRefFrame(obj, refFrame)
            obj.refFrames([obj.refFrames] == refFrame) = [];
        end

        function [listBoxStr, refFrames] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.refFrames))
                listBoxStr{end+1} = obj.refFrames(i).getListboxStr(); %#ok<AGROW>
            end

            refFrames = obj.refFrames;
        end

        function indCoordSyses = getRefFramesForInds(obj, inds)
            indCoordSyses = obj.refFrames(inds);
        end

        function inds = getIndsForRefFrames(obj, indRefFrames)
            inds = find(ismember(obj.refFrames, indRefFrames));
        end

        function indRefFrame = getRefFrameAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.refFrames))
                indRefFrame = obj.refFrames(ind);
            else
                indRefFrame = AbstractGeometricRefFrame.empty(1,0);
            end
        end

        function numRefFrames = getNumRefFrames(obj)
            numRefFrames = length(obj.refFrames);
        end

        function tf = usesGroundObj(obj, groundObj)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGroundObj(groundObj);
            end
        end

        function tf = usesGeometricPoint(obj, point)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricPoint(point);
            end
        end

        function tf = usesGeometricVector(obj, vector)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricVector(vector);
            end
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricCoordSys(coordSys);
            end
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricRefFrame(refFrame);
            end
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricAngle(angle);
            end
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = false;
            for(i=1:length(obj.refFrames))
                tf = tf || obj.refFrames(i).usesGeometricPlane(plane);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function activate_ap(speed, host, port)
% ACTIVATE_AP(SPEED, HOST, PORT)
%
% Activate Autopilot in Flightgear over socket connection (HOST,
% PORT) and set target speed = SPEED [mph].
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);

% activate autopilot
out.println(sprintf(...
    'set /autopilot/locks/altitude altitude-hold\r\n'));
out.println(sprintf(...
    'set /autopilot/locks/heading true-heading-hold\r\n'));
out.println(sprintf(...
    'set /autopilot/locks/speed speed-with-throttle\r\n'));

% set target speed (convert mph to knots)
out.println(sprintf(...
    'set /autopilot/settings/target-speed-kt %d\r\n',speed*0.869));

out.close();
fg_socket.close();

end


text:
function progressbar(n,N,w)

% progressbar - display a progress bar
%
%    progressbar(n,N,w);
%
% displays the progress of n out of N.
% n should start at 1.
% w is the width of the bar (default w=20).
%
%   Copyright (c) Gabriel Peyre 2006

if nargin<3
    w = 20;
end

% progress char
cprog = '.';
cprog1 = '*';
% begining char
cbeg = '[';
% ending char
cend = ']';

p = min( floor(n/N*(w+1)), w);

global pprev;
if isempty(pprev)
    pprev = -1;
end

if not(p==pprev)
    ps = repmat(cprog, [1 w]);
    ps(1:p) = cprog1;
    ps = [cbeg ps cend];
    if n>1
        % clear previous string
        fprintf( repmat('\b', [1 length(ps)]) );
    end
    fprintf(ps);
end
pprev = p;
if n==N
    fprintf('\n');
end

author: gpeyre

repo: matlab-toolboxes


text:
function best = findWeakRuleSamplesFeatures2(data,labels,dist,binVals,bins)

numDim = size(data,2);

best.dim = 1;
best.error = 0.5;
best.dir = 1;
best.tr = 0;

curBestErr = 0.5*ones(1,numDim);
binNo = ones(1,numDim);
bestDir = ones(1,numDim);

numS = 2500;
curSel = zeros(1,numS);
cumD = cumsum(dist);
for ndx = 1:numS
  curSel(ndx) = sum(cumD<rand(1,1))+1;
end

curBins = bins(curSel,:);
curLabels = labels(curSel);
curSelDim = find(rand(1,numDim)>0.9);

parfor dimNdx = 1:numel(curSelDim)
  dim = curSelDim(dimNdx);
  numBins = size(binVals,1)+1;
  binNdx =  curBins(:,dim)+numBins*(curLabels>0);
  allCount = histc(binNdx,.5:(2*numBins+0.5));
  allCount = allCount/sum(allCount);
  posCount = allCount(numBins+1:end-1);
  negCount = allCount(1:numBins);

  posLeft = 0;
  posRight = sum(posCount);
  negLeft = 0;
  negRight = sum(negCount);

  err = zeros(1,size(binVals,1));
  dir = ones(1,size(binVals,1));
  for ndx = 1:size(binVals,1)
    posLeft = posLeft + posCount(ndx);
    posRight = posRight - posCount(ndx);
    negLeft = negLeft + negCount(ndx);
    negRight = negRight - negCount(ndx);

    err(ndx) = posRight+negLeft-posLeft-negRight;
    if(err(ndx)<0)
      err(ndx) = - err(ndx);
      dir(ndx) = -1;
    end
  end
  err = 0.5-err/2;
  [curBestErr(dimNdx) binNo(dimNdx)]= min(err);
  bestDir(dimNdx) = dir(binNo(dimNdx));
end

[minError minDimNdx] = min(curBestErr);
minDim = curSelDim(minDimNdx);
best.error = minError;   best.dim = minDim;
best.dir = bestDir(minDimNdx);   best.tr = binVals(binNo(minDim),minDim);

end

%{
function [predError dir] = getError(data,label,tr, dist)

dir = 1;
predLabel = 2*(data>tr)-1;
predError = sum(( (predLabel.*label) ~=1).*dist);

if(predError>0.5); predError = 1-predError; dir = -1; end
end
%}

author: kristinbranson

repo: JAABA


text:
%DEMO_WAVELETS  Wavelet filter banks
%
%   This demo exemplifies the use of the wavelet filter bank trees. All
%   representations use "least asymmetric" Daubechies wavelet orthonormal
%   filters 'sym8' (8-regular, length 16).
%
%   .. figure::
%
%      DWT representation
%
%      The filter bank tree consists of 11 levels of iterated 2-band basic
%      wavelet filter bank, where only the low-pass output is further
%      decomposed. This results in 12 bands with octave resolution.
%
%   .. figure::
%
%      8-band DWT representation
%
%      The filter bank tree (effectively) consists of 3 levels of iterated
%      8-band basic wavelet filter bank resulting in 22 bands. Only the
%      low-pass output is decomposed at each level.
%
%   .. figure::
%
%      Full Wavelet filter bank tree representation
%
%      The filter bank tree depth is 8 and it is fully decomposed meaning
%      both outputs (low-pass and high-pass) of the basic filter bank is
%      plot further. This results in 256 bands linearly covering the
%      frequency axis.
%
%   .. figure::
%
%      Full Wavelet filter bank tree representation
%
%      The same case as before, but symmetric nearly orthogonal basic
%      filter bank is used.
%
%   .. figure::
%
%      Full Dual-tree Wavelet filter bank representation
%
%      This is a 2 times redundant representation using Q-shift dual-tree
%      wavelet filters.
%

% Read the test signal and crop it to the range of interest
[f,fs]=gspi;
f=f(10001:100000);
dr=50;

figure(1);
[c1,info]=fwt(f,'sym8',11);
plotwavelets(c1,info,fs,'dynrange',dr);

figure(2);
[c2,info]=wfbt(f,{'sym8',3,'quadband'});
plotwavelets(c2,info,fs,'dynrange',dr);

figure(3);
[c3,info]=wfbt(f,{'sym8',8,'full'});
plotwavelets(c3,info,fs,'dynrange',dr);

figure(4);
[c4,info]=wfbt(f,{'symorth3',8,'full'});
plotwavelets(c4,info,fs,'dynrange',dr);

figure(5);
[c5,info]=dtwfbreal(f,{'qshift5',8,'full','first','symorth3'});
plotwavelets(c5,info,fs,'dynrange',dr);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat