text:
function [G,E,R,data] = arap_gradient(V,F,U,varargin)
  % ARAP_GRADIENT Compute the Gradient of an as-rigid-as-possible for a mesh in
  % rest configuration (V,F) with vertices now at U, according to "A Simple
  % Geometric Model for Elastic Deformations" [Chao et al. 2010]
  %
  % G = arap_gradient(V,F,U)
  % [G,E,R] = arap_gradient(V,F,U,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   V  #V by dim rest vertex positions
  %   F  #F by simplex-size simplex indices into V
  %   U  #V by dim deformed vertex positions
  %     Optional:
  %       'Rotations' followed by R  dim by dim by #R list of best fit
  %         rotations
  %       'Energy' followed by either 'spokes','spokes-and-rims','elements'
  % Outputs:
  %   G  #V by dim list of gradient vectors per vertex
  %   E  arap energy
  %   R  dim by dim by #R list of best fit rotations
  %   data
  %
  % See also: arap, arap_hessian
  %
  % Example:
  % % Given a mesh (V,F) and deformed positions U0, flow to energy minimum
  % % using Newton's method.
  % clf;
  % hold on;
  %   tsurf(F,V,fphong,'FaceColor','r','SpecularStrength',0,'AmbientStrength',0.5);
  %   t = tsurf(F,U0,fphong,'FaceColor','b','SpecularStrength',0,'AmbientStrength',0.5);
  % hold off;
  % axis equal;
  % view(2);
  % camlight;
  % U = U0;
  % delta_t = 1e-1;
  % while true
  %   [G,E] = arap_gradient(V,F,U);
  %   U = U - delta_t * G;
  %   U = bsxfun(@plus,U,mean(V)-mean(U)+[max(V(:,1))-min(V(:,1)) 0 0]);
  %   t.Vertices = U;
  %   title(sprintf('E = %g\n',E),'FontSize',20);
  %   drawnow;
  % end
  %

  % default values
  switch size(F,2)
  case 4
    energy = 'elements';
  case 3
    energy = 'spokes-and-rims';
  end
  single_precision = true;
  R = [];
  data = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Energy','Rotations','Data','SinglePrecision'}, ...
    {'energy','R','data','single_precision'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % TODO: implement "flat" arap like `arap.m`, for now use placeholders with
  % non-flat defaults.
  ref_V = V;
  ref_F = F;
  dim = size(ref_V,2);
  flat = false;

  if isempty(data)
    ss = size(F,2);
    data.L = cotmatrix(V,F);
    if isempty(R)
      data.CSM = covariance_scatter_matrix(ref_V,ref_F,'Energy',energy);
    end
    [~,data.K] = arap_rhs(ref_V,ref_F,[],'Energy',energy);
  end

  % compute covariance matrix elements
  S = zeros(size(data.CSM,1),dim);
  S(:,1:dim) = data.CSM*repmat(U,dim,1);
  % dim by dim by n list of covariance matrices
  SS = permute(reshape(S,[size(data.CSM,1)/dim dim dim]),[2 3 1]);
  % fit rotations to each deformed vertex
  R = fit_rotations(SS,'SinglePrecision',single_precision);

  nr = size(R,3);
  Rcol = reshape(permute(R,[3 1 2]),nr*dim*dim,1);
  dV = data.K * Rcol;
  dV = reshape(dV,[size(V,1) dim]);

  G = -(data.L*U + dV);
  if nargout > 1
    E = trace(-U'*0.5*data.L*U - U'*dV - V'*0.5*data.L*V);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function z_callBack_virtualTeachPendant_J1_minus()
%% Soft realtime joints control Joint X button callback
% Copyright: Mohammad SAFEEA, 27th-July-2018
global virtual_smartPad_w_joints_command;
disp('Moving J1 in negative direction');
virtual_smartPad_w_joints_command(1)=-1;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function x = p34_start ( n )

%*****************************************************************************80
%
%% P34_START returns a starting point for optimization for problem 34.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ -1.5, 0.5 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function db_nc_audio_recon(file_name);
% Implements 4-channel tree-structured PR filter bank using 5th order
% Butterworth half-band filters constructed with 1st order allpass
% structures (ref. [Mit01] P10.40).  Uses double-buffer, non-causal IIR
% synthesis filtering technique (ref. [Cre96]) to process and play
% original multi-channel audio file, reconstructed audio file using uniform
% 16-bit quantization and PR structure described above, and error signal.

clc

B = 16; % Quantization bits
BU = 100; % Buffer length (even)
ST = 289; % Signal truncation parameter

% Process input data
[x1,fs] = wavread(file_name);
x1 = x1';
sz = size(x1);

if sz(2) > 8*BU*ST
    for n = 1:sz(1)
        x(n,:) = x1(n,1:8*BU*ST); % Truncate signal for simplicity, odd multiple of BU
    end
else
    x = x1;
end

x = [x zeros(sz(1),8*BU - rem(sz(2),8*BU))]; % Make input even length for all branches

% Perform four-channel, tree-structured filtering of multi-channel input
for m = 1:sz(1)

    % Analysis stage
    [H00,L00,Z1f00,Z0f00] = nciab(x(m,:),BU);
    [H01,L01,Z1f01,Z0f01] = nciab(L00,BU);
    [H10,L10,Z1f10,Z0f10] = nciab(L01,BU);

    % Perform quantization
    Q00 = uquan(H00,B);
    Q01 = uquan(H01,B);
    Q10 = uquan(H10,B);
    Q11 = uquan(L10,B);

    % Synthesis stage
    Y00 = ncisb(H10,L10,BU,Z1f10,Z0f10);
    Y01 = ncisb(H01,Y00,BU,Z1f01,Z0f01);
    y1(m,:) = ncisb(H00,Y01,BU,Z1f00,Z0f00);
end

% Produce error signal
szx = size(x);
szy = size(y1);
for n = 1:szy(1)
    y(n,:) = y1(n,1:szx(2)); % Truncate output signal to input signal length
end
z = x - y;

% Produce audio output
disp('Original audio file (PRESS ANY KEY)');
pause
wavplay(x',fs)
disp('Processed audio file (PRESS ANY KEY)');
pause
wavplay(y',fs)
disp('Error signal (PRESS ANY KEY)');
pause
wavplay(z',fs)

% Plot error signal
for t = 1:sz(1)
    figure(t)
    plot(z(t,:))
end

author: Sable

repo: mcbench-benchmarks


text:
function [thisFrame,movieFileID,inputMovieDims] = readFrame(inputMoviePath,frameNo,varargin)
	% Fast reading of frame from files on disk. This is an alternative to loadMovieList that is much faster when only a single frame needs to be read.
	% Biafra Ahanonu
	% started: 2020.10.19 [11:45:24]
	% inputs
		% inputMoviePath | Str: path to a movie. Supports TIFF, AVI, HDF5, NWB, or Inscopix ISXD.
		% frameNo | Int: frame number.
	% outputs
		%
	% Usage
		% For non-HDF5 file types, need to open a link to the
			% [thisFrame,movieFileID,inputMovieDims] = ciapkg.io.readFrame(inputMoviePath,frameNo);
			% Then for the second call, feed in movieFileID and inputMovieDims to improve read speed.
			% [thisFrame] = ciapkg.io.readFrame(inputMoviePath,frameNo,'movieFileID',movieFileID,'inputMovieDims',inputMovieDims);

	% changelog
		% 2021.06.30 [01:26:29] - Updated handling of no file path character input.
		% 2021.07.03 [09:02:14] - Updated to have backup read method for different tiff styles.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	[thisFrame,movieFileID,inputMovieDims] = ciapkg.io.readFrame(inputMoviePath,frameNo,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%ANGVEC2TR Convert angle and vector orientation to a homogeneous transform
%
% T = ANGVEC2TR(THETA, V) is a homogeneous transform matrix equivalent to a
% rotation of THETA about the vector V.
%
% Note::
% - The translational part is zero.
%
% See also EUL2TR, RPY2TR, ANGVEC2R.

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

function T = angvec2tr(theta, k)

    if nargin < 2
        error('RTB:angvec2tr:badarg', 'bad arguments');
    end

    T = r2t( angvec2r(theta, k) );

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [x,zo,xs]=estnoisem_noiseseg(yf,noise_seg,pv01,tz,pp)
%ESTNOISEM - estimate noise spectrum using minimum statistics
%
% Usage:    ninc=round(0.016*fs);   % frame increment [fs=sample frequency]
%           ovf=2;                  % overlap factor
%           f=rfft(enframe(s,hanning(ovf*ninc,'periodic'),ninc),ovf*ninc,2);
%           f=f.*conj(f);           % convert to power spectrum
%           x=estnoisem(f,ninc/fs); % estimate the noise power spectrum
%
% Inputs:
%   yf      input power spectra (one row per frame)
%   tz      frame increment in seconds
%           Alternatively, the input state from a previous call (see below)
%   pp      algorithm parameters [optional]
%
% Outputs:
%   x       estimated noise power spectra (one row per frame)
%   zo      output state
%   xs      estimated std error of x (one row per frame)
%           xs seems often to be an underestimate by a factor of 2 or 3
%
% The algorithm parameters are defined in reference [1] from which equation
% numbers are given in parentheses. They are as follows:
%
%        pp.taca      % (11): smoothing time constant for alpha_c [0.0449 seconds]
%        pp.tamax     % (3): max smoothing time constant [0.392 seconds]
%        pp.taminh    % (3): min smoothing time constant (upper limit) [0.0133 seconds]
%        pp.tpfall    % (12): time constant for P to fall [0.064 seconds]
%        pp.tbmax     % (20): max smoothing time constant [0.0717 seconds]
%        pp.qeqmin    % (23): minimum value of Qeq [2]
%        pp.qeqmax    % max value of Qeq per frame [14]
%        pp.av        % (23)+13 lines: fudge factor for bc calculation  [2.12]
%        pp.td        % time to take minimum over [1.536 seconds]
%        pp.nu        % number of subwindows to use [3]
%        pp.qith      % Q-inverse thresholds to select maximum noise slope [0.03 0.05 0.06 Inf ]
%        pp.nsmdb     % corresponding noise slope thresholds in dB/second   [47 31.4 15.7 4.1]
%
% Example use:      y=enframe(s,w,ni);                  % divide speech signal s(n) into
%                                                       % overlapping frames using window w(n)
%                   yf=rfft(y,nf,2);                    % take fourier transform
%                   dp=estnoisem(yf.*conj(yf),tinc);    % estimate the noise
%
% If convenient, you can call estnoisem in chunks of arbitrary size. Thus the following are equivalent:
%
%                   (a) dp=estnoisem(yp(1:300),tinc);
%
%                   (b) [dp(1:100),z]=estnoisem(yp(1:100),tinc);
%                       [dp(101:200),z]=estnoisem(yp(101:200),z);
%                       [dp(201:300),z]=estnoisem(yp(201:300),z);

% This is intended to be a precise implementation of [1] with Table III
% replaced by the updated table 5 from [2]. The only deliberate algorithm
% change is the introduction of a minimum value for 1/Qeq in equation (23).
% This change only affects the first few frames and improves the
% convergence of the algorithm. A minor improveemnt was reported in [3] but
% this has not yet been included.
%
% Refs:
%    [1] Rainer Martin.
%        Noise power spectral density estimation based on optimal smoothing and minimum statistics.
%        IEEE Trans. Speech and Audio Processing, 9(5):504-512, July 2001.
%    [2] Rainer Martin.
%        Bias compensation methods for minimum statistics noise power spectral density estimation
%        Signal Processing, 2006, 86, 1215-1229
%    [3] Dirk Mauler and Rainer Martin
%        Noise power spectral density estimation on highly correlated data
%        Proc IWAENC, 2006

%	   Copyright (C) Mike Brookes 2008
%      Version: $Id: estnoisem.m 1718 2012-03-31 16:40:41Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Modified code, Zheng-Hua Tan, 2012
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nr,nrf]=size(yf);          % number of frames and freq bins

x=zeros(nr,nrf);            % initialize output arrays
xs=zeros(nr,nrf);           % will hold std error in the future
if isempty(yf) && isstruct(tz)             % no real data
    zo=tz;              % just keep the same state
else
    if isstruct(tz)       % take parameters from a previous call
        nrcum=tz.nrcum;
        p=tz.p;          % smoothed power spectrum
        ac=tz.ac;               % correction factor (9)
        sn2=tz.sn2;              % estimated noise power
        pb=tz.pb;               % smoothed noisy speech power (20)
        pb2=tz.pb2;
        pminu=tz.pminu;
        actmin=tz.actmin;   % Running minimum estimate
        actminsub=tz.actminsub;           % sub-window minimum estimate
        subwc=tz.subwc;                   % force a buffer switch on first loop
        actbuf=tz.actbuf;  % buffer to store subwindow minima
        ibuf=tz.ibuf;
        lminflag=tz.lminflag;      % flag to remember local minimum
        tinc=tz.tinc;     % frame increment
        qq=tz.qq;         % parameter structure
    else
        tinc = tz;          % second argument is frame increment
        nrcum=0;            % no frames so far
        % default algorithm constants

        qq.taca=0.0449;    % smoothing time constant for alpha_c = -tinc/log(0.7) in equ (11)
        qq.tamax=0.392;    % max smoothing time constant in (3) = -tinc/log(0.96)
        qq.taminh=0.0133;    % min smoothing time constant (upper limit) in (3) = -tinc/log(0.3)
        qq.tpfall=0.064;   % time constant for P to fall (12)
        qq.tbmax=0.0717;   % max smoothing time constant in (20) = -tinc/log(0.8)
        qq.qeqmin=2;       % minimum value of Qeq (23)
        qq.qeqmax=14;      % max value of Qeq per frame
        qq.av=2.12;             % fudge factor for bc calculation (23 + 13 lines)
        qq.td=1.536;       % time to take minimum over
        qq.nu=8;           % number of subwindows
        qq.qith=[0.03 0.05 0.06 Inf]; % noise slope thresholds in dB/s
        qq.nsmdb=[47 31.4 15.7 4.1];

        if nargin>=4 && ~isempty(pp)
            qqn=fieldnames(qq);
            for i=1:length(qqn)
                if isfield(pp,qqn{i})
                    qq.(qqn{i})=pp.(qqn{i});
                end
            end
        end
    end

    % unpack parameter structure

    taca=qq.taca;    % smoothing time constant for alpha_c = -tinc/log(0.7) in equ (11)
    tamax=qq.tamax;    % max smoothing time constant in (3) = -tinc/log(0.96)
    taminh=qq.taminh;    % min smoothing time constant (upper limit) in (3) = -tinc/log(0.3)
    tpfall=qq.tpfall;   % time constant for P to fall (12)
    tbmax=qq.tbmax;   % max smoothing time constant in (20) = -tinc/log(0.8)
    qeqmin=qq.qeqmin;       % minimum value of Qeq (23)
    qeqmax=qq.qeqmax;      % max value of Qeq per frame
    av=qq.av;             % fudge factor for bc calculation (23 + 13 lines)
    td=qq.td;       % time to take minimum over
    nu=qq.nu;           % number of subwindows
    qith=qq.qith; % noise slope thresholds in dB/s
    nsmdb=qq.nsmdb;   % maximum permitted +ve noise slope in dB/s

    % derived algorithm constants

    aca=exp(-tinc/taca); % smoothing constant for alpha_c in equ (11) = 0.7
    acmax=aca;          % min value of alpha_c = 0.7 in equ (11) also = 0.7
    amax=exp(-tinc/tamax); % max smoothing constant in (3) = 0.96
    aminh=exp(-tinc/taminh); % min smoothing constant (upper limit) in (3) = 0.3
    bmax=exp(-tinc/tbmax); % max smoothing constant in (20) = 0.8
    snrexp = -tinc/tpfall;
    nv=round(td/(tinc*nu));    % length of each subwindow in frames
    if nv<4            % algorithm doesn't work for miniscule frames
        nv=4;
        nu=max(round(td/(tinc*nv)),1);
    end
    nd=nu*nv;           % length of total window in frames
    [md,hd]=mhvals(nd); % calculate the constants M(D) and H(D) from Table III
    [mv,hv]=mhvals(nv); % calculate the constants M(D) and H(D) from Table III
    nsms=10.^(nsmdb*nv*tinc/10);  % [8 4 2 1.2] in paper
    qeqimax=1/qeqmin;  % maximum value of Qeq inverse (23)
    qeqimin=1/qeqmax; % minumum value of Qeq per frame inverse

    if isempty(yf)      % provide dummy initialization
        ac=1;               % correction factor (9)
        subwc=nv;                   % force a buffer switch on first loop
        ibuf=0;
        p=x;          % smoothed power spectrum
        sn2=p;              % estimated noise power
        pb=p;               % smoothed noisy speech power (20)
        pb2=pb.^2;
        pminu=p;
        actmin=repmat(Inf,1,nrf);   % Running minimum estimate
        actminsub=actmin;           % sub-window minimum estimate
        actbuf=repmat(Inf,nu,nrf);  % buffer to store subwindow minima
        lminflag=zeros(1,nrf);      % flag to remember local minimum
    else

        if ~nrcum

            %Initialization excluding 3 left, current (noisy frame), 3 right frames
            By=[]; Cx=[]; dC=[]; Ax=[];
            xNc=find(noise_seg==1); % current noise frame indexs
            L_3f=[ [xNc-3]' [xNc-2]' [xNc-1]' ]; %3 frames(indexes) 'left' of 'xNc'
            R_3f=[ [xNc+1]' [xNc+2]' [xNc+3]' ]; % 3 frames(indexes)  'right' of 'xNc'

            dC=[ L_3f xNc' R_3f ];   SxY=sort(unique(dC(:))); % (frame indexes which should not be considered)

            By=SxY(find(SxY>0)); % (discard 0 or negative frame indexes)
            Ax=1:1:length(noise_seg); Cx = setdiff(Ax,By); % 'By' all unsatisfied frames (ascending sorted order)

            if isempty(Cx) == 0 % satisfied frame indexes
               p=yf(Cx(1),:);   %first satisfied frame from begining (as 'Cx' asced. sorted order)

            else    % otherwise

                Cx=find(noise_seg==0); %speech frames
                if length (Cx) > 3
                   Cx(1:3)=[];  % avoid first 3 speech frames - insense of "only slight better" 'init of p'
                end
                By=[];       % removed 3L-current-3R indexes (else impose to discard all frames in file)- during 'p' update [ ~isempty(intersect(t,By)) ]

                if isempty(Cx) == 1 % not exist 4th speech frame/only contents noisy segment

                      Cx=1;   p=yf(1,:); %no hope (as usual)
                      warning('File seems have only noisy segments ..  noise-power-spectra-estimation may not good!');
                    else
                      p=yf(Cx(1),:); %i.e. 4th speech frame
                end

            end %

            ne_min=p;
            ac=1;               % correction factor (9)
            sn2=p;              % estimated noise power
            pb=p;               % smoothed noisy speech power (20)
            pb2=pb.^2;
            pminu=p;
            actmin=repmat(Inf,1,nrf);   % Running minimum estimate
            actminsub=actmin;           % sub-window minimum estimate
            subwc=nv;                   % force a buffer switch on first loop
            actbuf=repmat(Inf,nu,nrf);  % buffer to store subwindow minima
            ibuf=0;
            lminflag=zeros(1,nrf);      % flag to remember local minimum
        end

        % loop for each frame
        if sum(pv01(1:10))>=1
            p=min(yf(1:min(50,nr),:),[],1);
            ne_min=p; sn2=p; pb=p; pb2=pb.^2; pminu=p;
        end

        for t=Cx(1):nr % start from the "initialization frame for 'p'  %% we use t instead of lambda in the paper
            yft=yf(t,:);        % noise speech power spectrum
            acb=(1+(sum(p)./sum(yft)-1).^2).^(-1);  % alpha_c-bar(t)  (9)
            ac=aca*ac+(1-aca)*max(acb,acmax);       % alpha_c(t)  (10)
            ah=amax*ac.*(1+(p./sn2-1).^2).^(-1);    % alpha_hat: smoothing factor per frequency (11)
            snr=sum(p)/sum(sn2);
            ah=max(ah,min(aminh,snr^snrexp));       % lower limit for alpha_hat (12)

            if noise_seg(t) ||  ~isempty(intersect(t,By)) || (t<11 && sum(pv01(1:10))>=1) %avoid the frames in "By" (3Left, curr. noise,3Right)
            else
                p=ah.*p+(1-ah).*yft;            % smoothed noisy speech power (3)
                b=min(ah.^2,bmax);              % smoothing constant for estimating periodogram variance (22 + 2 lines)
                pb=b.*pb + (1-b).*p;            % smoothed periodogram (20)
                pb2=b.*pb2 + (1-b).*p.^2;     	% smoothed periodogram squared (21)
            end

            qeqi=max(min((pb2-pb.^2)./(2*sn2.^2),qeqimax),qeqimin/(t+nrcum));   % Qeq inverse (23)
            qiav=sum(qeqi)/nrf;             % Average over all frequencies (23+12 lines) (ignore non-duplication of DC and nyquist terms)
            bc=1+av*sqrt(qiav);             % bias correction factor (23+11 lines)
            bmind=1+2*(nd-1)*(1-md)./(qeqi.^(-1)-2*md);      % we use the simplified form (17) instead of (15)
            bminv=1+2*(nv-1)*(1-mv)./(qeqi.^(-1)-2*mv);      % same expression but for sub windows
            kmod=bc*p.*bmind<actmin;        % Frequency mask for new minimum

            if any(kmod)
                actmin(kmod)=bc*p(kmod).*bmind(kmod);
                actminsub(kmod)=bc*p(kmod).*bminv(kmod);
            end
            if subwc>1 && subwc<nv              % middle of buffer - allow a local minimum
                lminflag=lminflag | kmod;    	% potential local minimum frequency bins
                pminu=min(actminsub,pminu);
                sn2=pminu;
            else
                if subwc>=nv                    % end of buffer - do a buffer switch
                    ibuf=1+rem(ibuf,nu);     	% increment actbuf storage pointer
                    actbuf(ibuf,:)=actmin;    	% save sub-window minimum
                    pminu=min(actbuf,[],1);
                    i=find(qiav<qith);
                    nsm=nsms(i(1));          	% noise slope max
                    lmin=lminflag & ~kmod & actminsub<nsm*pminu & actminsub>pminu;
                    if any(lmin)
                        pminu(lmin)=actminsub(lmin);
                        actbuf(:,lmin)=repmat(pminu(lmin),nu,1);
                    end
                    lminflag(:)=0;
                    actmin(:)=Inf;
                    subwc=0;
                end
            end
            subwc=subwc+1;
            x(t,:)=sn2;
            %ne_min=min(ne_min, sn2);
            qisq=sqrt(qeqi);
            % empirical formula for standard error based on Fig 15 of [2]
            xs(t,:)=sn2.*sqrt(0.266*(nd+100*qisq).*qisq/(1+0.005*nd+6/nd)./(0.5*qeqi.^(-1)+nd-1));
        end
    end

    if nargout>1    % we need to store the state for next time
        zo.nrcum=nrcum+nr;      % number of frames so far
        zo.p=p;          % smoothed power spectrum
        zo.ac=ac;               % correction factor (9)
        zo.sn2=sn2;              % estimated noise power
        zo.pb=pb;               % smoothed noisy speech power (20)
        zo.pb2=pb2;
        zo.pminu=pminu;
        zo.actmin=actmin;   % Running minimum estimate
        zo.actminsub=actminsub;           % sub-window minimum estimate
        zo.subwc=subwc;                   % force a buffer switch on first loop
        zo.actbuf=actbuf;  % buffer to store subwindow minima
        zo.ibuf=ibuf;
        zo.lminflag=lminflag;      % flag to remember local minimum
        zo.tinc=tinc;     % must be the last one
        zo.qq=qq;
    end
    if ~nargout
        clf;
        subplot(212);
        plot((1:nr)*tinc,10*log10([sum(yf,2) sum(x,2)]))
        ylabel('Frame Energy (dB)');
        xlabel(sprintf('Time (s)   [%d ms frame incr]',round(tinc*1000)));
        axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
        subplot(211);
        plot(1:nrf,10*log10([sum(yf,1)'/nr sum(x,1)'/nr]))
        ylabel('Power (dB)');
        xlabel('Frequency bin');
        axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
    end
end

function [m,h,d]=mhvals(d)
% Values are taken from Table 5 in [2]
%[2] R. Martin,"Bias compensation methods for minimum statistics noise power
%               spectral density estimation", Signal Processing Vol 86, pp1215-1229, 2006.

% approx: plot(d.^(-0.5),[m 1-d.^(-0.5)],'x-'), plot(d.^0.5,h,'x-')
persistent dmh
if isempty(dmh)
    dmh=[
        1   0       0;
        2   0.26    0.15;
        5   0.48    0.48;
        8   0.58    0.78;
        10  0.61    0.98;
        15  0.668   1.55;
        20  0.705   2;
        30  0.762   2.3;
        40  0.8     2.52;
        60  0.841   3.1;
        80  0.865   3.38;
        120 0.89    4.15;
        140 0.9     4.35;
        160 0.91    4.25;
        180 0.92    3.9;
        220 0.93    4.1;
        260 0.935   4.7;
        300 0.94    5];
end

if nargin>=1
    i=find(d<=dmh(:,1));
    if isempty(i)
        i=size(dmh,1);
        j=i;
    else
        i=i(1);
        j=i-1;
    end
    if d==dmh(i,1)
        m=dmh(i,2);
        h=dmh(i,3);
    else
        qj=sqrt(dmh(i-1,1));    % interpolate using sqrt(d)
        qi=sqrt(dmh(i,1));
        q=sqrt(d);
        h=dmh(i,3)+(q-qi)*(dmh(j,3)-dmh(i,3))/(qj-qi);
        m=dmh(i,2)+(qi*qj/q-qj)*(dmh(j,2)-dmh(i,2))/(qi-qj);
    end
else
    d=dmh(:,1);
    m=dmh(:,2);
    h=dmh(:,3);
end

author: zhenghuatan

repo: rVAD

save_path: github-repos/MATLAB/zhenghuatan-rVAD


text:
function [Dimg,D,X,err,Z] = perform_dictionary_signature_learning(Y, k, options)

% perform_dictionary_signature_learning - compute an image signature
%
%   [Dimg,D,X,err,options.Z] = perform_dictionary_signature_learning(Y, k, w, options);
%
%   The algorithm is described in
%       M. Aharon and M. Elad,
%       " Sparse and Redundant Modeling of Image Content Using an Image-Signature-Dictionary",
%       Submitted.
%
%   Y is a set of exemplar, each Y(:,i)=p(:) is a w*w vector where p should
%       be some 2D patch of size (w,w) extracted from some image(s).
%   k is the width of the dictionary signature.
%
%   The number of iterations is options.niter.
%   The sparse coder used is set in options.sparse_coder to either 'omp' or
%       'mp'.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;

% width of the patches
ww = size(Y,1);
w = sqrt(ww);
% number of exemplar
m = size(Y,2);
d = k;
dd = d^2; % number of atoms in the dictionary
kk = k^2;

if m<2*dd
    warning('You should increase the number of exemplars.');
end

% option for OMP
if not(isfield(options, 'nbr_max_atoms'))
    options.nbr_max_atoms = 4;
end
if not(isfield(options, 'use_mex'))
    options.use_mex = 1;
end

Dimg = getoptions(options, 'Dimg', randn(k));
niter = getoptions(options, 'niter', 40);
use_cg = getoptions(options, 'use_cg', 1);
options.niter_max = getoptions(options, 'niter_cg', 10);
sparse_coder = getoptions(options, 'sparse_coder', 'omp');

Dimg = Dimg/sqrt( sum(Dimg(:).^2) );

if isfield(options, 'strict_sparsity')
    options.nbr_atoms_max = options.strict_sparsity;
end

% remove mean
mu = sum(Y);
Y = Y - repmat( mu/ww, [ww 1] );

% compute the dictionary extraction sparse matrix
if isfield(options, 'Z')
    Z = options.Z;
else
    Z = sparse(ww*dd,kk);
    num = 0;
    for x=1:d
        for y=1:d
            num = num+1;
            selx = x:x+w-1;
            sely = y:y+w-1;
            selx = mod(selx-1,k)+1;
            sely = mod(sely-1,k)+1;
            [Ya,Xa] = meshgrid(sely,selx);
            sel = Xa(:) + k*(Ya(:)-1);
            Z( (1:ww)' + (num-1)*ww + (sel-1)*dd*ww ) = 1;
        end
    end
end
err = []; Dimg = Dimg(:);
for i=1:niter
    disp(['--> Iteration ' num2str(i) '/' num2str(niter) '.']);
    % extract dictionary and precompute the matrices
    D = reshape( Z*Dimg, ww, dd);
    % remove zero component and normalize
    D = D - repmat( mean(D), [ww 1] );
    D = D ./ repmat( sqrt( sum(D.^2) ), [ww 1] );
    % sparse coding
    disp('-> Sparse coding.');
    if strcmp(sparse_coder, 'omp')
        X = perform_omp(D,Y,options);
    else
        X = perform_mp(D,Y,options);
    end
    X = full(X);
    err(end+1) = norm(Y-D*X, 'fro');
    % update dictionary
    if not(use_cg)
        % pseudo inverion
        D = Y * pinv(X);
    else
        options.x = D';
        [D,err1] = perform_conjugate_gradient(X*X',X*Y',options);
        D = D';
    end
    % reconstruct
    Dimg = 1/ww * Z'*D(:);
end
Dimg = reshape(Dimg,k,k);
D = reshape( Z*Dimg(:), ww, dd);
D = D - repmat( mean(D), [ww 1] );
D = D ./ repmat( sqrt( sum(D.^2) ), [ww 1] );

% add the low pass DC
D = [ones(ww,1)/sqrt(ww) D];
X = [mu/sqrt(ww); X];

author: gpeyre

repo: matlab-toolboxes


text:
% plotcurve() - plot curve(s) with optional significance highlighting.
%
% Usage: >> plotcurve(times, data, 'key1', 'val1', 'key2', val2' ...);
%
% Required inputs:
%   times - [float] vector of time indices
%   data  - [float] data array, size of [n x times]. If n>1 several
%           curves are plotted (unless 'plotmean' option is used).
%
% Optional inputs:
%  'maskarray' = Input bootstrap limits. Can be 1-D [min max], 2-D (min
%                and max for all ordinate) or 3-D (min and max for all
%                ordinates at all time points).
%  'val2mask'  = Value to use for generating mask. By default use data.
%  'plotmean' = ['on'|'off'] For 'curve' plots only. Average all
%                frequencies given as input. Default: 'on'.
%  'highlightmode'  = ['background'|'bottom'] For 'curve' plots only,
%                display significant time regions either in the plot
%                background or underneatht the curve.
%  'xlabel'    = [string] x label
%  'ylabel'    = [string] y label
%  'legend'    = [cell] legend. Cell array of string, with one string
%                per curve.
%  'ylim'      = [min max] or [min] limits for the ordinate axis.
%  'title'     = Optional figure title. If two conditions are given
%                as input, title can be a cell array with two text
%                string elements {none}
%  'vert'      = Latencies to mark with a dotted vertical line   {none}
%  'linewidth' = Line width for marktimes traces (thick=2, thin=1) {2}
%  'chanlocs'  = channel location structure.
%  'plottopo'  = [min max] plot topography within the time limits defined
%                in this function. If several lines are given as input, one
%                scalp map is plot for each line.
%
% Authors: Arnaud Delorme, 2004, Bhaktivedanta Institute

% Copyright (C) 2004  Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function plotcurve( times, R, varargin);

	if nargin < 2
        help plotcurve;
        return;
	end;
	g = finputcheck( varargin, { 'maskarray' ''        []       [];
                             'val2mask'      'real'    []           R;
                             'highlightmode' 'string'  { 'background' 'bottom' } 'background';
                             'plotmean'      'string'  { 'on' 'off' }            'off';
                             'plotindiv'     'string'  { 'on' 'off' }            'on';
                             'logpval'       'string'  { 'on' 'off' }            'off';
                             'title'         'string'  []                        '';
                             'xlabel'        'string'  []                        '';
                             'plotmode'      'string'  {'single' 'topo'}         'single';
                             'ylabel'        'string'  []                        '';
                             'legend'        'cell'    []                        {};
                             'colors'        'cell'    []                        {};
                             'plottopotitle' 'cell'    []                        {};
                             'chanlocs'      'struct'  []                        struct;
                             'ylim'          'real'    []                        [];
                             'vert'          'real'    []                        [];
                             'plottopo'      'real'    []                        [];
                             'linewidth'     'real'    []                        2;
                             'marktimes'     'real'    []                        [] });
   if isstr(g), error(g); end;
  % keyboard;
   if isempty(g.colors), g.colors = { 'r' 'g' 'b' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' ...
                   'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' }; end;
   if strcmpi(g.plotindiv, 'off'), g.plotmean = 'on'; end;

   if ~any(length(times) == size(R))
       try,
           R = reshape(R, length(times), length(R)/length(times))';
       catch, error('Size of time input and array input does not match');
       end;
   end;

   % regions of significance
   % -----------------------
   if ~isempty(g.maskarray)
       if length(unique(g.maskarray)) < 4
          Rregions = g.maskarray;
       else
          Rregions = ones(size(g.val2mask));
          switch dims(g.maskarray)
               case 3, Rregions  (find(g.val2mask > g.maskarray(:,:,1) & (g.val2mask < g.maskarray(:,:,2)))) = 0;
               case 2, if size(g.val2mask,2) == size(g.maskarray,2)
                           Rregions  (find(g.val2mask < g.maskarray)) = 0;
                       elseif size(g.val2mask,1) == size(g.maskarray,1)
                           Rregions(find((g.val2mask > repmat(g.maskarray(:,1),[1 length(times)])) ...
                             & (g.val2mask < repmat(g.maskarray(:,2),[1 length(times)])))) = 0;
                       else
                           Rregions(find((g.val2mask > repmat(g.maskarray(:,1),[1 length(times)])) ...
                             & (g.val2mask < repmat(g.maskarray(:,2),[1 length(times)])))) = 0;
                       end;
               case 1, Rregions  (find(g.val2mask < repmat(g.maskarray(:),[1 size(g.val2mask,2)]))) = 0;
           end;
           Rregions = sum(Rregions,1);
       end;
   else
       Rregions = [];
   end

  % plotting
  % --------
  if size(R,1) == length(times), R = R'; end;
  if strcmpi(g.plotmean, 'on') | strcmpi(g.plotindiv, 'off')
      if strcmpi(g.plotindiv, 'on')
          R = [ R; mean(R,1) ];
      else
          R = mean(R,1);
      end;
  end;
  ax = gca;
  if ~isempty(g.maskarray) & strcmpi(g.highlightmode, 'bottom')
      pos = get(gca, 'position');
      set(gca, 'position', [ pos(1)+pos(3)*0.1 pos(2)+pos(4)*0.1 pos(3)*0.9 pos(4)*0.85 ]);
  end;

  % plot topographies
  % -----------------
  if ~isempty(g.plottopo)
      tmpax = gca;
      pos = get(gca, 'position');
      set(gca, 'position', [ pos(1) pos(2) pos(3) pos(4)/2 ]);

      for index = 1:size(g.plottopo)
          axes('position', [ (index-1)*pos(3)/size(g.plottopo,1)+pos(1) pos(2)+pos(4)/2 pos(3)/size(g.plottopo,1) pos(4)/2 ]);
          %topoplot(g.plottopo(index,:), g.chanlocs, 'maplimits', 'minmax');
          topoplot(g.plottopo(index,:), g.chanlocs);
          if ~isempty(g.plottopotitle)
              title(g.plottopotitle{index}, 'interpreter', 'none');
          end;
      end;

      axes(tmpax);
  end;

  for ind = 1:size(R,1)
      if ind == size(R,1) & strcmpi(g.plotmean, 'on') & size(R,1) > 1
           plot(times,R(ind,:), 'k', 'linewidth', 2);
      elseif ~isempty(g.colors),
           tmp = plot(times,R(ind,:), 'k');
           set(tmp, 'color', g.colors{mod(ind-1, length(g.colors))+1});

           % change the line style when number of plots exceed number of colors in g.colors
           lineStyles = {'-', '--',':','-.'};
           set(tmp,'LineStyle',lineStyles{min(ceil(ind/length(g.colors)),length(lineStyles))});

           hold on;
      else plot(times,R(ind,:));
      end;
  end;

  % ordinate limits
  % ---------------
  if isempty(g.ylim),
      yll = min(reshape(R, [1 prod(size(R))]));
      ylh = max(reshape(R, [1 prod(size(R))]));
      yll2 = yll - (ylh-yll)/10;
      ylh2 = ylh + (ylh-yll)/10;
      if ~isnan(yll), g.ylim = [yll2 ylh2]; end;
  end;
  if ~isempty(g.ylim) & length(g.ylim) == 2
      if any(g.ylim)
          ylim(g.ylim);
      else
          ylim([0 1]);
          axis off;
          box off;
      end;
  elseif ~isempty(g.ylim)
      yl = ylim;
      ylim([g.ylim yl(2)]);
  end
  yl = ylim;

  % highlight regions
  % -----------------
  if ~isempty(g.maskarray)
      axsignif = highlight(ax, times, Rregions, g.highlightmode, g.xlabel);

      % replot data (above highlighted regions)
      % ---------
      axes(ax);
      for ind = 1:size(R,1)
          if ind == size(R,1) & strcmpi(g.plotmean, 'on') & size(R,1) > 1
               plot(times,R(ind,:), 'k', 'linewidth', 2);
          elseif ~isempty(g.colors),
              tmp = plot(times,R(ind,:), 'k'); set(tmp, 'color', g.colors{mod(ind-1, length(g.colors))+1} ); hold on;
          else plot(times,R(ind,:));
          end;
      end;
      if strcmpi(g.highlightmode, 'bottom'), xlabel(''); set(ax, 'xtick', []); end;
  end;
  box on;

  ylim(yl);
  if strcmpi(g.logpval, 'on')
      set(gca, 'ytickmode', 'manual', 'yticklabel', round(10.^-get(gca, 'ytick')*1000)/1000, 'ydir', 'reverse');
  end;

  % vertical lines
  % --------------
  hold on
  xl = xlim;
  if ~isnan(g.marktimes) % plot marked time
      for mt = g.marktimes(:)'
          plot([mt mt],[yl(1) yl(2)],'--k','LineWidth',g.linewidth);
      end
  end
  hold off
  if ~isempty(g.vert)
      for index = 1:length(g.vert)
          line([g.vert(index), g.vert(index)], [yl(1) yl(2)], 'linewidth', 1, 'color', 'm');
      end;
  end;
  xlim([times(1) times(end)]);

  % title and legend
  % ----------------
  if strcmpi(g.plotmode, 'topo') % plot in scalp array
      NAME_OFFSETX = 0.1;
      NAME_OFFSETY = 0.2;
      xx = xlim; xmin = xx(1); xdiff = xx(2)-xx(1); xpos = double(xmin+NAME_OFFSETX*xdiff);
      yy = ylim; ymax = yy(2); ydiff = yy(2)-yy(1); ypos = double(ymax-NAME_OFFSETY*ydiff);
      t=text(xpos, ypos,g.title);
      axis off;
      line([0 0], [yl(1) yl(2)], 'linewidth', 1, 'color', 'k');
      line([xl(1) xl(2)], [0 0], 'linewidth', 1, 'color', 'k');
      set(ax, 'userdata', { g.xlabel g.ylabel g.legend });
  else
      title(g.title, 'interpreter', 'none')
      if ~isempty(g.legend)
          hh = legend(g.legend(:));
          set(hh, 'unit', 'pixels', 'interpreter', 'none')
      end;
      if isempty(g.maskarray)
          xlabel(g.xlabel);
      end;
      ylabel(g.ylabel)
  end;

% -----------------
% highlight regions
% -----------------
function axsignif = highlight(ax, times, regions, highlightmode, myxlabel);
color1 = [0.75 0.75 0.75];
color2 = [0 0 0];
yl  = ylim;
yl(1) = yl(1)-max(abs(yl));
yl(2) = yl(2)+max(abs(yl));

if ~strcmpi(highlightmode, 'background')
    pos = get(ax, 'position');
    set(gca, 'xtick', []);
    axsignif = axes('position', [pos(1) pos(2)-pos(4)*0.05 pos(3) pos(4)*0.05 ]);
    plot(times, times, 'w');
    set(axsignif, 'ytick', []);
    yl2 = ylim;
    yl2(1) = yl2(1)-max(abs(yl2));
    yl2(2) = yl2(2)+max(abs(yl2));
    xlim([times(1) times(end)]);
    xlabel(myxlabel);
else
    axsignif = [];
    xlabel(myxlabel);
end;

if ~isempty(regions)
    axes(ax);
    in_a_region = 0;
    for index=1:length(regions)
        if regions(index) & ~in_a_region
            tmpreg(1) = times(index);
            in_a_region = 1;
        end;
        if (~regions(index) | index == length(regions)) & in_a_region
            tmpreg(2) = times(index);
            in_a_region = 0;
            if strcmpi(highlightmode, 'background')
                tmph = patch([tmpreg(1) tmpreg(2) tmpreg(2) tmpreg(1)], ...
                    [yl(1) yl(1) yl(2) yl(2)], color1); hold on;
                set(tmph, 'edgecolor', color1);
            else
                oldax = ax;
                axes(axsignif);
                tmph = patch([tmpreg(1) tmpreg(2) tmpreg(2) tmpreg(1)], ...
                    [yl2(1) yl2(1) yl2(2) yl2(2)], color2); hold on;
                set(tmph, 'edgecolor', color2);
                axes(oldax);
            end;
        end;
    end;
    ylim(yl);
end;

  function res = dims(array)
    res = min(ndims(array), max(size(array,2),size(array,3)));

author: PatternRecognition

repo: OpenBMI


text:
function x = real2int_vbmc(x,trinfo,integervars)
%REAL2INT_VBMC Convert to integer-valued representation.

if ~any(integervars); return; end

xtemp = warpvars_vbmc(x,'inv',trinfo);
xtemp(:,integervars) = round(xtemp(:,integervars));
xtemp = warpvars_vbmc(xtemp,'d',trinfo);

x(:,integervars) = xtemp(:,integervars);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
%DELLABLIST Delete a label list from dataset
%
%	B = DELLABLIST(A,LABLISTNAME)
%	B = DELLABLIST(A,LABLISTNUMBER)
%
% INPUT
%   A              - Dataset
%   LABLISTNAME    - String to identify the label list to be deleted
%   LABLISTNUMBER  - Number to identify the label list to be deleted
%
% OUTPUT
%   B              - Dataset
%
% DESCRIPTION
% In the multi-label system for datasets, additional labellings can be
% added by ADDLABELS. This is stored in the LABLIST and NLAB fields of
% the dataset. By this command, DELLABLIST such a labeling can be removed.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MULTI_LABELING, ADDLABELS, ADDLABLIST, CHANGELABLIST, CURLABLIST

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% 5g ldpc encoding
% input: s, bit sequence of dimension K * 1
% output: encoded_bits, bit sequence
% reference: 3GPP TS 38.212 section 5.3.2
% author: Xiao, Shaoning 萧少宁
% license: MIT

function [encoded_bits, H, Z_c, encoded_bits_original] = ldpc_encode(s, base_graph_index)

K = length(s);

encoded_bits = zeros(3*K, 1);

if base_graph_index == 1
    a = 4;
    b = 22;
    c = 26;
    d = 42;
    e = 46;
    z = K/b;
    Z_c = z;
    N = 66 * Z_c;
    z = K/b;
    set_index = lifting_size_table_lookup(z);
    load parity_check_matrices_protocol_1
    BG = parity_check_matrices_protocol_1(:, :, set_index); %#ok<NODEF>
elseif base_graph_index == 2
    a = 4;
    b = 10;
    c = 14;
    d = 38;
    e = 42;
    z = K/b;
    Z_c = z;
    N = 50 * Z_c;
    set_index = lifting_size_table_lookup(z);
    load parity_check_matrices_protocol_2
    BG = parity_check_matrices_protocol_2(:, :, set_index); %#ok<NODEF>
else
  error('wrong base graph index in ldpc encoding.');
end

BG(BG ~= -1) = mod(BG(BG ~= -1), Z_c);

for k = (2*Z_c):(K-1)
  if s(k+1) ~= -1
    encoded_bits(k-2*Z_c+1) = s(k+1);
  else
    s(k+1) = 0;
    encoded_bits(k-2*Z_c+1) = -1;
  end
end

% set_index = lifting_size_table_lookup(Z_c);
% BG = parity_check_matrices_protocol(:, :, set_index);

A_prime = BG(1:a, 1:b);
B_prime = BG(1:a, (b+1):c);
C_prime = BG((a+1):e, 1:b);
D_prime = BG((a+1):e, (b+1):c);

z = Z_c;

A = spalloc(a*z, b*z, nnz(A_prime + ones(size(A_prime))));

for row_index = 1:a
    for column_index = 1:b
        if A_prime(row_index, column_index) ~= -1
            A((row_index-1)*z+1:row_index*z, (column_index-1)*z+1:column_index*z) = sparse(1:z, [(mod(A_prime(row_index, column_index), z)+1):z, 1:mod(A_prime(row_index, column_index), z)], ones(1, z), z, z);
        end
    end
end

B = spalloc(a*z, a*z, nnz(B_prime + ones(size(B_prime))));

for row_index = 1:a
    for column_index = 1:a
        if B_prime(row_index, column_index) ~= -1
            B((row_index-1)*z+1:row_index*z, (column_index-1)*z+1:column_index*z) = sparse(1:z, [(mod(B_prime(row_index, column_index), z)+1):z, 1:mod(B_prime(row_index, column_index), z)], ones(1, z), z, z);
        end
    end
end

C = spalloc(d*z, b*z, nnz(C_prime + ones(size(C_prime))));

for row_index = 1:d
    for column_index = 1:b
        if C_prime(row_index, column_index) ~= -1
            C((row_index-1)*z+1:row_index*z, (column_index-1)*z+1:column_index*z) = sparse(1:z, [(mod(C_prime(row_index, column_index), z)+1):z, 1:mod(C_prime(row_index, column_index), z)], ones(1, z), z, z);
        end
    end
end

D = spalloc(d*z, a*z, nnz(D_prime + ones(size(D_prime))));

for row_index = 1:d
    for column_index = 1:a
        if D_prime(row_index, column_index) ~= -1
            D((row_index-1)*z+1:row_index*z, (column_index-1)*z+1:column_index*z) = sparse(1:z, [(mod(D_prime(row_index, column_index), z)+1):z, 1:mod(D_prime(row_index, column_index), z)], ones(1, z), z, z);
        else
            D((row_index-1)*z+1:row_index*z, (column_index-1)*z+1:column_index*z) = spalloc(z, z, 0);
        end
    end
end

B_inv = spalloc(a*z, a*z, 20*z);

if (base_graph_index == 1) && (set_index ~= 7)
    B_inv(1:z, 1:z)             = speye(z);
    B_inv(1:z, 1+z:2*z)         = speye(z);
    B_inv(1:z, 1+2*z:3*z)       = speye(z);
    B_inv(1:z, 1+3*z:4*z)       = speye(z);
    B_inv(1+z:2*z, 1:z)         = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+z:2*z)     = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+2*z:3*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+3*z:4*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1:z)       = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+z:2*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+2*z:3*z) = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1:z)       = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+z:2*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+2*z:3*z) = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
elseif (base_graph_index == 2) && ((set_index ~= 4) && (set_index ~= 8))
    B_inv(1:z, 1:z)             = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+z:2*z)         = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+2*z:3*z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+3*z:4*z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1:z)         = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+z:2*z)     = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+2*z:3*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+3*z:4*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1:z)       = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+z:2*z)   = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+2*z:3*z) = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+3*z:4*z) = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1:z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+z:2*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+2*z:3*z) = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
elseif (base_graph_index == 1) && (z == 208)
    B_inv(1:z, 1:z)             = sparse(circshift(eye(208), 105));
    B_inv(1:z, 1+z:2*z)         = sparse(circshift(eye(208), 105));
    B_inv(1:z, 1+2*z:3*z)       = sparse(circshift(eye(208), 105));
    B_inv(1:z, 1+3*z:4*z)       = sparse(circshift(eye(208), 105));
    B_inv(1+z:2*z, 1:z)         = speye(208) + sparse(circshift(eye(208), 105));
    B_inv(1+z:2*z, 1+z:2*z)     = sparse(circshift(eye(208), 105));
    B_inv(1+z:2*z, 1+2*z:3*z)   = sparse(circshift(eye(208), 105));
    B_inv(1+z:2*z, 1+3*z:4*z)   = sparse(circshift(eye(208), 105));
    B_inv(1+2*z:3*z, 1:z)       = sparse(circshift(eye(208), 105));
    B_inv(1+2*z:3*z, 1+z:2*z)   = sparse(circshift(eye(208), 105));
    B_inv(1+2*z:3*z, 1+2*z:3*z) = speye(208) + sparse(circshift(eye(208), 105));
    B_inv(1+2*z:3*z, 1+3*z:4*z) = speye(208) + sparse(circshift(eye(208), 105));
    B_inv(1+3*z:4*z, 1:z)       = sparse(circshift(eye(208), 105));
    B_inv(1+3*z:4*z, 1+z:2*z)   = sparse(circshift(eye(208), 105));
    B_inv(1+3*z:4*z, 1+2*z:3*z) = sparse(circshift(eye(208), 105));
    B_inv(1+3*z:4*z, 1+3*z:4*z) = speye(208) + sparse(circshift(eye(208), 105));
elseif (base_graph_index == 1) && ((z ~= 208) && (set_index == 7))
    B_inv(1:z, 1:z)             = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+z:2*z)         = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+2*z:3*z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1:z, 1+3*z:4*z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1:z)         = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+z:2*z)     = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+2*z:3*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+3*z:4*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1:z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+z:2*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+2*z:3*z) = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1:z)       = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+z:2*z)   = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+2*z:3*z) = sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [z, 1:(z-1)], ones(1, z), z, z);
elseif (base_graph_index == 2) && ((set_index == 4) || (set_index == 8))
    B_inv(1:z, 1:z)             = speye(z);
    B_inv(1:z, 1+z:2*z)         = speye(z);
    B_inv(1:z, 1+2*z:3*z)       = speye(z);
    B_inv(1:z, 1+3*z:4*z)       = speye(z);
    B_inv(1+z:2*z, 1:z)         = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+z:2*z)     = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+2*z:3*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+z:2*z, 1+3*z:4*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1:z)       = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+z:2*z)   = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+2*z:3*z) = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+2*z:3*z, 1+3*z:4*z) = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1:z)       = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+z:2*z)   = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+2*z:3*z) = sparse(1:z, [2:z, 1], ones(1, z), z, z);
    B_inv(1+3*z:4*z, 1+3*z:4*z) = speye(z) + sparse(1:z, [2:z, 1], ones(1, z), z, z);
end

s = s(:);

p_1 = mod(B_inv * (A * s), 2);
p_2 = mod(C * s + D * p_1, 2);

w = [p_1; p_2];

for k = K:(N+2*Z_c-1)
   encoded_bits(k-2*Z_c+1) = w(k-K+1);
end

H = [A, B, spalloc(a*z, d*z, 0); C, D, speye(d*z)];

encoded_bits_original = [s; w];

% mod(H * [s; p_1; p_2]) = 0

clear A
clear B
clear C
clear D
clear B_inv
clear A_prime
clear B_prime
clear C_prime
clear D_prime

end

author: xiaoshaoning


text:
function pass = test_PTdecomposition( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e5*pref.techPrefs.chebfuneps;

% Example 1 :
p = ballfun(@(r,lam,th)cos(r.^2.*sin(th).^2.*cos(lam).*sin(lam)), 'spherical');
t = ballfun(@(r,lam,th)sin(r.^2.*sin(th).*cos(th).*sin(lam)), 'spherical');
V = ballfunv.PT2ballfunv(p,t);
[p2, t2] = PTdecomposition(V);
pass(1) = norm(diff(p,2,'spherical')-diff(p2,2,'spherical'))<tol;
pass(2) = norm(diff(p,3,'spherical')-diff(p2,3,'spherical'))<tol;
pass(3) = norm(diff(t,2,'spherical')-diff(t2,2,'spherical'))<tol;
pass(4) = norm(diff(t,3,'spherical')-diff(t2,3,'spherical'))<tol;

% Example 2 :
p = ballfun(@(x,y,z)x.^2+y.*z);
t = ballfun(@(x,y,z)x.*y.*z);
V = ballfunv.PT2ballfunv(p,t);
[p2, t2] = PTdecomposition(V);
pass(5) = norm(diff(p,2,'spherical')-diff(p2,2,'spherical'))<tol;
pass(6) = norm(diff(p,3,'spherical')-diff(p2,3,'spherical'))<tol;
pass(7) = norm(diff(t,2,'spherical')-diff(t2,2,'spherical'))<tol;
pass(8) = norm(diff(t,3,'spherical')-diff(t2,3,'spherical'))<tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [MSE, SNR, PSNR] = Calc_MSE_SNR(I1,I2,b)
% This function computes MSE, SNR and PSNR between two gray
% images excluding a border of width b on all four sides.
% The two gray images must have the same size.
%
%          I1(:,:)      first input original image, range 0~255
%          I2(:,:)      second input observed image, range 0~255
%          b            width of border on each side of the image to exclude
%                       from error calculations
%          MSE          output, mean square error
%          SNR          output, signal noise ratio
%          PSNR         output, peak signal noise ratio
%
% Copyright (c) Apr., 2006. Dengwen Zhou. All rights reserved.
% Department of Computer Science & Technology
% North China Electric Power University(NCEPU)
%
% Last time modified: Oct. 11, 2012
%

% Exclude b border pixels
I1 = double(I1); I2 = double(I2);
I1 = I1(b+1:end-b,b+1:end-b);
I2 = I2(b+1:end-b,b+1:end-b);

% Compute the errors
num = numel(I1);
d = sum((I1(:)-I2(:)).^2);
s = sum(I1(:).^2);
MSE = d/num;
SNR = 10*log10(s/d);
PSNR = 10*log10(255*255/MSE);


text:
function fixed = p06_fixed_points ( m, fixed_num )

%*****************************************************************************80
%
%% P06_FIXED_POINTS returns the fixed points in problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the fixed points.
%
  fixed = [];

  return
end


text:
%% 2D Heat Transfer D2Q4
% by Manuel Diaz
clc; clear; %close all;

%% Domain Grid
L = 30;    % Length
H = 100;    % Height
n = 30;    % Nodes in x
m = 100;    % Nodes in y
dx = L/n; dy = H/m; dt = 1;
x = 1:dx:L; y = 1:dy:H;

%% Initial Condition
f1 = zeros(m,n);
f2 = zeros(m,n);
f3 = zeros(m,n);
f4 = zeros(m,n);
rho= zeros(m,n); % initial value of the dependent variable rho = T(x,t)
feq= zeros(m,n);

%% Constants
csq   = (dx^2)/(dt^2);
alpha = 0.25;
omega = 1/(2*alpha/(dt*csq)+0.5);
w     = [1/4,1/4,1/4,1/4];
cx    = [  1, -1,  0,  0];
cy    = [  0,  0,  1, -1];
link  = [  1,  2,  3,  4];
tEnd  = 400; % time steps
tPlot = 10;
frames= tEnd/tPlot;
count = 0;

%% Movie Parameters
tPlot = 10; frames= tEnd/tPlot; count = 0;
figure(1)
colordef white %black

%% Boundary Conditions
twall = 1;

%% Main Loop
for cycle = 2:tEnd;
    % Collision process
    rho = f1 + f2 + f3 + f4;

    % for this case k1=k2=k3=k4=1/4 then feq1=...feq4=feq
    feq = 0.25 * rho;

    f1 = w(1) * rho;
    f2 = w(2) * rho;
    f3 = w(3) * rho;
    f4 = w(4) * rho;

    f1 = omega * feq + (1-omega) * f1;
    f2 = omega * feq + (1-omega) * f2;
    f3 = omega * feq + (1-omega) * f3;
    f4 = omega * feq + (1-omega) * f4;

    % Streaming process
    f1 = stream2d( f1 , [cy(1),cx(1)]);
    f2 = stream2d( f2 , [cy(2),cx(2)]);
    f3 = stream2d( f3 , [cy(3),cx(3)]);
    f4 = stream2d( f4 , [cy(4),cx(4)]);

    % Boundary conditions
    f1(:,1) = twall*w(1) + twall*w(2) - f2(:,1); %Dirichlet BC
    f3(:,1) = twall*w(3) + twall*w(4) - f4(:,1); %Dirichlet BC

    f1(m,:) = 0;        %Dirichlet BC
    f2(m,:) = 0;        %Dirichlet BC
    f3(m,:) = 0;        %Dirichlet BC
    f4(m,:) = 0;        %Dirichlet BC

    f1(:,n) = 0;        %Dirichlet BC
    f2(:,n) = 0;        %Dirichlet BC
    f3(:,n) = 0;        %Dirichlet BC
    f4(:,n) = 0;        %Dirichlet BC

    f1(1,:) = f1(2,:);  %Neumann BC
    f2(1,:) = f2(2,:);  %Neumann BC
    f3(1,:) = f3(2,:);  %Neumann BC
    f4(1,:) = f4(2,:);  %Neumann BC

    % Visualization
    if mod(cycle,tPlot) == 1
        count = count + 1;
        contourf(rho)
        colormap hot
        colorbar('location','southoutside')
        M(count)=getframe;
    end

    % Animated gif file
    if mod(cycle,tPlot) == 1
        F = getframe;
        if count == 1
            [im,map] = rgb2ind(F.cdata,256,'nodither');
            im(1,1,1,tEnd/tPlot) = 0;
        end
        im(:,:,1,count) = rgb2ind(F.cdata,map,'nodither');
    end
end

%% Make Movie
movie(M,1,10); % movie(M,n,fps)

%% Export to Gif
imwrite(im,map,'heat_eq_d2q4.gif','DelayTime',0,'LoopCount',3)

repo: Aero-matlab


text:
function res = is_atom(x)
% Determine whether an expression can be further subdivided into parts
% literal matlab expressions (including matrices, most structs, strings, function_handles) are considered atomic
%
% The exceptions are:
%  * struct with .head / .part fields      (canonically represented expression)
%  * struct with .tracking.expression field  (impure expression, i.e. expression with value)
%
% Notes:
%  in the future, more data structures may be made accessible to the expression system
%
% See also:
%   exp_eval, exp_head, exp_parts
%
%                                                     Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                                     2010-04-15

if isfield(x,{'head','parts'})
    res = true;
elseif isfield(x,'tracking') && isfield(x.tracking,'expression')
    res = is_atom(x.tracking.expression);
else
    res = false;
end

repo: REST


text:
% return a get batch function
% -------------------------------------------------------------------------
function fn = getBatchWrapper_augVOC2012(opts)
fn = @(images, mode) getBatch_dict(images, mode, opts) ;
end

function [imBatch, semanticMaskBatch, instanceMaskBatch, weightBatch] = getBatch_dict(images, mode, opts)
[imBatch, semanticMaskBatch, instanceMaskBatch, weightBatch] = getImgBatch_augVOC2012(images, mode, opts, 'prefetch', nargout == 0) ;
end

author: aimerykong


text:
function rrd = jpleph_mice (et, ntarg, ncent)

% reads the jpl planetary ephemeris and gives the position and velocity
% of the point 'ntarg' with respect to point 'ncent' using MICE routines

% input

%   et    = TDB julian date at which interpolation is wanted

%   ntarg = integer number of 'target' point

%   ncent = integer number of center point

%   the numbering convention for 'ntarg' and 'ncent' is:

%        1 = mercury           8 = neptune
%        2 = venus             9 = pluto
%        3 = earth            10 = moon
%        4 = mars             11 = sun
%        5 = jupiter
%        6 = saturn
%        7 = uranus

% output

%   rrd = output 6-word array containing position and velocity
%         of point 'ntarg' relative to 'ncent'. the units are
%         determined by the value of km passed via global.

% global

%   iephem  = initialization flag (1 = initialize)
%   ephname = name of ephemeris binary data file (de436.bsp, etc.)
%   km      = state vector units flag (1 = km & km/sec, 0 = au & au/day)
%   au      = numerical value of astronomical unit (kilometers)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global iephem ephname km au

if (iephem == 1)

    % load binary ephemeris data file

    go_dir = Core.getGlobalConfig().getLocalStorageDir();
    cspice_furnsh(fullfile(go_dir, ephname));

    % reset initialization flag

    iephem = 0;

end

% set name of target body

switch ntarg

    case (1)

        targ = 'mercury';

    case (2)

        targ = 'venus';

    case (3)

        targ = 'earth';

    case (4)

        targ = 'mars';

    case(5)

        targ = 'jupiter';

    case (6)

        targ = 'saturn';

    case (7)

        targ = 'neptune';

    case (8)

        targ = 'uranus';

    case (9)

        targ = 'pluto';

    case (10)

        targ = 'moon';

    case (11)

        targ = 'sun';

end

% set name of central body

switch ncent

    case (1)

        obs = 'mercury';

    case (2)

        obs = 'venus';

    case (3)

        obs = 'earth';

    case (4)

        obs = 'mars';

    case (5)

        obs = 'jupiter';

    case (6)

        obs = 'saturn';

    case (7)

        obs = 'neptune';

    case (8)

        obs = 'uranus';

    case (9)

        obs = 'pluto';

    case (10)

        obs = 'moon';

    case (11)

        obs = 'sun';

end

% compute time, expressed as TDB seconds past J2000 TDB (2451545.0)

etime = 86400.0d0 * (et - 2451545.0d0);

% compute position and velocity vectors in eme2000 system (no corrections)

starg = mice_spkezr(targ, etime, 'J2000', 'NONE', obs);

% provide output in user-requested units

if (km == 1)

    % state is kilometers and kilometers/second

    rrd = starg.state;

else

    % state is au's and au's/day

    rrd(1:3) = starg.state(1:3) / au;

    rrd(4:6) = 86400.0 * starg.state(4:6) / au;

    rrd = rrd';

end

author: goGPS-Project


text:
function vardist = vardistCreate(X, Q, type, constraint)

% VARDISTCREATE description.

% VARGPLVM

%
% creates the structure of the variational distirbution over the latent
% values in the GP-LVM
%
% The variational distribution is assumed to be factorized over the
% the latent variables. Each factor each a Gaussain N(x_n|mu_n,S_n)
% with a diagonal covariance matrix
%
% The structure of the vardist is similar to the structure of a kernel

if nargin == 3
    constraint = optimiDefaultConstraint('positive');
else
    constraint = 'identity';
end

vardist.type = 'vardist';
vardist.vartype = type;

vardist.numData = size(X,1);
vardist.latentDimension = Q;
vardist.nParams = 2*vardist.numData*vardist.latentDimension;
%  number of training points
N = size(X,1);

vardist.transforms(1).index = [(N*Q+1):vardist.nParams];
vardist.transforms(1).type = constraint;

% initialize the parameters
vardist.means  = X;
vardist.covars = 0.1*ones(N,Q) + 0.001*randn(N,Q);
vardist.covars(vardist.covars<0.05) = 0.05;
%vardist.covars = (eps^2)*ones(size(vardist.covars));
%vardist.means = randn(N,Q);
%pmeans = randn(Q,N);
%pcovars = randn(Q,N);
%params = [pmeans(:)' pcovars(:)'];
%vardist = vardistExpandParam(vardist, params);

author: SheffieldML

repo: vargplvm


text:
function desc = calGradient(rgb_im, seg, numRegion)

    G = fspecial('gaussian',[4 4],2);
    Ig = imfilter(rgb_im,G,'same');
    [Gmag,Gdir] = imgradient(rgb2gray(Ig));

    binNum = 20;
    inter = max(max(Gmag)) / binNum;
    Gmag = int32(Gmag / inter);

    binVal = 1:binNum;
    desc = zeros([numRegion binNum]);

    cnt = 0;
    ind={};
    for iReg=1:numRegion
        ind{iReg} = seg(:)==iReg;
    end

    for bin=1:binNum
        cnt = cnt + 1;
        I =  (Gmag(:)==binVal(bin)) ;
        for iReg=1:numRegion
            desc(iReg,cnt) = sum(I(ind{iReg}));
        end
    end

    tmp = sum(desc, 2);
    desc = desc ./ repmat(tmp(:,:), [1 size(desc,2)]);

end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
%%% Replay the set of solution vectors:

if ~exist('param_list'),
   if ~exist('solution');
      fprintf(1,'Error: Need to calibrate first\n');
      return;
   else
      param_list = solution;
   end;
end;

N_iter = size(param_list,2);

if N_iter == 1,
   fprintf(1,'Warning: There is a unique state in the list of parameters.\n');
end;

%M = moviein(N_iter);

for nn = 1:N_iter,

   solution = param_list(:,nn);

   extract_parameters;
   comp_error_calib;

   ext_calib;

   drawnow;

%   Mnn = getframe(gcf);

%   M(:,nn) = Mnn;

end;

%fig = gcf;

%figure(fig+1);
%close;
%figure(fig+1);

%clf;
%movie(M,20);

repo: ekfmonoslam


text:
function [ mO ] = ApplyLocalLinearKernel( mI, filterRadius, regFctr )
% ----------------------------------------------------------------------------------------------- %
%[ mOutputImage ] = ApplyLocalLinearFilter( mInputImage, filterRadius, regFctr )
% Applying Linear Edge Preserving Smoothing Filter based on Local Linear (Affine) model.
% Input:
%   - mI            -   Input Image.
%                       Structure: Image Matrix (Single Channel).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
%   - filterRadius  -   Filter Radius.
%                       The filter radius.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ...}.
%   - regFctr       -   Regularization Factor.
%                       Regularize the local covariance (Variance).
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
% Output:
%   - mO            -   Output Image.
%                       Structure: Image Matrix (Single Channel).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
% References:
%   1.  "Guided Image Filtering".
% Remarks:
%   1.  This is basically estimating the Linear Function (Affine)
%       parameters for the Local Window. Namely, the ouput is a Linear
%       combination of the input Window and a DC Factor. The final step is
%       aggregation (Uniform) off all estimations of the parameters.
%   2.  Prefixes:
%       -   'v' - Vector.
%       -   'm' - Matrix.
%       -   't' - Tensor (Multi Dimension Matrix)
%       -   's' - Struct.
%       -   'c' - Cell Array.
%   3.  The calculation of the Local Variance might be negative due to
%       numerical diffuculties. If artifacts appear, this might be the
%       casue. Usually using matrices of type 'double' solves it.
%   4.  This implemnetation is `ApplyGuidedFilter` where the Guiding Image
%       is the Input Image.
%   5.  Speed otimizzation can be achived by wiser use of 'mNumEffPixels'.
%       Instead of dividing by it calculate its reciprocal once. Moreover,
%       it can be used only once in the aggregation process.
% TODO:
%   1.  Create Multi Variable Linear Model.
%   2.  Some speed potimization could be made (Taking advantage of 'mNumEffPixels').
% Release Notes:
%   -   1.0.000     05/01/2016  Royi Avital
%       *   First release version
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF = 0;
ON  = 1;

BORDER_TYPE_CONSTANT    = 1;
BORDER_TYPE_CIRCULAR    = 2;
BORDER_TYPE_REPLICATE   = 3;
BORDER_TYPE_SYMMETRIC   = 4;

numRows = size(mI, 1);
numCols = size(mI, 2);

borderType      = BORDER_TYPE_CONSTANT;
borderValue     = 0;
normalizeFlag   = OFF;

mNumEffPixels = ApplyBoxFilter(ones([numRows, numCols]), filterRadius, borderType, borderValue, normalizeFlag);

mLocalMean          = ApplyBoxFilter(mI, filterRadius, borderType, borderValue, normalizeFlag) ./ mNumEffPixels;
mLocalMeanSquare    = ApplyBoxFilter((mI .* mI), filterRadius, borderType, borderValue, normalizeFlag) ./ mNumEffPixels;
mLocalCovariance    = mLocalMeanSquare - (mLocalMean .* mLocalMean);

mO = zeros([numRows, numCols]);

for jj = 1:numCols
    for ii = 1:numRows
        % The coordinate i, j are set.

        sumWeights = 0;

        % Running on the Neighborhood of the pixel
        for ll = -filterRadius:filterRadius
            for kk = -filterRadius:filterRadius

                jRowIdx = ii + kk;
                jColIdx = jj + ll;

                if((jColIdx >= 1) && (jColIdx <= numCols) && (jRowIdx >= 1) && (jRowIdx <= numRows))
                    % Valid Pixel
                    numPixels       = 0;
                    jPixelWeight    = 0;

                    % Running on Wk Window
                    for nn = -filterRadius:filterRadius
                        for mm = -filterRadius:filterRadius

                            % K Index
                            kRowIdx = jRowIdx + mm;
                            kColIdx = jColIdx + nn;

                            if((kColIdx >= 1) && (kColIdx <= numCols) && (kRowIdx >= 1) && (kRowIdx <= numRows) && (abs(kColIdx - jj) <= filterRadius) && (abs(kRowIdx - ii) <= filterRadius))

                                numPixels = numPixels + 1;

                                % Weight
                                jPixelWeight = jPixelWeight + 1 + (((mI(ii, jj) - mLocalMean(kRowIdx, kColIdx)) * (mI(jRowIdx, jColIdx) - mLocalMean(kRowIdx, kColIdx))) / (mLocalCovariance(kRowIdx, kColIdx) + regFctr));
                            end
                        end
                    end

                    % jPixelWeight = jPixelWeight / (numPixels * numPixels);
                    mO(ii, jj) = mO(ii, jj) + (jPixelWeight * mI(jRowIdx, jColIdx));

                    sumWeights = sumWeights + jPixelWeight;

                end
            end
        end

        mO(ii, jj) = mO(ii, jj) / sumWeights;

    end
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%function [W]=myCTwriteVMC(inputCT,Xnew, Ynew, Znew);
function [W]=myCTwriteVMC(inputCT,scaleX, scaleY, scaleZ, str);

% Create CT input file *.ct for VMC++ engine from CERR MatLab CT matrix
% scaleX - X voxel size in cm
% scaleY - Y voxel size in cm
% scaleZ - Z voxel size in cm
% numSlice - number of slices in CT matrix
% In this code CERR CT data converted to density from HF (water has 1024).
% CT coordinate system is following that real (Xmin,Ymin,Zmin) is (0,0,0) in all cases here
% Output - *.ct file and dimension
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

CT=inputCT;
W=size(CT);

%[newfile,newpath] = uiputfile('*.ct','Enter CT file name for VMC++ input ');
%str=strcat(newpath,newfile);

fid=fopen(str,'w');

count=fwrite(fid,W(1),'int32');
count=fwrite(fid,W(2),'int32');
count=fwrite(fid,W(3),'int32');

%X=[0 Xnew];
%Y=[0 Ynew];
%Z=[0 Znew];

X(1)=0;
for i=1:W(1),
    X(i+1)=scaleX*i;
end
Y(1)=0;
for i=1:W(2),
    Y(i+1)=scaleY*i;
end
Z(1)=0;
for i=1:W(3),
    Z(i+1)=scaleZ*i;
end

count=fwrite(fid,X,'float32');
count=fwrite(fid,Y,'float32');
count=fwrite(fid,Z,'float32');
scale=1e3;
%h = waitbar(0,'Please wait for CT file writing');
density_max=4;
for k=1:W(3),
    temp=double(CT(:,:,k))/1024; % round up to 2 decimal places
    temp(temp>density_max)=density_max-1e-6; % cut values greater than 3.
    count=fwrite(fid,temp,'float32');
%    waitbar(k/W(3));
end

%close(h);
fclose(fid);

numSlice=W(3);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [x,obj,err,iter] = fusedl1(A,b,lambda,opts)

% Solve the fused Lasso (Fused L1) minimization problem by ADMM
%
% min_x ||x||_1 + lambda*\sum_{i=2}^p |x_i-x_{i-1}|,
%   s.t. Ax=b
%
% ---------------------------------------------
% Input:
%       A       -    d*n matrix
%       b       -    d*1 vector
%       lambda  -    >=0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       x       -    n*1 vector
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 20/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

[d,n] = size(A);
x = zeros(n,1);

z = x;
Y1 = zeros(d,1);
Y2 = x;

Atb = A'*b;
I = eye(n);
invAtAI = (A'*A+I)\I;

% parameters for "flsa" (from SLEP package)
tol2 = 1e-10;      % the duality gap for termination
max_step = 50;     % the maximal number of iterations
x0 = zeros(n-1,1); % the starting point

iter = 0;
for iter = 1 : max_iter
    xk = x;
    zk = z;
    % update x.
    % flsa solves min_x 1/2||x-v||_2^2+lambda1*||x||_1+lambda2*\sum_{i=2}^p |x_i-x_{i-1}|
    x = flsa(z-Y2/mu,x0,1/mu,lambda/mu,n,max_step,tol2,1,6);
    % update z
    z = invAtAI*(-A'*Y1/mu+Atb+Y2/mu+x);
    dY1 = A*z-b;
    dY2 = x-z;
    chgx = max(abs(xk-x));
    chgz = max(abs(zk-z));
    chg = max([chgx chgz max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = comp_fusedl1(x,1,lambda);
            err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y1 = Y1 + mu*dY1;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
obj = comp_fusedl1(x,1,lambda);
err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);

function f = comp_fusedl1(x,lambda1,lambda2)
% compute f = lambda1*||x||_1 + lambda2*\sum_{i=2}^p |x_i-x_{i-1}|.
% x - p*1 vector
f = 0;
p = length(x);
for i = 2 : p
   f = f+abs(x(i)-x(i-1));
end
f = lambda1*norm(x,1)+lambda2*f;

author: canyilu

repo: LibADMM-toolbox


text:
%
% Bayesian Optimization of Combinatorial Structures
%
% Copyright (C) 2018 R. Baptista & M. Poloczek
%
% BOCS 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.
%
% BOCS 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 BOCS.  If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) 2018 MIT & University of Arizona
% Authors: Ricardo Baptista & Matthias Poloczek
% E-mails: rsb@mit.edu & poloczek@email.arizona.edu
%

% Script runs discrete optimization algorithms for the
% contamination control simulation problem. The results are
% compared for different values of the \lambda tuning
% parameter.

clear; close all; clc
addpath(genpath('../algorithms'))
addpath(genpath('../stat_model'))
addpath(genpath('../test_problems/ContStudy'))
addpath(genpath('../tools'))

%% Setup parameters

% Setup fixed parameters
n_vars  = 30;
n_proc  = 1;
test_name = 'contamination';

% Number of runs and optimization iterations
n_func     = 10;
n_runs     = 10;
n_init     = 20;
evalBudget = 270;

% problem parameters (Monte Carlo samples)
mcSamps = 1e2;

% Variance prior parameters (Inverse Gamma)
aPr    = 2;
bPr    = 1;

% Regularization parameters
lambda_vals = [0, 1e-4, 1e-2, 1];
lambda_str  = {'0', '1em4', '1em2', '1'};

% Set additive regularization function
reg_term = @(x) sum(x,2);

%% Generate Test Cases

% setup objective functions
inputs_all = cell(n_func, n_runs);

for t1=1:n_func

    fprintf('Setting up test function %d\n', t1);

    % Generate random case study
    seed = randi(10000,1);

    for t2=1:n_runs

        % Set inputs struct for each problem
        inputs_all{t1,t2} = struct;
        inputs_all{t1,t2}.n_vars      = n_vars;

        % Save other definitions
        inputs_all{t1,t2}.evalBudget  = evalBudget;
        inputs_all{t1,t2}.n_runs      = n_runs;
        inputs_all{t1,t2}.n_init      = n_init;
        inputs_all{t1,t2}.lambda_vals = lambda_vals;

        % Set priors for estimator
        inputs_all{t1,t2}.aPr         = aPr;
        inputs_all{t1,t2}.bPr         = bPr;

        % Save objective function and regularization term
        inputs_all{t1,t2}.model = @(x) contamination_prob(x, mcSamps, seed);
        inputs_all{t1,t2}.reg_term = @(x) reg_term(x);

        % Generate initial samples for statistical models
        inputs_all{t1,t2}.x_vals = sample_models(n_init, n_vars);
        inputs_all{t1,t2}.y_vals = inputs_all{t1,t2}.model(inputs_all{t1,t2}.x_vals);

    end
end

inputs_all = reshape(inputs_all, n_func*n_runs, 1);

% Make folder
mkdir(['../results/' test_name])

% Save test cases
save(['../results/' test_name '/all_tests'])

% Run cases
run_cases(inputs_all, lambda_vals, test_name, n_proc);

% -- END OF FILE --

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function shoreline_test ( )

%*****************************************************************************80
%
%% SHORELINE_TEST tests the SHORELINE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHORELINE_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SHORELINE library.\n' );

  circle_centered_test ( );
  circle_offcenter_test ( );
  ellipse_centered_test ( );
  ellipse_slanted_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHORELINE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function circle_centered_test ( )

%*****************************************************************************80
%
%% CIRCLE_CENTERED_TEST uses a circle centered at the origin.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2012
%
%  Author:
%
%    John Burkardt
%
  global c
  global r

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_CENTERED_TEST\n' );
%
%  C = [0,0],  R = 1.
%
  c = [ 0.0, 0.0 ];
  r = 1.0;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [0,0],  R = 1/2.
%
  c = [ 0.0, 0.0 ];
  r = 0.5;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [0,0],  R = 1/4.
%
  c = [ 0.0, 0.0 ];
  r = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [0,0],  R = 1/8.
%
  c = [ 0.0, 0.0 ];
  r = 0.125;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );

  return
end
function circle_offcenter_test ( )

%*****************************************************************************80
%
%% CIRCLE_OFFCENTER_TEST uses a circle not centered at the origin.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2012
%
%  Author:
%
%    John Burkardt
%
  global c
  global r

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_OFFCENTER_TEST\n' );
%
%  C = [-1,0],  R = 1/2.
%
  c = [-1.0, 0.0 ];
  r = 0.50;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 2.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1), c(2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [0.5,0.5],  R = 1/4.
%
  c = [ 0.0, 0.0 ];
  r = 0.5;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [0.75,0.50],  R = 1/4.
%
  c = [ 0.75, 0.50 ];
  r = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );

  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  C = [1.0,1.0],  R = 1/2.
%
  c = [ 1.0, 1.0 ];
  r = 0.50;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r^2 / 4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R = %f, C = [%f,%f], Area = %f\n', r, c(1:2), area );
  shoreline ( @circle, m, n, x_min, x_max, y_min, y_max, step_num );

  return
end
function ellipse_centered_test ( )

%*****************************************************************************80
%
%% ELLIPSE_CENTERED_TEST uses an ellipse centered at the origin.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2012
%
%  Author:
%
%    John Burkardt
%
  global angle
  global c
  global r1
  global r2

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSE_CENTERED_TEST\n' );
%
%  Angle = 0.0, C = [0,0],  R1 = 0.75, R2 = 0.50;
%
  angle = 0.0;
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.50;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  Angle = 0.0, C = [0,0],  R1 = 0.75, R2 = 0.25
%
  angle = 0.0;
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  Angle = 0.0, C = [0,0],  R1 = 0.75, R2 = 0.125
%
  angle = 0.0;
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.125;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  Angle = 0.0, C = [0,0],  R1 = 0.75, R2 = 0.0625
%
  angle = 0.0;
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.0625;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );

  return
end
function ellipse_slanted_test ( )

%*****************************************************************************80
%
%% ELLIPSE_SLANTED_TEST uses a slanted ellipse centered at the origin.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2012
%
%  Author:
%
%    John Burkardt
%
  global angle
  global c
  global r1
  global r2

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSE_SLANTED_TEST\n' );
%
%  ANGLE = 0.0, C = [0,0],  R1 = 0.75, R2 = 0.25.
%
  angle = 0.0 * ( pi / 180.0 );
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  ANGLE = 30, C = [0,0],  R1 = 0.75, R2 = 0.25
%
  angle = 30.0 * ( pi / 180.0 );
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );
%
%  ANGLE = 45, C = [0,0],  R1 = 0.75, R2 = 0.125
%
  angle = 45.0 * ( pi / 180.0 );
  c = [ 0.0, 0.0 ];
  r1 = 0.75;
  r2 = 0.25;
  m = 11;
  n = 11;
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;
  step_num = 4;

  area = pi * r1 * r2;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Angle = %f, R = [%f,%f], C = [%f,%f], Area = %f\n', ...
    angle, r1, r2, c(1:2), area );
  shoreline ( @ellipse, m, n, x_min, x_max, y_min, y_max, step_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% LFCalDispEstPoses - Visualize pose estimates associated with a calibration info file
%
% Usage:
%     LFCalDispEstPoses( InputPath, CalOptions, DrawFrameSizeMult, BaseColour )
%
% Draws a set of frames, one for each camera pose, in the colour defined by BaseColour. All inputs
% except InputPath are optional. Pass an empty array "[]" to omit a parameter. See
% LFUtilCalLensletCam for example usage.
%
% Inputs:
%
%     InputPath : Path to folder containing CalInfo file.
%
%     [optional] CalOptions : struct controlling calibration parameters
%                 .CalInfoFname : Name of the file containing calibration estimate; note that this
%                                 parameter is automatically set in the CalOptions struct returned
%                                 by LFCalInit
%
%     [optional] DrawFrameSizeMult : Multiplier on the displayed frame sizes
%
%     [optional] BaseColour : RGB triplet defining a base drawing colour, RGB values are in the
%                             range 0-1
%
%
% User guide: LFToolbox.pdf
% See also:  LFUtilCalLensletCam

% Copyright (c) 2013-2020 Donald G. Dansereau

function LFCalDispEstPoses( InputPath, CalOptions, DrawFrameSizeMult, BaseColour )

%---Defaults---
CalOptions = LFDefaultField( 'CalOptions', 'CalInfoFname', 'CalInfo.json' );
DrawFrameSizeMult = LFDefaultVal( 'DrawFrameSizeMult', 1 );
BaseColour = LFDefaultVal( 'BaseColour', [1,1,0] );

%---
Fname = fullfile(CalOptions.CalInfoFname);
EstCamPosesV = LFStruct2Var( LFReadMetadata(fullfile(InputPath, Fname)), 'EstCamPosesV' );

% Establish an appropriate scale for the drawing
AutoDrawScale = EstCamPosesV(:,1:3);
AutoDrawScale = AutoDrawScale(:);
AutoDrawScale = (max(AutoDrawScale) - min(AutoDrawScale))/10;

%---Draw cameras---
BaseFrame = ([0 1 0 0 0 0; 0 0 0 1 0 0; 0 0 0 0 0 1]);
DrawFrameMed = [AutoDrawScale * DrawFrameSizeMult * BaseFrame; ones(1,6)];
DrawFrameLong = [AutoDrawScale*3/2 * DrawFrameSizeMult * BaseFrame; ones(1,6)];

for( PoseIdx = 1:size(EstCamPosesV, 1) )
    CurH = eye(4);
    CurH(1:3,1:3) = rodrigues( EstCamPosesV(PoseIdx,4:6) );
    CurH(1:3,4) = EstCamPosesV(PoseIdx,1:3);

    CurH = CurH^-1;

    CamFrame = CurH * DrawFrameMed(:,1:4);
    plot3( CamFrame(1,:), CamFrame(2,:), CamFrame(3,:), '-', 'color', BaseColour.*0.5, 'LineWidth',2 );
    hold on

    CamFrame = CurH * DrawFrameLong(:,5:6);
    plot3( CamFrame(1,:), CamFrame(2,:), CamFrame(3,:), '-', 'color', BaseColour, 'LineWidth',2 );
end

axis equal
grid on
xlabel('x [m]');
ylabel('y [m]');
zlabel('z [m]');
title('Estimated camera poses');

repo: LFToolbox


text:
classdef H1Projector_toP1Discontinuous < Projector

    methods (Access = public)

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

        function xProj = project(obj, x)
            LHS = obj.computeLHS();
            RHS = obj.computeRHS(x);
            f = LHS\RHS;
            fVals = obj.reshapeFValues(f, x.ndimf);
            s.mesh    = obj.mesh;
            s.fValues = fVals;
            xProj = P1DiscontinuousFunction(s);
        end

    end

    methods (Access = private)

        function LHS = computeLHS(obj)
            LHSM    = obj.computeMassMatrix();
            LHSK    = obj.computeStiffnessMatrix();
            epsilon = obj.mesh.computeMeanCellSize();
            LHS     = LHSM + epsilon^2*LHSK;
        end

        function LHSM = computeMassMatrix(obj)
            s.fun  = P1DiscontinuousFunction.create(obj.mesh, 1);
            s.mesh = obj.mesh;
            s.type = 'MassMatrix';
            s.quadratureOrder = 'QUADRATIC';
            lhs = LHSintegrator.create(s);
            LHSM = lhs.compute();
        end

        function LHSK = computeStiffnessMatrix(obj)
            s.type  = 'StiffnessMatrix';
            s.mesh  = obj.mesh;
            s.fun   = P1DiscontinuousFunction.create(obj.mesh,1);
            s.quadratureOrder = 'CONSTANT';
            lhs = LHSintegrator.create(s);
            LHSK = lhs.compute();
        end

        function RHS = computeRHS(obj,fun)
            quad = obj.createRHSQuadrature(fun);
            xV = quad.posgp;
            dV = obj.mesh.computeDvolume(quad);
            obj.mesh.interpolation.computeShapeDeriv(xV);
            shapes = permute(obj.mesh.interpolation.shape,[1 3 2]);
            conne = obj.createDiscontinuousConnectivity();

            nGaus = quad.ngaus;
            nFlds = fun.ndimf;
            nElem = obj.mesh.nelem;
            nNode = size(conne,2);
            nDofs = nElem*nNode;

            fGaus = fun.evaluate(xV);
            f     = zeros(nDofs,nFlds);
            for iField = 1:nFlds
                for igaus = 1:nGaus
                    dVg(:,1) = dV(igaus, :);
                    fG = squeeze(fGaus(iField,igaus,:));
                    for inode = 1:nNode
                        dofs = conne(:,inode);
                        Ni = shapes(inode,igaus);
                        int = Ni.*fG.*dVg;
                        f(dofs,iField) = f(dofs,iField) + int;
                    end
                end
            end
            RHS = f;
        end

        function conn = createDiscontinuousConnectivity(obj)
            nElem   = obj.mesh.nelem;
            nNodeEl = obj.mesh.nnodeElem;
            nDofs   = nElem*nNodeEl;
            dofs = 1:nDofs;
            conn = reshape(dofs, [nNodeEl,nElem])';
        end

        function q = createRHSQuadrature(obj, fun)
            ord = obj.determineQuadratureOrder(fun);
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature(ord);
        end

        function fVals = reshapeFValues(obj, x, nFlds)
            nElem = obj.mesh.nelem;
            nNode = obj.mesh.nnodeElem;
            fVals = reshape(x',nFlds, nNode, nElem);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%MAX2d	Maximum of image
%
%	[r,c] = max2d(image)
%
%	Return the interpolated coordinates (r,c) of the greatest peak in image.
%
% SEE ALSO:	ihough xyhough

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

function [r,c,m] = max2d(im)

	ncols = numcols(im);
	nrows = numrows(im);

	[mx,where] = max(im(:));

    [r,c] = ind2sub(size(im), where);

    m = mx;

	%[r,c,mx2]
	% now try to interpolate the peak over a 3x3 window

	% can't interpolate if against an edge
	if (c>1) & (c<ncols) & (r>1) & (r<nrows),
		dx = [
			c-1 c c+1
			c-1 c c+1
			c-1 c c+1];
		dy = [
			r-1 r-1 r-1
			r   r  r
			r+1   r+1  r+1];

		p = im(r-1:r+1,c-1:c+1);
		c = sum(sum(dx.*p)) / sum(sum(p));
		r = sum(sum(dy.*p)) / sum(sum(p));
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function op = prox_hinge( q , r, y)

%PROX_HINGE    Hinge-loss function.
%    OP = PROX_HINGE( q , r, y ) implements the nonsmooth function
%        OP(X) = q * sum( max( r - y.*x, 0 ) ).
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    then it must be a positive real scalar.
%    R is also optional; if omitted, R = 1 is assumed. R may be any real number.
%    Y is also optional; if omitted, Y = 1 is assumed. Y may be any scalar
%       or vector of the same size as X
% Dual: prox_hingeDual.m
%
% See also PROX_HINGEDUAL

if nargin < 3
    y = [];
elseif ~isempty(y) && ( ~isnumeric(y) || ~isreal(y) )
    error( 'Argument 3 must be a real vector');
end
if nargin < 2
	r = 1;
elseif ~isnumeric( r ) || ~isreal( r ) %|| numel( r ) ~= 1
	error( 'Argument 2 must be real.' );
end
if nargin < 1
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || numel( q ) ~= 1 || q <= 0,
	error( 'Argument 1 must be positive.' );
end

if isempty(y)
    op = tfocs_prox( @(x)hinge(x,r,q), @(x,t)prox_hinge(x,t,r,q) );
else
    ry  = r./y;
    op = tfocs_prox( @(x)hingeY(x,r,q,y), @(x,t)prox_hingeY(x,t,r,q,ry,y) );
end

% -- the actual functions --

function v = hingeY(x,r,q,y)
    if ~isscalar(r), assert( numel(r) == numel(x),'r is wrong size' ); end
    if ~isscalar(y), assert( numel(y) == numel(x),'y is wrong size'); end
    v = q*sum( max( r(:) - y(:).*x(:), 0 ) );
end
function v = hinge(x,r,q)
    if ~isscalar(r), assert( numel(r) == numel(x),'r is wrong size' ); end
    v = q*sum( max( r(:) - x(:), 0 ) );
end

%   PROX_F( Y, t ) = argmin_X  F(X) + 1/(2*t)*|| X - Y ||^2
function x = prox_hinge(x,t,r,q)
    tq = t * q;
    x = r + (x-r).*( x > r ) + (x + tq - r).*( x + tq < r  );
end
%{
    in the q = r = t = 1 case, the prox is:
prox(x) = {     x, if x > 1
                1, if x <= 1, and x > 0
                x + q, if x <= 0
%}

function x = prox_hingeY(x,t,r,q,ry,y)
    tq = t * q;
    x = ry + (x-ry).*( y.*x > r ) + (x + tq*y - ry).*( y.*(x + tq*y) < r );
end

end

% Added Feb, 2011; modified Dec, 2011

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS