text:
function [ss,gg,tt,ff,zo]=ssubmmse(si,fsz,pp)
%SSUBMMSE performs speech enhancement using mmse estimate of spectral amplitude or log amplitude [SS,ZO]=(S,FSZ,P)
%
% Usage: y=ssubmmse(x,fs);   % enhance the speech using default parameters
%
% Inputs:
%   si      input speech signal
%   fsz     sample frequency in Hz
%           Alternatively, the input state from a previous call (see below)
%   pp      algorithm parameters [optional]
%
% Outputs:
%   ss        output enhanced speech
%   gg(t,f,i) selected time-frequency values (see pp.tf below)
%   tt        centre of frames (in seconds)
%   ff        centre of frequency bins (in Hz)
%   zo        output state (or the 2nd argument if gg,tt,ff are omitted)
%
% The algorithm operation is controlled by a small number of parameters:
%
%        pp.of          % overlap factor = (fft length)/(frame increment) [2]
%        pp.ti          % desired frame increment [0.016 seconds]
%        pp.ri          % set to 1 to round ti to the nearest power of 2 samples [0]
%        pp.ta          % time const for smoothing SNR estimate [0.396 seconds]
%        pp.gx          % maximum posterior SNR as a power ratio [1000 = +30dB]
%        pp.gn          % min posterior SNR as a power ratio when estimating prior SNR [1 = 0dB]
%        pp.gz          % min posterior SNR as a power ratio [0.001 = -30dB]
%        pp.xn          % minimum prior SNR [0]
%        pp.xb          % bias compensation factor for prior SNR [1]
%        pp.lg          % MMSE target: 0=amplitude, 1=log amplitude [1]
%        pp.ne          % noise estimation: 0=min statistics, 1=MMSE [0]
%        pp.bt          % threshold for binary gain or -1 for continuous gain [-1]
%        pp.mx          % input mixture gain [0]
%        pp.rf          % round output signal to an exact number of frames [0]
%        pp.tf          % selects time-frequency planes to output in the gg() variable ['g']
%                           'i' = input power spectrum
%                           'I' = input complex spectrum
%                           'n' = noise power spectrum
%                           'z' = "posterior" SNR (i.e. (S+N)/N )
%                           'x' = "prior" SNR (i.e. S/N )
%                           'g' = gain
%                           'o' = output power spectrum
%                           'O' = output complex spectrum
%
% The applied gain is mx+(1-mx)*optgain where optgain is calculated according to [1] or [2].
% If pp.bt>=0 then optgain is first thresholded with pp.bt to produce a binary gain 0 or 1.
%
% The default parameters implement the original algorithm in [1,2].
%
% Several parameters relate to the estimation of xi, the so-called "prior SNR",
%
%             xi=max(a*pp.xb*xu+(1-a)*max(gami-1,pp.gn-1),pp.xn);
%
% This is estimated as a smoothed version of 1 less than gami, the "posterior SNR"
% which is the noisy speech power divided by the noise power. This is
% clipped to a min of (pp.gn-1), smoothed using a factor "a" which corresponds to a
% time-constant of pp.ta and then clipped to a minimum of pp.xn. The
% previous value is taken to be pp.xb*xu where xu is the ratio of the
% estimated speech amplitude squared to the noise power.
%
% In addition it is possible to specify parameters for the noise estimation algorithm
% which implements reference [3] or [7] according to the setting of pp.ne
%
% Minimum statistics noise estimate [3]: pp.ne=0
%        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]
%
% MMSE noise estimate [7]: pp.ne=1
%        pp.tax      % smoothing time constant for noise power estimate [0.0717 seconds](8)
%        pp.tap      % smoothing time constant for smoothed speech prob [0.152 seconds](23)
%        pp.psthr    % threshold for smoothed speech probability [0.99] (24)
%        pp.pnsaf    % noise probability safety value [0.01] (24)
%        pp.pspri    % prior speech probability [0.5] (18)
%        pp.asnr     % active SNR in dB [15] (18)
%        pp.psini    % initial speech probability [0.5] (23)
%        pp.tavini   % assumed speech absent time at start [0.064 seconds]
%
% If convenient, you can call specsub in chunks of arbitrary size. Thus the following are equivalent:
%
%                   (a) y=ssubmmse(s,fs);
%
%                   (b) [y1,z]=ssubmmse(s(1:1000),fs);
%                       [y2,z]=ssubmmse(s(1001:2000),z);
%                       y3=ssubmmse(s(2001:end),z);
%                       y=[y1; y2; y3];
%
% If the number of output arguments is either 2 or 5, the last partial frame of samples will
% be retained for overlap adding with the output from the next call to ssubmmse().
%
% See also specsub() for an alternative gain function
%
% Refs:
%    [1] Ephraim, Y. & Malah, D.
%        Speech enhancement using a minimum-mean square error short-time spectral amplitude estimator
%        IEEE Trans Acoustics Speech and Signal Processing, 32(6):1109-1121, Dec 1984
%    [2] Ephraim, Y. & Malah, D.
%        Speech enhancement using a minimum mean-square error log-spectral amplitude estimator
%        IEEE Trans Acoustics Speech and Signal Processing, 33(2):443-445, Apr 1985
%    [3] 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.
%    [4] O. Cappe.
%        Elimination of the musical noise phenomenon with the ephraim and malah noise suppressor.
%        IEEE Trans Speech Audio Processing, 2 (2): 345349, Apr. 1994. doi: 10.1109/89.279283.
%    [5] J. Erkelens, J. Jensen, and R. Heusdens.
%        A data-driven approach to optimizing spectral speech enhancement methods for various error criteria.
%        Speech Communication, 49: 530541, 2007. doi: 10.1016/j.specom.2006.06.012.
%    [6] R. Martin.
%        Statistical methods for the enhancement of noisy speech.
%        In J. Benesty, S. Makino, and J. Chen, editors,
%        Speech Enhancement, chapter 3, pages 4364. Springer-Verlag, 2005.
%    [7] Gerkmann, T. & Hendriks, R. C.
%        Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay
%        IEEE Trans Audio, Speech, Language Processing, 2012, 20, 1383-1393

% Bugs/suggestions:
%   (1) sort out behaviour when si() is a matrix rather than a vector
%
%      Copyright (C) Mike Brookes 2004-2011
%      Version: $Id: ssubmmse.m 2460 2012-10-29 22:20:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if numel(si)>length(si)
    error('Input speech signal must be a vector not a matrix');
end
if isstruct(fsz)
    fs=fsz.fs;
    qq=fsz.qq;
    qp=fsz.qp;
    ze=fsz.ze;
    s=zeros(length(fsz.si)+length(si(:)),1); % allocate space for speech
    s(1:length(fsz.si))=fsz.si;
    s(length(fsz.si)+1:end)=si(:);
else
    fs=fsz;     % sample frequency
    s=si(:);
    % default algorithm constants

    qq.of=2;        % overlap factor = (fft length)/(frame increment)
    qq.ti=16e-3;    % desired frame increment (16 ms)
    qq.ri=0;        % round ni to the nearest power of 2
    qq.ta=0.396;    % Time const for smoothing SNR estimate = -tinc/log(0.98) from [1]
    qq.gx=1000;     % maximum posterior SNR = 30dB
    qq.gn=1;        % min posterior SNR as a power ratio when estimating prior SNR [1]
    qq.gz=0.001;    % min posterior SNR as a power ratio [0.001 = -30dB]
    qq.xn=0;        % minimum prior SNR = -Inf dB
    qq.xb=1;        % bias compensation factor for prior SNR [1]
    qq.lg=1;        % use log-domain estimator by default
    qq.ne=0;        % noise estimation: 0=min statistics, 1=MMSE [0]
    qq.bt=-1;       % suppress binary masking
    qq.mx=0;        % no input mixing
    qq.tf='g';      % output the gain time-frequency plane by default
    qq.rf=0;
    if nargin>=3 && ~isempty(pp)
        qp=pp;      % save for estnoisem call
        qqn=fieldnames(qq);
        for i=1:length(qqn)
            if isfield(pp,qqn{i})
                qq.(qqn{i})=pp.(qqn{i});
            end
        end
    else
        qp=struct;  % make an empty structure
    end
end
% derived algorithm constants
if qq.ri
    ni=pow2(nextpow2(ti*fs*sqrt(0.5)));
else
    ni=round(qq.ti*fs);    % frame increment in samples
end
tinc=ni/fs;         % true frame increment time
a=exp(-tinc/qq.ta); % SNR smoothing coefficient
gx=qq.gx;           % max posterior SNR as a power ratio
gz=qq.gz;           % min posterior SNR as a power ratio
kk=sqrt(2*pi);      % sqrt(8)*Gamma(1.5) - required constant
xn=qq.xn;           % floor for prior SNR, xi
ne=qq.ne;           % noise estimation: 0=min statistics, 1=MMSE [0]
gn1=max(qq.gn-1,0); % floor for posterior SNR when estimating prior SNR
xb=qq.xb;
tf=qq.tf;
rf=qq.rf || nargout==2 || nargout==5;            % round down to an exact number of frames

% calculate power spectrum in frames

no=round(qq.of);                  	% integer overlap factor
nf=ni*no;                           % fft length
w=sqrt(hamming(nf+1))'; w(end)=[];  % for now always use sqrt hamming window
w=w/sqrt(sum(w(1:ni:nf).^2));       % normalize to give overall gain of 1
if rf>0
    rfm='';                         % truncated input to an exact number of frames
else
    rfm='r';
end
[y,tt]=enframe(s,w,ni,rfm);
tt=tt/fs;                           % frame times
yf=rfft(y,nf,2);
yp=yf.*conj(yf);                    % power spectrum of input speech
[nr,nf2]=size(yp);                  % number of frames
ff=(0:nf2-1)*fs/nf;
if isstruct(fsz)
    if ne>0
        [dp,ze]=estnoiseg(yp,ze);       % estimate the noise using MMSE
    else
        [dp,ze]=estnoisem(yp,ze);       % estimate the noise using minimum statistics
    end
    ssv=fsz.ssv;
    xu=fsz.xu;                      % saved unsmoothed SNR
else
    if ne>0
        [dp,ze]=estnoiseg(yp,tinc,qp);	% estimate the noise using MMSE
    else
        [dp,ze]=estnoisem(yp,tinc,qp);	% estimate the noise using minimum statistics
    end
    ssv=zeros(ni*(no-1),1);        	% dummy saved overlap
    xu=1;                           % dummy unsmoothed SNR from previous frame
end
if ~nr                             	% no data frames
    ss=[];
    gg=[];
else
    gam=max(min(yp./dp,gx),gz);     % gamma = posterior SNR
    g=zeros(nr,nf2);                % create space for gain matrix
    x=zeros(nr,nf2);                % create space for prior SNR
    if qq.lg                        % use log domain estimator
        for i=1:nr
            gami=gam(i,:);
            xi=max(a*xb*xu+(1-a)*max(gami-1,gn1),xn);  % prior SNR
            xir=xi./(1+xi);
            gi=xir.*exp(0.5*expint(xir.*gami));
            g(i,:)=gi;             	% save gain for later
            x(i,:)=xi;              % save prior SNR
            xu=gami.*gi.^2;         % unsmoothed prior SNR
        end
    else
        for i=1:nr
            gami=gam(i,:);
            xi=max(a*xb*xu+(1-a)*max(gami-1,gn1),xn);  % prior SNR
            v=0.5*xi.*gami./(1+xi);	% note that this is 0.5*vk in [1]
            gi=(0.277+2*v)./gami; 	% accurate to 0.02 dB for v>0.5
            mv=v<0.5;
            if any(mv)
                vmv=v(mv);
                gi(mv)=kk*sqrt(vmv).*((0.5+vmv).*besseli(0,vmv)+vmv.*besseli(1,vmv))./(gami(mv).*exp(vmv));
            end
            g(i,:)=gi;              % save gain for later
            x(i,:)=xi;              % save prior SNR
            xu=gami.*gi.^2;         % unsmoothed prior SNR
        end
    end
    if qq.bt>=0
        g=g>qq.bt;
    end
    g=qq.mx+(1-qq.mx)*g;                    % mix in some of the input
    se=(irfft((yf.*g).',nf).').*repmat(w,nr,1);     % inverse dft and apply output window
    ss=zeros(ni*(nr+no-1),no);                      % space for overlapped output speech
    ss(1:ni*(no-1),end)=ssv;
    for i=1:no
        nm=nf*(1+floor((nr-i)/no));         % number of samples in this set
        ss(1+(i-1)*ni:nm+(i-1)*ni,i)=reshape(se(i:no:nr,:)',nm,1);
    end
    ss=sum(ss,2);
    if nargout>2 && ~isempty(tf)
        gg=zeros(nr,nf2,length(tf));  % make space
        for i=1:length(tf)
            switch tf(i)
                case 'i'            % 'i' = input power spectrum
                    gg(:,:,i)=yp;
                case 'I'            % 'i' = input power spectrum
                    gg(:,:,i)=yf;
                case 'n'            % 'n' = noise power spectrum
                    gg(:,:,i)=dp;
                case 'z'            % 'z' = posterior SNR (i.e. (S+N)/N )
                    gg(:,:,i)=gam;
                case 'x'            % 'x' = prior SNR
                    gg(:,:,i)=x;
                case 'g'            % 'g' = gain
                    gg(:,:,i)=g;
                case 'o'            % 'o' = output power spectrum
                    gg(:,:,i)=yp.*g.^2;
                case 'O'            % 'o' = output power spectrum
                    gg(:,:,i)=yf.*g;
            end
        end
    end
end
if nargout==2 || nargout==5
    if nr
        zo.ssv=ss(end-ni*(no-1)+1:end);    % save the output tail for next time
        ss(end-ni*(no-1)+1:end)=[];        % only output the frames that are completed
    else
        zo.ssv=ssv;  %
    end
    zo.si=s(length(ss)+1:end);      % save the tail end of the input speech signal
    zo.fs=fs;                       % save sample frequency
    zo.qq=qq;                       % save local parameters
    zo.qp=qp;                       % save estnoisem parameters
    zo.ze=ze;                       % save state of noise estimation
    zo.xu=xu;
    if nargout==2
        gg=zo;                      % 2nd of two arguments is zo
    end
elseif rf==0
    ss=ss(1:length(s));             % trim to the correct length if not an exact number of frames
end
if ~nargout && nr>0
    ffax=ff/1000;
    ax=zeros(4,1);
    ax(1)=subplot(223);
    imagesc(tt,ffax,20*log10(g)');
    colorbar;
    axis('xy');
    title(sprintf('Filter Gain (dB): ta=%.2g',qq.ta));
    xlabel('Time (s)');
    ylabel('Frequency (kHz)');

    ax(2)=subplot(222);
    imagesc(tt,ffax,10*log10(yp)');
    colorbar;
    axis('xy');
    title('Noisy Speech (dB)');
    xlabel('Time (s)');
    ylabel('Frequency (kHz)');

    ax(3)=subplot(224);
    imagesc(tt,ffax,10*log10(yp.*g.^2)');
    colorbar;
    axis('xy');
    title('Enhanced Speech (dB)');
    xlabel('Time (s)');
    ylabel('Frequency (kHz)');

    ax(4)=subplot(221);
    imagesc(tt,ffax,10*log10(dp)');
    colorbar;
    axis('xy');
    title('Noise Estimate (dB)');
    xlabel('Time (s)');
    ylabel('Frequency (kHz)');
    linkaxes(ax);
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
 function [err,f1,f2] = nufft2_err_mm(om, N1, N2, J1, J2, K1, K2, alpha, beta)
%function [err,f1,f2] = nufft2_err_mm(om, N1, N2, J1, J2, K1, K2, alpha, beta)
% Compute worst-case error for each input frequency for min-max 2D NUFFT.
% in:
%	om	[M,2]	digital frequency omega in radians
%	N1,N2		signal length
%	J1,J2		# of neighbors used per frequency location
%	K1,K2		FFT size (should be > N1)
%	alpha	[L,1]	Fourier series coefficients of scaling factors
%	beta		scale gamma=2pi/K by this in Fourier series
%				typically is K/N (me) or 0.5 (Liu)
% out:
%	err	[M,1]	worst-case error over unit-norm signals
%
% Copyright 2001-12-7, Jeff Fessler, The University of Michigan

% if no arguments, give an example
if nargin < 4
	help(mfilename)
	N1 = 1; K1 = 2*N1; J1 = 7; gam1 = 2*pi/K1;
	N2 = 1; K2 = 2*N2; J2 = 6; gam2 = 2*pi/K2;
	alpha = 1;
	alpha = 'best';
	[err,f1,f2] = nufft2_err_mm('all', N1, N2, J1, J2, K1, K2, alpha);
	mesh(f1, f2, err)
	plot(f1, err)
	xlabel '\omega_1 / \gamma', ylabel 'E_{max}(\omega_1,\omega_2)'
	clear err
return
end

if ~isvar('alpha') | isempty(alpha)
	alpha = [1];	% default Fourier series coefficients of scaling factors
end
if ~isvar('beta') | isempty(beta)
	beta = 0.5;	% default is Liu version for now
end
alpha1 = alpha;
alpha2 = alpha;
beta1 = beta;
beta2 = beta;

%
% trick to look at all relevant om's
%
if ischar(om)
	if ~streq(om, 'all'), error 'unknown om argument', end
	gam1 = 2*pi/K1;
	gam2 = 2*pi/K2;
	[f1,f2] = ndgrid(linspace(0,1,31), linspace(0,1,33));
	om = [gam1*f1(:) gam2*f2(:)];
end

%
% see if 'best' alpha is desired
%
if ischar(alpha)
	if ~streq(alpha, 'best'), error 'unknown alpha argument', end
	[alpha1, beta1, ok] = nufft_best_alpha(J1, 0, K1/N1);
	if ~ok, error 'unknown J,K/N', end
	[alpha2, beta2, ok] = nufft_best_alpha(J2, 0, K2/N2);
	if ~ok, error 'unknown J,K/N', end
end

tol = 0;
T1 = nufft_T(N1, J1, K1, tol, alpha1, beta1);		% [J,J]
T2 = nufft_T(N2, J2, K2, tol, alpha2, beta2);		% [J,J]
r1 = nufft_r(om(:,1), N1, J1, K1, alpha1, beta1);	% [J,M]
r2 = nufft_r(om(:,2), N2, J2, K2, alpha2, beta2);	% [J,M]
T = kron(T2, T1);	% J1*J2 x J1*J2

% kron for each om
M = size(om,1);
r = zeros(J1*J2,M);
for ii=1:M
	r(:,ii) = kron(r2(:,ii), r1(:,ii));
end

%
% worst-case error at each frequency
%
Tr = T * r;			% [J,M]
err = sum(conj(r) .* Tr).';	% [M,1]
err = reale(err);
err = min(err, 1);
err = sqrt(1 - err);

if isvar('f1')
	err = reshape(err, size(f1));
end

author: thomaskuestner


text:
function pth = sessionGUI_roiPath(roi);
% return a full path to an ROI file, given the name of an ROI in the
% current ROI list, or the numeric index. If a full path is already
% provided, return it unchanged.
mrGlobals2;

% get the parent directory for ROIs
switch get(GUI.controls.roiType, 'Value')
	case 1, parent = roiDir(INPLANE{1},1); % inplane
	case 2, parent = roiDir(VOLUME{1},0); % volume shared
	case 3, parent = roiDir(VOLUME{1},1); % gray local
	case 4, parent = fullfile(HOMEDIR, 'Volume', 'ROIs');  % volume local
	case 5,  % browse for directory
		parent = get(GUI.controls.roiType, 'String');
		parent = parent{end};
end

% use the roi value and the parent directory to construct a full ROI path.
if ischar(roi) & exist(roi, 'file')
    pth = fullpath(roi);

elseif isnumeric(roi)
    roiNames = get(GUI.controls.roi, 'String');
    pth = fullfile(parent, roiNames{roi(1)});

elseif ischar(roi)
    pth = fullfile(parent, roi);

else
    error('Invalid ROI specification. ')

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function S = symmetrize(M)

% SYMMETRIZE Make a matrix symmetric.
%   SYMMETRIZE(M) for square matrix M, returns a symmetric version of M by
%   doing (M+M')/2. It issues an error for non-square matrices.

S = (M+M')/2;

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function SO3VF = exp(F)
% overloads |exp(SO3VF)|
%
% Syntax
%   SO3VF = exp(SO3VF)
%
% Input
%  SO3VF - @SO3VectorField
%
% Output
%  SO3VF - @SO3VectorField
%

SO3VF = SO3VectorFieldHandle(@(rot) g(rot),F.SRight,F.SLeft);

function g = g(rot)
g = F.eval(rot);
g = vector3d(exp(g.xyz.'));
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [y,exitCode]=geometricMedian(x,eta,AbsTol,RelTol,maxIter,test4Duplicates)
%%GEOMETRICMEDIAN Compute the weighted geometric median of the (possibly
%          multivariate) points x. This finds the vector y that minimizes
%          sum_i eta(i)*norm(y-x(:,i),1) where all eta(i)>0.
%          The geometric median is also known as the L1-median and the
%          spatial median.
%
%INPUTS: x The xDimXN set of N real vectors whose geometric median is
%          desired.
%      eta The NX1 or 1XN set of all positive weights. If this parameter is
%          omitted or an empty matrix is passed, then uniform weights are
%          used.
% AbsTol, RelTol Absolute and relative tolerances for determining
%          convergence. If diff is the difference between the estimate at
%          the current and the previous iteration, then convergence is
%          declared if diff<AbsTol||diff<RelTol*norm(y). These same
%          tolerances are used in a similar manner to determine whether y
%          equals any of the points in x during an iterative update and
%          whether any of the vectors in x are duplicates of each other.
%          The defaults if omitted or empty matrices are passed are
%          RelTol=1e-9, AbsTol=1e-12.
%  maxIter The maximum number of iterations to perform. The default if
%          ommitted or an empty matrix is passed is max(1000,20*xDim).
% test4Duplicates The algorithms assume that none of the vectors in x are
%          duplicated of each other. test4Duplicates can be set to false if
%          it is known that there are no duplicated vectors in x. The
%          default if omitted or an empty matrix is passed is true.
%
%OUTPUTS: y The xDimX1 geometric median of the vectors in x.
%  exitCode A value indicating how the iterative algorithm terminated.
%           Possible values are:
%           0 The algorithm converged.
%           1 The algorithm did not convergence withing maxiter iteration.
%
%This function uses the modified Weiszfeld algorithm, which is described in
%[1]. The algorithm is inistialized with the Convergence is linear. A
%quasi-Newton algorithm would presumably be better for very large problems.
%
%EXAMPLE:
%Here, we have a set of data with an outlier. one will see that the mean is
%pulled towards the outlier, whereas the geometric median is more robust.
% x=[-26,   5,   1,  -2,  13,  -1,   0,   3,  -2, -12,  120;
%      0,  12,  12,  19,  -8, -16,  -5,   3,   4,   6,  100];
% yGM=geometricMedian(x);
% yM=mean(x,2);
% costGM=sum(sqrt(sum((bsxfun(@minus,yGM,x)).^2,1)))
% costM=sum(sqrt(sum((bsxfun(@minus,yM,x)).^2,1)))
% figure(1)
% clf
% hold on
% scatter(x(1,:),x(2,:),'ok')
% scatter(yGM(1),yGM(2),'r','filled')
% scatter(yM(1),yM(2),'b','filled')
% legend('Points','Geometric Median','Mean')
% h1=xlabel('x');
% h2=ylabel('y');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] Y. Vardi and C.-H. Zhang, "The multivariate l1-median and associated
%    data depth," Proceedings of the National Academy of Sciences of the
%    United States of America, vol. 97, no. 4, pp. 1423-1426, 23 Feb. 2000.
%
%August 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(x,2);
xDim=size(x,1);

if(nargin<6||isempty(test4Duplicates))
    test4Duplicates=true;
end

if(nargin<5||isempty(maxIter))
    maxIter=max(1000,20*xDim);
end

if(nargin<4||isempty(RelTol))
    RelTol=1e-9;
end

if(nargin<3||isempty(AbsTol))
    AbsTol=1e-12;
end

if(nargin<2||isempty(eta))
    eta=ones(1,n);
end

%Get rid of duplicate values and adjust the weights.
if(test4Duplicates)
    xUnique=zeros(xDim,n);
    xUnique(:,1)=x(:,1);
    etaUnique=zeros(1,n);
    etaUnique(1)=eta(1);
    numUnique=1;

    xMag=sqrt(sum(x.*x,1));
    for curX=2:n

        d=x(:,curX)-xUnique(:,1:numUnique);
        dMag=sqrt(sum(d.^2,1));
        idx=find(dMag<AbsTol|dMag<RelTol*xMag(curX),1);

        if(isempty(idx))%It is unique.
          numUnique=numUnique+1;
          xUnique(:,numUnique)=x(:,curX);
          etaUnique(numUnique)=eta(curX);
        else%It is not unique.
            etaUnique(idx)=etaUnique(idx)+eta(curX);
        end
    end

    x=xUnique;
    eta=etaUnique;
    n=numUnique;
end

%The special case when given a single value.
if(n==1)
  y=x;
  exitCode=0;
  return
end

eta=eta(:).';
xMag=sqrt(sum(x.*x,1));%Magnitude

%Start with the weighted mean as the initial estimate.
y=sum(bsxfun(@times,(eta/sum(eta)),x),2);

exitCode=1;
for curIter=1:maxIter
  yPrev=y;

  %Duplicates were removed (weights combined) prior to beginning iterations, so
  %there should be at most 1 point coincident with x. We will use absolute and
  %relative tolerances to determine equality.
  d=bsxfun(@minus,y,x);
  dMag=sqrt(sum(d.^2,1));
  idx=find(dMag<AbsTol|dMag<RelTol*xMag,1);

  if(isempty(idx))
    %If no points coincide with x, then use Weiszfeld's original iteration.

    %Defined after Equation 2.3
    w=eta./dMag;
    w=w/sum(w);

    %Equation 2,3
    y=sum(bsxfun(@times,w,x),2);
  else
    %If any point in x coincides with y, then use the modified Weiszfeld
    %iteration.

    %Equation 2,5. The only nonzero weight is the one  where the values
    %coincide.
    etay=eta(idx);

    %w as in Equation 23, but not normalized.
    w=eta./dMag;
    w(idx)=0;
    %The normalization constant.
    wNorm=sum(w);

    %Equation 2,4
    Ty=sum(bsxfun(@times,(w/wNorm),x),2);

    %Equation 2,7
    ry=norm(sum(bsxfun(@times,w,-d),2),1);

    %Equation 2.6
    y=max(0,(1-etay/ry))*Ty+min(1,etay/ry)*y;
  end

  yNorm=norm(y);
  diff=norm(y-yPrev);

  if(diff<AbsTol||diff<RelTol*yNorm)
      exitCode=0;
      break;
  end
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
 function yy = ir_imfill1(xx)
%function yy = ir_imfill1(xx)
%|
%| 1D version of imfill 'holes' that works along 1st dimension of input
%|
%| in
%| xx [N (L)]	logical array
%| out
%| yy [N (L)]	logical array
%|
%| 2015-08-13, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(xx, 'test'), ir_imfill1_test, return, end

if ~islogical(xx), fail('input must be logical'), end

yy = cummax(xx,1) & flipdim(cummax(flipdim(xx,1),1),1);

function ir_imfill1_test
x = ellipse_im(64) > 0;
x(end/2,end/4+[1:5]) = 0;
y1 = ir_imfill1(x);
y2 = ir_imfill1(x')';
im plc 2 2
im(1, x)
im(2, y1)
im(3, y2)
im(4, y2-y1)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef OrthogonalLeastSquares < handle
    %CSOLSAPPROX Implements OLS algorithm for sparse approximation

    properties
        % Maximum residual norm
        MaxResNorm = 1e-4
        % Indicates if we should stop on exceeding residual norm
        StopOnResidualNorm = true
        % Indicates if we should stop when residual norm stops improving
        StopOnResNormStable = true
        % Maximum number of iterations for approximation
        MaxIters
        Verbose = false
        % Minimum Sparsity
        MinK = 4
        % Ignored atom (which won't be considered in identification step)
        IgnoredAtom = -1
    end

    properties(SetAccess=private)
        % The dictionary
        Dict
        % Ambient signal dimensions
        N
        % Number of atoms in dictionary
        D
        % Sparsity level of representations (may be negative)
        K
        % Result of a solver
        result
    end

    methods
        function self  = OrthogonalLeastSquares(Dict, K)
            % We assume that all the columns in dictionary are normalized.
            if isa(Dict, 'spx.dict.Operator')
                self.Dict = Dict;
            elseif ismatrix(Dict)
                self.Dict = spx.dict.MatrixOperator(Dict);
            else
                error('Unsupported operator.');
            end
            [self.N, self.D] = size(Dict);
            if ~exist('K', 'var')
                % No sparsity level has been pre-specified.
                K = -1;
            end
            self.K = K;
            % Maximum number of iterations
            maxIter = self.N;
            if K > 0
                % We have to consider pre-specified sparsity level
                maxIter = K;
            end
            self.MaxIters = maxIter;
            if K > 0 && self.MinK >= K
                self.MinK = 0;
            end
        end

        function result  = solve(self,y)
            % Uses QR factorization to implement the OLS

            % Initialization
            % Solves approximation problem using OLS
            d = self.D;
            n = self.N;
            r = y;
            min_k  = self.MinK;
            dict = self.Dict;
            % Active indices
            omega = [];
            if self.StopOnResNormStable
                oldResNorm = norm(r);
            end
            maxIter = self.MaxIters;
            % Create space for storing the Q R factors
            Q = zeros(n, maxIter);
            R = zeros(maxIter);
            % zr stores projections of y along the directions
            % in Q.
            % Since vectors in Q are orthonormal,
            % hence computing projections is as easy as
            % taking inner product.
            % This also simplifies the process of updating
            % the residual as the residual is orthogonal
            % to the vectors selected in Q so far.
            zr = [];

            ignored_atom = self.IgnoredAtom;
            % Convert the dictionary into a matrix
            dict = double(dict);

            for iter=1:maxIter
                % The Q part of previous iteration
                QLast = Q(:, 1:iter-1); % n * (k - 1)
                % The Orthogonal Projector
                OrthoProjector = eye(n) -  QLast * QLast';
                % Orthogonalize the remaining atoms of the dictionary
                qdict = OrthoProjector * dict;
                % Normalize the columns
                qdict = spx.norm.normalize_l2(qdict);
                % Compute inner products
                innerProducts = qdict' * r;
                % Mark the inner products of already selected columns as 0.
                innerProducts(omega) = 0;
                innerProducts = abs(innerProducts);
                if ignored_atom > 0
                    % forcefully ignore this atom
                    innerProducts(ignored_atom) = 0;
                end
                % Find the highest inner product
                [~, index] = max(innerProducts);
                % Add this index to support
                omega = [omega, index];
                % Pick the new atom from the dictionary
                new_atom = dict(:, index);
                % Orthogonalize the new atom
                %Compute projections to previously selected vectors in Q
                projections = QLast'* new_atom;
                % Remove the projection
                new_q = new_atom - QLast * projections;
                % Normalize
                norm_q = norm(new_q);
                new_q = new_q / norm_q;
                % Place it
                Q(:, iter) = new_q;
                % Update R
                R(1:iter-1, iter) = projections;
                R(iter, iter) = norm_q;
                % Compute the projection of y on new q.
                zr(iter) = new_q' * y;
                % Let us update the residual.
                r = r - zr(iter) * new_q;
                if self.StopOnResidualNorm || self.StopOnResNormStable
                    resNorm = norm(r);
                    if resNorm < self.MaxResNorm && iter > min_k
                        break;
                    end
                    if self.StopOnResNormStable
                        change = abs(oldResNorm  - resNorm);
                        if change/oldResNorm < .01  && iter > min_k
                            % No improvement
                            break;
                        end
                    end
                end
            end
            % Estimate
            z = zeros(d, 1);
            % We need to use back-substitution  with R to get z from zr.
            opts.UT = true;
            tmp = linsolve(R(1:iter,1:iter), zr(1:iter)', opts);
            z(omega)= tmp;
            % Solution vector
            result.z = z;
            % Residual obtained
            result.r = r;
            % residual norm
            result.rnorm = resNorm;
            % Number of iterations
            result.iterations = iter;
            % Solution support
            result.support = omega;
            self.result = result;
        end

    end

end

author: indigits

repo: sparse-plex


text:
% Features computation for the Harmonic Model + Phase Distortion (HMPD)
%
% Octave compatible
%
% Please read the README.txt file for general information about HMPD before
% using it.
%
% Inputs
%  frames : [Nxstruct] N structures containing the sinusoidal parameters,
%           as provided by hmpd_analysis_harmonic.m or another sinusoidal
%           analysis.
%  fs     : [Hz] The sampling rate of the analyzed waveform
%  [opt]  : Additional options (see code below)
%           Even though this argument is optional, it is highly recommended to
%           adapt the f0min and f0max options if no f0s argument is provided.
%
% Outputs
%  f0s   : [s, Hz] [Nx2] A time/data column vector, as above, containing the
%          analysis instants and the f0 curve.
%  AE    : [NxD] A matrix containing the amplitude envelope.
%          D is either opt.dftlen/2+1 (from hmpd_features_compute.m), or
%          opt.amp_order+1, depending if compression is disabled or enabled.
%  PDM   : [NxD] A matrix containing the Phase Distortion Mean.
%          D is either opt.dftlen/2+1 or opt.pdm_order+1, depending if
%          compression is disabled or enabled.
%  PDD   : [NxD] A matrix containing the Phase Distortion Deviation.
%          D is either opt.dftlen/2+1 or opt.pdd_order+1, depending if
%          compression is disabled or enabled.
%  opt   : The options used for the analysis.
%
% Copyright (c) 2013 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [f0s, AE, PDM, PDD, opt] = hmpd_analysis_features(frames, fs, opt)
%  keyboard
    if opt.debug>0; disp('HMPD Vocoder: Features computation ...'); end

    irregf0s = [[frames.t]', [frames.f0]']; % Retrieve f0 from the frames

    if opt.debug>0; disp('    Estimate amplitude envelope ...'); end
    [AE, frames] = hmpd_amplitude_envelope_estimate(frames, fs, opt);

    if opt.debug>0; disp('    Estimate phase envelope ...'); end
    rpspdopt = opt;
    rpspdopt.pd_method = 1;
    rpspdopt.pd_vtf_rm =false;% Already done in hmpd_amplitude_envelope_estimate
    rpspdopt.harm2freq = true;
    rpspdopt.usemex = opt.usemex;
    PE = phase_rpspd(frames, fs, rpspdopt);

    if opt.debug>0; disp('    Compute phase statistics ...'); end

    % Phase's mean
    PDM = hmpd_phase_mean(PE, opt.pdm_nbper*opt.sin_nbat); % Below 3, doubles the voice, why?

    % Phase's standard-deviation
    % Remove first the trend from PE, otherwise the std measure is overestimated
    PEtrend = hmpd_phase_smooth(PE, opt.pdd_nbper*opt.sin_nbat); % Compute the trend
    PE = PE - PEtrend;

    PDD = hmpd_phase_deviation(PE, opt.pdd_nbper*opt.sin_nbat);

    if opt.debug>0; disp('    Compress features ...'); end
    [AE, PDM, PDD] = hmpd_features_compress(irregf0s, AE, PDM, PDD, fs, opt);

    if opt.debug>0; disp('    Resample features with uniform sampling ...'); end
    uT = (0:opt.regularstepsize:irregf0s(end,1))'; % New uniform time instants

    % First, resample the f0 curve using uniform step size
    uf0s = exp(interp1(irregf0s(:,1), log(irregf0s(:,2)), uT, 'linear', NaN));

    % Drop the time instants outside of the features
    idx = find(~isnan(uf0s));
    uT = uT(idx);
    f0s = [uT, uf0s(idx)];

    % Then, resample amplitude, PD's mean and PD's deviation
    AE = irregsampling2uniformsampling(irregf0s(:,1), AE, uT, [], [], 'linear', NaN, opt.usemex);
    PDM = irregsampling2uniformsampling(irregf0s(:,1), PDM, uT, @unwrap, @wrap, 'linear', 0, opt.usemex);
    PDD = irregsampling2uniformsampling(irregf0s(:,1), PDD, uT, [], [], 'linear', 0, opt.usemex);

return


text:
function area = sphere01_triangle_vertices_to_area ( v1, v2, v3 )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_VERTICES_TO_AREA: area of a spherical triangle on unit sphere.
%
%  Discussion:
%
%    A unit sphere centered at 0 in 3D satisfies the equation:
%
%      X * X + Y * Y + Z * Z = 1
%
%    A spherical triangle is specified by three points on the surface
%    of the sphere.
%
%    The area formula is known as Girard's formula.
%
%    The area of a spherical triangle on the unit sphere is:
%
%      AREA = ( A + B + C - PI )
%
%    where A, B and C are the (surface) angles of the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real V1(3), V2(3), V3(3), the vertices of the triangle.
%
%    Output, real AREA, the area of the spherical triangle.
%

%
%  Compute the lengths of the sides of the spherical triangle.
%
  [ as, bs, cs ] = sphere01_triangle_vertices_to_sides ( v1, v2, v3 );
%
%  Get the spherical angles.
%
  [ a, b, c ] = sphere01_triangle_sides_to_angles ( as, bs, cs );
%
%  Get the area.
%
  area = sphere01_triangle_angles_to_area ( a, b, c );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Y1 h11]=derotated_dtcwt(I1,n,biot,Qshift)
%     X -> 2D real matrix/Image
%
%     nlevels -> No. of levels of wavelet decomposition
%
%     biort ->  'antonini'   => Antonini 9,7 tap filters.
%               'legall'     => LeGall 5,3 tap filters.
%               'near_sym_a' => Near-Symmetric 5,7 tap filters.
%               'near_sym_b' => Near-Symmetric 13,19 tap filters.
%
%     qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,
%                              (only 6,6 non-zero taps).
%               'qshift_a' =>  Q-shift 10,10 tap filters,
%                              (with 10,10 non-zero taps, unlike qshift_06).
%               'qshift_b' => Q-Shift 14,14 tap filters.
%               'qshift_c' => Q-Shift 16,16 tap filters.
%               'qshift_d' => Q-Shift 18,18 tap filters.
%
%
%     Y1     -> The real lowpass image from the final level
%     h11     -> A cell array containing the 6 complex highpass subimages
%     for each level.
%clear all;
%clc;
[Y1,h1] = dtwavexfm2(I1,n,biot,Qshift);
%[Y2,h2] = dtwavexfm2(I2,n,biot,Qshift);

h11{n}=h1{n};
for k=n:-1:2
    for m=1:6
        xp=imresize(h1{k}(:,:,m),2);%
        argxp=angle(xp);
        argx=angle(h1{k-1}(:,:,m));
        argx=argx-2.*argxp;
        absx=abs(h1{k-1}(:,:,m));
        xa=absx.*cos(argx);
        xb=absx.*sin(argx);
        h11{k-1}(:,:,m)=complex(xa,xb);
    end
end
%figure;
%cimage5(h11{1}(:,:,4));
%figure;
%cimage5(h1{1}(:,:,4));
%figure;
%cimage5(h11{2}(:,:,4));
%figure;
%cimage5(h1{2}(:,:,4));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%
% Score a MotorProgram
%   log-probability score as in Equation 1 in the main text
%
% Input
%  M: motor program
%  libclass: object of type LibraryClass
%
%  By default, it includes all components of the log-score.
%  To choose only specific components, you can use this syntax:
%    'token',true,'type',true,'image',true,'stat',true
%
%     For the token, type, and image model respectively set to true or
%       false. 'stat' is only needed for alphabet model, which is not
%       included here and should not be used.
%
%  Output
%   score : score
%   out : struct that breaks down the score by component
%
function [score,out] = scoreMP(M,libclass,varargin)

    % should we include type and/or token-level variables?
    assert(isa(libclass,'Library'));
    n = length(varargin);
    if mod(n,2)~=0
       error('the number of optional arguments must be even');
    end

    defaultinc = (n == 0 || (n==2 && strcmp(varargin{1},'strokes')));
    inc_stat = false;
    if defaultinc
        inc_image = true;
        inc_token = true;
        inc_type = true;
    else % specified in the arguments
        inc_image = nan;
        inc_token = nan;
        inc_type = nan;
    end
    list_sid = 1:M.ns;
    for i=1:2:n-1
       switch varargin{i}
           case 'strokes'
               list_sid = varargin{i+1};
               assert(isnumeric(list_sid));
           case 'token'
               inc_token = varargin{i+1};
               assert(islogical(inc_token));
           case 'type'
               inc_type = varargin{i+1};
               assert(islogical(inc_type));
           case 'image'
               inc_image = varargin{i+1};
               assert(islogical(inc_image));
           case 'stat'
               inc_stat = varargin{i+1};
               assert(islogical(inc_stat));
           otherwise
               error('invalid input parameters to scoreMP');
       end
    end

    % check that we have specified all the terms to include/exclude
    if ~defaultinc
        if any(isnan([inc_image inc_token inc_type]))
           error('score terms not fully specified');
        end
    end
    PM = M.parameters;
    out = struct;

    % check that we have relations instantiated
    if inc_token || inc_type
        assert(M.has_relations(list_sid));
    end

    % type-level score for number of strokes
    if inc_type
        out.number = CPD.score_number(libclass,M.ns);
    end

    out.S = cell(M.ns,1);
    for i=list_sid
       MS = M.S{i};
       out.S{i} = struct;

       % type-level scores
       if inc_type
           out.S{i}.sequence = CPD.score_sequence(libclass,M.ns,MS.ids);
           out.S{i}.scales_type = CPD.score_invscale_type(libclass,MS.invscales_type,MS.ids);
           out.S{i}.relation = CPD.score_relation_type(libclass,MS.R);
           if ~isempty(MS.shapes_type)
              out.S{i}.shape_type = CPD.score_shape_type(libclass,MS.shapes_type,MS.ids);
           end
       end

       % token-level scores
       if inc_token
           if isempty(MS.shapes_type)
               out.S{i}.shape_marginalize = CPD.score_shape_marginalize(libclass,MS.shapes_token,MS.ids);
           else
               out.S{i}.shape_token = CPD.score_shape_token(libclass,MS.shapes_token,MS.shapes_type);
           end
           out.S{i}.scales_token = CPD.score_invscale_token(libclass,MS.invscales_token,MS.invscales_type);
           out.S{i}.position = CPD.score_position(libclass,MS.pos_token,MS.R,M.S(1:i-1));
           if strcmp(MS.R.type,'mid')
               if isempty(MS.R.eval_spot_type)
                   out.S{i}.relation_token = CPD.score_relation_token_approx_marginalize(libclass,MS.R.eval_spot_token);
               else
                   out.S{i}.relation_token = CPD.score_relation_token(libclass,MS.R.eval_spot_token,MS.R.eval_spot_type);
               end
           end
       end

    end

    % token-level image variables
    if inc_token
        out.image_blur = CPD.score_image_blur(M.blur_sigma,PM);
        out.image_noise = CPD.score_image_noise(M.epsilon,PM);
        out.A = CPD.score_affine(libclass,M.A);
    end

    % include the stat penalty term
    if inc_stat
        error('This term should not be used in the score.');
        out.stat = CPD.score_stat(libclass,M);
    end

    % score the image
    if inc_image
        out.image = CPD.score_image(M.I,M.pimg);
    end

    score = sumfields(out);
    if isnan(score)
       assert(false);
    end
end

function total = sumfields(S)
%
% Add up all the scalar, numeric fields
% in a structure S. Recursively call
% on cell arrays of fields
%
    if isempty(S)
       total = 0;
       return
    end
    assert(isstruct(S));
    total = 0;
    fds = fieldnames(S);
    for i=1:length(fds)
       field = S.(fds{i});
       if ~isempty(field)

          if isnumeric(field) % if number
              total = total + sum(field(:));

          elseif iscell(field) % if cell, call recursively
              for j=1:length(field)
                total = total + sumfields(field{j});
              end

          else
              error('invalid field type');
          end

       end
    end
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [  net,res,opts ] = rmsprop(  net,res,opts )
% Modified RMSProp using second-order information.
%   1.Tieleman, T. and Hinton, G. Lecture 6.5 - RMSProp, COURSERA: Neural Networks for Machine Learning.
%   Technical report, 2012.
%   2.Ye, C., Yang, Y., Fermuller, C., & Aloimonos, Y. (2017).
%   On the Importance of Consistency in Training Deep Neural Networks. arXiv preprint arXiv:1708.00631.

    if ~isfield(opts.parameters,'second_order')
        opts.parameters.second_order=0;
    end
    if opts.parameters.second_order
        [  net,res,opts ] = gradient_decorrelation(  net,res,opts );
    end

    if ~isfield(opts.parameters,'weightDecay')
        opts.parameters.weightDecay=1e-4;
    end

    if ~isfield(opts.parameters,'clip')
        opts.parameters.clip=1e0;
    end

    if ~isfield(opts.parameters,'eps')
        opts.parameters.eps=1e-6;
    end

    if ~isfield(net,'iterations')||(isfield(opts,'reset_mom')&&opts.reset_mom==1)
        net.iterations=0;
    end

    net.iterations=net.iterations+1;

    mom_factor=(1-opts.parameters.mom.^net.iterations);

    for layer=1:numel(net.layers)
        if isfield(net.layers{layer},'weights')
            if ~isfield(net.layers{layer},'momentum')||(isfield(opts,'reset_mom')&&opts.reset_mom==1)
                net.layers{layer}.momentum{1}=zeros(size(net.layers{layer}.weights{1}),'like',net.layers{layer}.weights{1});
                net.layers{layer}.momentum{2}=zeros(size(net.layers{layer}.weights{2}),'like',net.layers{layer}.weights{2});

            end

            net.layers{layer}.momentum{1}=opts.parameters.mom.*net.layers{layer}.momentum{1}+(1-opts.parameters.mom).*(res(layer).dzdw.^2);
            normalized_grad=res(layer).dzdw./(net.layers{layer}.momentum{1}.^0.5+opts.parameters.eps)./mom_factor;
            if isfield(opts.parameters,'clip')&&opts.parameters.clip>0
                mask=abs(normalized_grad)>opts.parameters.clip;
                normalized_grad(mask)=sign(normalized_grad(mask)).*opts.parameters.clip;
            end
            net.layers{layer}.weights{1}=net.layers{layer}.weights{1}-opts.parameters.lr*normalized_grad- opts.parameters.weightDecay * net.layers{layer}.weights{1};

            net.layers{layer}.momentum{2}=opts.parameters.mom.*net.layers{layer}.momentum{2}+(1-opts.parameters.mom).*(res(layer).dzdb.^2);
            normalized_grad=res(layer).dzdb./(net.layers{layer}.momentum{2}.^0.5+opts.parameters.eps)./mom_factor;
            if isfield(opts.parameters,'clip')&&opts.parameters.clip>0
                mask=abs(normalized_grad)>opts.parameters.clip;
                normalized_grad(mask)=sign(normalized_grad(mask)).*opts.parameters.clip;
            end
            net.layers{layer}.weights{2}=net.layers{layer}.weights{2}-opts.parameters.lr*normalized_grad;
        end
    end

   if ~isfield(opts,'reset_mom')||opts.reset_mom==1
        opts.reset_mom=0;
    end
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function varargout = grSimplifyBranches(nodes, edges)
%GRSIMPLIFYBRANCHES Replace branches of a graph by single edges
%
%   [NODES2 EDGES2] = grSimplifyBranches(NODES, EDGES)
%   renvoie une version simplifiee d'un graphe, en ne gardant que les
%   points multiples et les aretes reliant les points multiples.
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY :
%   10/02/2004 : doc
%   17/01/2006 : uses faster method to find neighbour edges
%   18/01/2006 : replace call to subfunctions by inlining -> faster

Mnodes = [];    % size Nn*2 -> nodes coordinates
Sedges = [];    % size Ne*2 -> indices of nodes
Mpoints = [];   % size Nn*1 -> indices of Multiple points
                %     in nodes input array

branch = [];    % size Nb*2 (variable)

Nn = 0;
Ne = 0;
Nb = 0;

% look for the first multiple point
p = 1;
while length(find(edges(:,1) == p | edges(:,2) == p)) < 3
    p = p + 1;
end

Mpoints(1) = p;
Mnodes(1, 1:2) = nodes(p, 1:2);
Nn = Nn + 1;

% add the branches of the first multiple point
neighbours = find(edges(:,1)==p | edges(:,2)==p);
for b = 1:length(neighbours)
    Nb = Nb+1;
    edge = edges(neighbours(b),:);
    if edge(1) == p
        branch(Nb, 1:2) = [p edge(2)]; %#ok<AGROW>
    else
        branch(Nb, 1:2) = [p edge(1)]; %#ok<AGROW>
    end
end

% process each branch, until there is no more branch to process.
% b is index of current branch
b = 0;
while b < length(branch)
    b = b+1;

    % check if the branch is valid
    if branch(b, 1) == 0
        continue;
    end

    % initialize iteration
    pNode = branch(b, 1);
    node = branch(b,2);
    neighbours = find(edges(:,1) == node | edges(:,2) == node);

%     disp(sprintf('node %3d (%03d ; %03d) -> %3d (%03d ; %03d)', ...
%         Mnode, nodes(Mnode, 1), nodes(Mnode, 2), ...
%         node, nodes(node, 1), nodes(node, 2)));

    while length(neighbours) < 3
        % look for the next point on the current branch
        next = 0;
        for n = 1:length(neighbours)
            edge = edges(neighbours(n), :);
            if edge(1)~= node && edge(1)~= pNode
                next = edge(1);
                break;
            end
            if edge(2)~= node && edge(2)~= pNode
                next = edge(2);
                break;
            end
        end

        pNode = node;
        node = next;
        neighbours = find(edges(:,1) == node | edges(:,2) == node);
    end

    % node is now  the next multiple point, and pNode contains the last
    % point of the branch.

    % check if the multiple point has already been processed
    index = find(Mpoints==node);
    if ~isempty(index)
        % find the branch starting with node, and with pNode has
        % second point, and set it to [0 0] to avoid it to be
        % processed again
        %disp('remove branch');
        for b2 = 1:Nb
            if branch(b2, 1) == node && branch(b2, 2) == pNode
                %disp('find branch');
                branch(b2, 1:2) = 0; %#ok<AGROW>
                break;
            end
        end

    else
        % add the multiple point to the list of points
        %disp('add point');
        Nn = Nn+1;
        Mnodes(Nn, 1:2) = nodes(node, 1:2); %#ok<AGROW>
        index = Nn;
        Mpoints(Nn) = node; %#ok<AGROW>

        % add each neighbour of the new multiple point (but not
        % the neighbour containing pNode) to the list of branches
        for n = 1:length(neighbours)
            edge = edges(neighbours(n), :);
            if edge(1) ~= pNode && edge(2) ~= pNode
                %disp('add a branch');
                Nb = Nb + 1;
                if edge(1) == node
                    branch(Nb, 1:2) = [node edge(2)]; %#ok<AGROW>
                else
                    branch(Nb, 1:2) = [node edge(1)]; %#ok<AGROW>
                end
            end
        end
    end

    %disp('add new Edge');
    Ne = Ne + 1;
    Sedges(Ne, 1:2) = [find(Mpoints == branch(b,1)), index]; %#ok<AGROW>
end

% process output depending on how many arguments are needed
if nargout == 1
    out{1} = Mnodes;
    out{2} = Sedges;
    varargout{1} = out;
end

if nargout == 2
    varargout{1} = Mnodes;
    varargout{2} = Sedges;
end

return;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function r = resolution(a)

r = a.resolution;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [p2, R2] = SE3exp(j, dt)
global uLINK
% see Murray, Li, Sastry p.42

norm_w = norm(uLINK(j).w);
if norm_w < eps
    p2 = uLINK(j).p + dt * uLINK(j).vo;
    R2 = uLINK(j).R;
else
    th = norm_w*dt;
    wn = uLINK(j).w/norm_w;		% normarized vector
    vo = uLINK(j).vo/norm_w;

    w_wedge = [0 -wn(3) wn(2);wn(3) 0 -wn(1);-wn(2) wn(1) 0];
    drot = w_wedge * sin(th) + w_wedge^2 * (1-cos(th));
    rot  = eye(3) + drot;

    p2 = rot * uLINK(j).p -drot*cross(wn, vo) + wn * wn' * vo * th;
    R2 = rot * uLINK(j).R;
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function rss=get_rss_from_file(url)
csi_trace = read_bf_file(url);
csi_trace = get_csi_trace_sqeezed(csi_trace);
rss = get_rss_from_csitrace(csi_trace);
end

author: linteresa

repo: WiAR

save_path: github-repos/MATLAB/linteresa-WiAR


text:
%% ASCII Art
%  The ASCII art is an interesting art of representing images using
% characters. In the past a great number of powerful libs and tools are
% develped, the most famous is aalib that allow to paint in real time a
% stream of images (a video) containing only ASCII characters (on a TTY).
% So for example the using the mplayer under Linux it is possible to see a
% DVD film completely in ASCII art, or to play Quake with it's aalib version
% (tty-quake). Here a simple ASCII art converter function is given using
% multiple palettes depending on the gradient of the image, something
% different with respecto to the other AA tools.

%% Reading an image
%  Reading and preparing an image for the conversion: an image must be
% gayscale to be used with gabAscii.

% Getting a standard image:
img = imread('peppers.png');

% Converting to grayscale:
img = rgb2gray(img);

%% Generating the ASCII-Art version of the image
%  Here the ASCII-Art version of the image is generated using the default
% parameters, in this tool the size of the axels and the palettes can be
% defined as parameters.

% Generation of the output:
asciiImg = gabAscii(img),

author: Sable

repo: mcbench-benchmarks


text:
function [varargout]=wfdbloadlib(varargin)
%
% [isloaded,config]=wfdbloadlib(debugLevel,networkWaitTime)
%
% Loads the WDFDB libarary if it has not been loaded already into the
% MATLAB classpath. And optionally prints configuration environment and debug information
% regarding the settings used by the classes in the JAR file.
%
% Inputs:
%
% debugLevel
%       (Optional) 1x1 integer between 0 and 5 represeting the level of debug information to output from
%       Java class when output configuration information. Level 0 (no debug information),
%       level =5 is maximum level of information output by the class (logger set to finest). Default is level 0.
%
% networkWaitTime
%       (Optional) 1x1 integer representing the longest time in
%       milliseconds  for which the JVM should wait for a data stream from
%       PhysioNet (default is =1000  , ie one second). If you need to change this time to a
%       longer value across the entire toolbox, it is better modify to default value in the source
%       code below and restart MATLAB.
%
%
% Written by Ikaro Silva, 2013
%         Last Modified: April 7, 2015
% Since 0.0.1
%
%

%endOfHelp
mlock;
persistent isloaded wfdb_path wfdb_native_path config

%%%%% SYSTEM WIDE CONFIGURATION PARAMETERS %%%%%%%
%%% Change these values for system wide configuration of the WFDB binaries

%If you are using your own custom version of the WFDB binaries, set this to true
%NOTE: this parameter is completely ignored if the 'WFDB_COMMAND_PATH' parameter
%described above is set (i.e.: the library will used the WFDB commands located
% according to the path in 'WFDB_COMMAND_PATH').
%You will need to restart MATLAB/Octave if to sync the changes.
%The default is to used commands shipped with the toolbox, this location can be obtained by running the command:
%[~,config]=wfdbloadlib; config.WFDB_NATIVE_BIN
WFDB_CUSTOMLIB=0;

%WFDB_PATH: If empty, will use the default given config.WFDB_PATH
%this is where the toolbox searches  for data files (*.dat, *.hea etc).
%When unistalling the toolbox, you may wish to clear this directory to save space.
%See http://www.physionet.org/physiotools/wag/setwfd-1.htm for more details.
WFDB_PATH=[];

%WFDBCAL: If empty, will use the default giveng confing.WFDBCAL
%The WFDB library require calibration data in order to convert between sample values
%(expressed in analog-to-digital converter units, or adus) and physical units.
%See http://www.physionet.org/physiotools/wag/wfdbca-5.htm for more details.
WFDBCAL=[];

%CACHE: If CACHE==1, the toolbox will attemp to download data from
%CACHE_SOURCE to CACHE_DEST if the record is not found no the standard
%WFDB PATH. Change CACHE_DEST path to a PhysioNet mirror, if you wish to
%use a server closer to your geographical location. It is safe to delete files on
%CACHE_SOURCE, as they can be re-downloaded if need me.

CACHE=1; %Default is to use the cache system
CACHE_SOURCE=[]; %If empty, defaults to last element of WFDB_PATH
CACHE_DEST=[]; %If empty, defaults to WFDB_JAVA_HOME/../database

%debugLevel: Ouput JVM information while running commands
debugLevel=0;

%networkWaitTime: Setting maximum waiting period for fetching data from
%PhysioNet servers (default location: http://physionet.org/physiobank).
networkWaitTime=1000;

%%%% END OF SYSTEM WIDE CONFIGURATION PARAMETERS

%If called without arguments, keep the current configuration.
%Explicit arguments override the default settings above.
if(~isempty(config))
    debugLevel=config.DEBUG_LEVEL;
    networkWaitTime=config.NETWORK_WAIT_TIME;
end

inputs={'debugLevel','networkWaitTime'};
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

inOctave=is_octave;
fsep=filesep;
if(ispc && inOctave)
	fsep=['\\']; %Need to escape '\' for regexp in Octave and Windows
end
if(isempty(isloaded))
    jar_path=which('wfdbloadlib');
    cut=strfind(jar_path,'wfdbloadlib.m');
    wfdb_path=jar_path(1:cut-1);

    if(~inOctave)
        ml_jar_version=version('-java');
    else
        %In Octave
        ml_jar_version=javaMethod('getProperty','java.lang.System','java.version');
        ml_jar_version=['Java ' ml_jar_version];
    end
    %Check if path has not been added yet
    wfdb_path=[wfdb_path 'wfdb-app-JVM7-0-10-0.jar'];
    javaaddpath(wfdb_path);
    isloaded=1;

    %Check if there are any empty space on the path directory, and
    %issue a warning if there is
    warnMe=strfind(wfdb_path,' ');
    if(~isempty(warnMe))
       warning('Your WFDB Toolbox installation  path contain white spaces!! This may cause issues with the WFDB Toolbox!');
       warning(['The installation path is set to: ' wfdb_path]);
    end

    if(isunix)
        % work around the fact that glibc requires thread stack size
        % to include space for thread-local variables
        % (https://sourceware.org/bugzilla/show_bug.cgi?id=11787)
        % and JVM "process reaper" threads use a tiny stack by default
        javaMethod('setProperty','java.lang.System','jdk.lang.processReaperUseDefaultStackSize','true');
    end
end

%set configuration
if(isempty(config))
        config.MATLAB_VERSION=version;
        config.inOctave=inOctave;
        if(inOctave)
            javaWfdbExec=javaObject('org.physionet.wfdb.Wfdbexec','wfdb-config',WFDB_CUSTOMLIB);
            javaWfdbExec.setLogLevel(debugLevel);
            config.WFDB_VERSION=char(javaMethod('execToStringList',javaWfdbExec,{'--version'}));
        else
            javaWfdbExec=org.physionet.wfdb.Wfdbexec('wfdb-config',WFDB_CUSTOMLIB);
            javaWfdbExec.setLogLevel(debugLevel);
            config.WFDB_VERSION=char(javaWfdbExec.execToStringList('--version'));
        end
        env=regexp(char(javaWfdbExec.getEnvironment),',','split');
        for e=1:length(env)
            tmpstr=regexp(env{e},'=','split');
            varname=strrep(tmpstr{1},'[','');
            varname=strrep(varname,' ','');
            varname=strrep(varname,']','');
            eval(['config.' varname '=''' tmpstr{2} ''';']);
        end
        config.MATLAB_PATH=strrep(which('wfdbloadlib'),'wfdbloadlib.m','');
        wver=regexp(wfdb_path,fsep,'split');
        config.WFDB_JAVA_VERSION=wver{end};
        config.MATLAB_ARCH=computer('arch');
        %Remove empty spaces from arch name
        del=strfind(config.osName,' ');
        config.osName(del)=[];

        %Define WFDB Environment variables
        if(isempty(WFDB_PATH))
            tmpCache=[config.MATLAB_PATH '..' filesep 'database' filesep];
            WFDB_PATH=['. ' tmpCache ' http://physionet.org/physiobank/database/'];
        end
        if(isempty(WFDBCAL))
            WFDBCAL=[config.WFDB_JAVA_HOME fsep 'database' fsep 'wfdbcal'];
        end
        config.WFDB_PATH=WFDB_PATH;
        config.WFDBCAL=WFDBCAL;
        config.WFDB_CUSTOMLIB=WFDB_CUSTOMLIB;
            warnMe=strfind(wfdb_path,' ');
    if(~isempty(warnMe))
       warning('Your WFDB Toolbox installation  path contain white spaces!! This may cause issues with the WFDB Toolbox!');
    end

    %Set CACHE configurations
    if(isempty(CACHE_SOURCE) && CACHE)
        ind=strfind(config.WFDB_PATH,'http');
        if(~isempty(ind))
            CACHE_SOURCE=config.WFDB_PATH(ind:end);
        else
            warning(['Could not set CACHE, CACHE_SOURCE invalid']);
            CACHE=0;
        end
    end
    config.CACHE_SOURCE=CACHE_SOURCE;

    if(isempty(CACHE_DEST) && CACHE)
        CACHE_DEST=[config.MATLAB_PATH '..' filesep 'database' filesep];
        if(~isdir(CACHE_DEST))
            mkdir(CACHE_DEST);
        end
        if(~isdir(CACHE_DEST))
            warning(['Could not set CACHE, CACHE_DEST directory does not exist: ' CACHE_DEST]);
            CACHE=0;
        end
    end
    config.CACHE_DEST=CACHE_DEST;
    config.CACHE=CACHE;

    %Set enviroment variables used by WFBD
    setenv('WFDB',config.WFDB_PATH);
    setenv('WFDBCAL',config.WFDBCAL);
end

%Set dynamic configuration parameters
config.DEBUG_LEVEL=debugLevel;
config.NETWORK_WAIT_TIME=networkWaitTime;

outputs={'isloaded','config'};
for n=1:nargout
    eval(['varargout{n}=' outputs{n} ';']);
end

%% subfunction that checks if we are in octave
function r = is_octave ()
    r = exist ('OCTAVE_VERSION', 'builtin')>0;

author: ikarosilva

repo: wfdb-app-toolbox


text:
% Demonstration of generative model functions.
%
% See GENERATIVE_MODEL and EVALUATE_GENERATIVE_MODEL for further details
% and interpretation.

clear
close all
clc

data = load('demo_generative_models_data');
A     = data.A;
Aseed = data.Aseed;
D     = data.D;

% get cardinality of network
n = length(A);

% set model type
modeltype = 'matching';

% set whether the model is based on powerlaw or exponentials
modelvar = [{'powerlaw'},{'powerlaw'}];

% choose some model parameters
nparams = 100;
params = [unifrnd(-10,0,nparams,1), unifrnd(-1,1,nparams,1)];

% generate synthetic networks and energy for the neighbors model;
[B,E,K] = evaluate_generative_model(Aseed,A,D,modeltype,modelvar,params);
X = [E,K];

% show scatterplot of parameter values versus energy and KS statistics
names = [...
    {'energy'},...
    {'degree'},...
    {'clustering'},...
    {'betweenness'},...
    {'edge length'}];

f = figure(...
    'units','inches',...
    'position',[2,2,4,4]);
for i = 1:size(X,2)
    subplot(3,2,i);
    scatter(params(:,1),params(:,2),100,X(:,i),'filled');
    set(gca,...
        'clim',[0,1]);
    colormap(jet);
    xlabel('geometric parameter, \eta');
    ylabel('topological parameter, \gamma');
    title(names{i});
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Results, pz_d] = CD_PLSA(TrainData, TestData, TrainLabel, TestLabel, numCluster, numIter, numSource, numTrain, numTarget, numTest)

% function [Results, pz_d] = CD_PLSA(Train_Data,Test_Data,Parameter_Setting)

% The common program for CD_PLSA, which can deal with multiple classes,
% multiple source domains and multiple target domains

%%%% Input:
% The parameter Train_data stores the file pathes of training data and the
% corresponding labels
% The parameter Test_data stores the file pathes of test data and the
% corresponding labels
% The parameter Parameterfile stores the parameter setting information

%%%% Output
% The variable Results is a matrix with size numIteration x numTarget, where
% numIteration is the number of iterations, numTarget is the number of
% target domains. Results record the detailed results of each iteration.

% The variable pz_d is a matrix with size n x c, where n is the number of
% instances in all target domains, specifically, n = n_1 + ... + nt (n_t is
% the number of instances in t-th target domain), c is the number of
% classes
%
% Note that if you want to deal with large data set, you should set larget
% memory for Matlab. You can set it in the file C:\boot.ini (This may not
% be true in your system), change '/fastdetect' to '/fastdetect /3GB'.
%
% Be good luck for your research, if you have any questions, you can
% contact the email: zhuangfz@ics.ict.ac.cn

TrainX = TrainData;
TrainY = TrainLabel;
TestX = TestData;
TestY = TestLabel;
numK =  numCluster;
numIteration = numIter;
labelset = union(TrainLabel,[]);
numC = length(labelset);
pyz = rand(numK,numC);
%pyz = ones(numK,numC);
pyz = pyz/sum(sum(pyz));

start = 1;
if start == 1
    DataSetX = [TrainX TestX];
    Learn.Verbosity = 1;
    Learn.Max_Iterations = 200;
    Learn.heldout = .1; % for tempered EM only, percentage of held out data
    Learn.Min_Likelihood_Change = 1;
    Learn.Folding_Iterations = 20; % for TEM only: number of fiolding in iterations
    Learn.TEM = 0; %tempered or not tempered
    [Pw_z,Pz_d,Pd,Li,perp,eta] = pLSA(DataSetX,[],numK,Learn); %start PLSA
    %xlswrite(strcat('pwz_','common_selected','.xls'),Pw_z);
end

%pwy = xlsread(strcat('pwz_','common_selected','.xls'));

pwy = Pw_z;

pw_yc_s = [];
for i = 1:numSource
    pw_yc_s = [pw_yc_s, pwy];
end
pw_yc_t = [];
for i = 1:numTarget
    pw_yc_t = [pw_yc_t, pwy];
end
clear pwy;

pd_zc_s = [];
for i = 1:numSource
    A = zeros(numTrain(i),numC);
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTrain(t);
        end
    end
    if i == 1
        pos = 0;
    end
    for j = 1:numTrain(i)
        for k = 1:numC
            if TrainY(pos+j) == labelset(k)
                A(j,k) = 1;
            end
        end
    end
    for k = 1:numC
        A(:,k) = A(:,k)/sum(A(:,k));
    end
    pd_zc_s = [pd_zc_s;A];
end

% random initialization for pd_z_t
% In our paper, pd_z_t is assigned as the predicted results by supervised classifiers
pd_zc_t = [];
for i = 1:numTarget
    A = ones(numTest(i),numC);
    for k = 1:numC
        A(:,k) = A(:,k)/sum(A(:,k));
    end
    pd_zc_t = [pd_zc_t;A];
end

pc = zeros(1,numSource+numTarget);
numAll = sum(numTrain) + sum(numTest);
for i = 1:numSource
    pc(i) = numTrain(i)/numAll;
end
for i = 1:numTarget
    pc(i+numSource) = numTest(i)/numAll;
end

iter_results = [];

O_s = TrainX;
O_t = TestX;

for iterID = 1:numIteration

    temp_pw_yc_s = [];
    temp_pw_yc_t = [];
    temp_pd_zc_s = [];
    temp_pd_zc_t = [];
    temp_pyz = zeros(size(pyz));
    temp_pc = zeros(size(pc));

    stepLen = 3000;
    numStep = fix(size(TrainX,1)/stepLen);
    step = fix(size(TrainX,1)/numStep);

    for i = 1:numSource
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTrain(t);
            end
        end
        if i == 1
            pos = 0;
        end

        A = [];
        D = zeros(numK,numTrain(i));
        for stepID = 1:numStep
            if stepID < numStep
                tempsum2_s = pw_yc_s((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK)*pyz*pd_zc_s(pos+1:pos+numTrain(i),:)';
                tempsum2_s = tempsum2_s*pc(i);
                [xs ys] = find(tempsum2_s == 0);
                for q = 1:size(xs,1)
                    tempsum2_s(xs(q,1),ys(q,1)) = 1;
                end
                B = O_s((stepID-1)*step+1:stepID*step,pos+1:pos+numTrain(i))./tempsum2_s;
                A = [A; pw_yc_s((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK).*(B*pd_zc_s(pos+1:pos+numTrain(i),:)*pyz')];
                D = D + pw_yc_s((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK)'*B;
                clear B;
            end
            if stepID == numStep
                tempsum2_s = pw_yc_s((stepID-1)*step+1:size(O_s,1),(i-1)*numK+1:i*numK)*pyz*pd_zc_s(pos+1:pos+numTrain(i),:)';
                tempsum2_s = tempsum2_s*pc(i);
                [xs ys] = find(tempsum2_s == 0);
                for q = 1:size(xs,1)
                    tempsum2_s(xs(q,1),ys(q,1)) = 1;
                end
                B = O_s((stepID-1)*step+1:size(O_s,1),pos+1:pos+numTrain(i))./tempsum2_s;
                A = [A; pw_yc_s((stepID-1)*step+1:size(O_s,1),(i-1)*numK+1:i*numK).*(B*pd_zc_s(pos+1:pos+numTrain(i),:)*pyz')];
                D = D + pw_yc_s((stepID-1)*step+1:size(O_s,1),(i-1)*numK+1:i*numK)'*B;
                clear B;
            end
        end

        A = A*pc(i);
        for l = 1:numK
            A(:,l) = A(:,l)/sum(A(:,l));
        end
        temp_pw_yc_s = [temp_pw_yc_s A];
        clear A;

        temp_pyz = temp_pyz + (pyz.*(D*pd_zc_s(pos+1:pos+numTrain(i),:)))*pc(i);
        clear D;
        temp_pc(i) = sum(sum(TrainX(:,pos+1:pos+numTrain(i))));
    end

    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end

        A = [];
        C = zeros(numTest(i),numK);
        D = zeros(numK,numTest(i));
        for stepID = 1:numStep
            if stepID < numStep
                tempsum2_t = pw_yc_t((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK)*pyz*pd_zc_t(pos+1:pos+numTest(i),:)';
                tempsum2_t = tempsum2_t*pc(numSource+i);
                [xs ys] = find(tempsum2_t == 0);
                for q = 1:size(xs,1)
                    tempsum2_t(xs(q,1),ys(q,1)) = 1;
                end
                B = O_t((stepID-1)*step+1:stepID*step,pos+1:pos+numTest(i))./tempsum2_t;
                A = [A; pw_yc_t((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK).*(B*pd_zc_t(pos+1:pos+numTest(i),:)*pyz')];
                C = C + B'*pw_yc_t((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK);
                D = D + pw_yc_t((stepID-1)*step+1:stepID*step,(i-1)*numK+1:i*numK)'*B;
                clear B;
            end
            if stepID == numStep
                tempsum2_t = pw_yc_t((stepID-1)*step+1:size(O_t,1),(i-1)*numK+1:i*numK)*pyz*pd_zc_t(pos+1:pos+numTest(i),:)';
                tempsum2_t = tempsum2_t*pc(numSource+i);
                [xs ys] = find(tempsum2_t == 0);
                for q = 1:size(xs,1)
                    tempsum2_t(xs(q,1),ys(q,1)) = 1;
                end
                B = O_t((stepID-1)*step+1:size(O_t,1),pos+1:pos+numTest(i))./tempsum2_t;
                A = [A; pw_yc_t((stepID-1)*step+1:size(O_t,1),(i-1)*numK+1:i*numK).*(B*pd_zc_t(pos+1:pos+numTest(i),:)*pyz')];
                C = C + B'*pw_yc_t((stepID-1)*step+1:size(O_t,1),(i-1)*numK+1:i*numK);
                D = D + pw_yc_t((stepID-1)*step+1:size(O_t,1),(i-1)*numK+1:i*numK)'*B;
                clear B;
            end
        end

        A = A*pc(numSource+i);
        for l = 1:numK
            A(:,l) = A(:,l)/sum(A(:,l));
        end
        temp_pw_yc_t = [temp_pw_yc_t A];
        clear A;

        A = pd_zc_t(pos+1:pos+numTest(i),:).*(C*pyz);
        clear C;
        A = A*pc(numSource+i);
        for k = 1:numC
            A(:,k) = A(:,k)/sum(A(:,k));
        end
        temp_pd_zc_t = [temp_pd_zc_t;A];

        temp_pyz = temp_pyz + (pyz.*(D*pd_zc_t(pos+1:pos+numTest(i),:)))*pc(numSource+i);
        clear D;
        temp_pc(numSource+i) = sum(sum(TestX(:,pos+1:pos+numTest(i))));
    end

    temp_pyz = temp_pyz/sum(sum(temp_pyz));
    temp_pc = temp_pc/sum(temp_pc);

    pw_yc_s = temp_pw_yc_s;
    pw_yc_t = temp_pw_yc_t;
    pd_zc_t = temp_pd_zc_t;
    pyz = temp_pyz;
    pc = temp_pc;

    pz_d = [];
    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end

        pzd = zeros(numTest(i),numC);
        for j = 1:size(pzd,1)
            for k = 1:size(pzd,2)
                pzd(j,k) = pd_zc_t(pos+j,k)*sum(pyz(:,k));
            end
        end
        pz_d = [pz_d; pzd];
        nCorrect = 0;
        for j = 1:size(pzd,1)
            [va vi] = max(pzd(j,:));
            if labelset(vi) == TestY(pos+j)
                nCorrect = nCorrect + 1;
            end
        end
        acc = [double(iterID) nCorrect/(numTest(i))]
        iter_results(iterID,i) = nCorrect/(numTest(i));
    end
    for i = 1:size(pz_d,1)
         pz_d(i,:) = pz_d(i,:)/sum( pz_d(i,:));
    end

    O_s = TrainX;
    O_t = TestX;
end

Results = iter_results;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function make_frames ( frame_total )

%*****************************************************************************80
%
%% MAKE_FRAMES makes the frames of an animation.
%
%  Discussion:
%
%    This program makes FRAME_TOTAL + 1 JPEG files, which constitute
%    individual frames of an animation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FRAME_TOTAL, the total number of frames.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_FRAMES:\n' );
  fprintf ( 1, '  Make a (potential) animation by creating a sequence\n' );
  fprintf ( 1, '  of JPG files, each containing one frame.\n' );

  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    frame_total = input ( 'Enter the number of frames to create:  ' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program has been asked to create\n' );
  fprintf ( 1, '  frames 0 through %d.\n', frame_total );
  fprintf ( 1, '  The first frame will be stored in "%s",\n', ...
    sprintf ( 'frame%04d.jpg', 0 ) );
  fprintf ( 1, '  The last frame will be stored in "%s".\n', ...
    sprintf ( 'frame%04d.jpg', frame_total ) );

  for frame_id = 0 : frame_total
    make_frame ( frame_id, frame_total );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_FRAMES:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function [db,f]=v_lpccc2db(cc,np,nc,c0)
%V_LPCCC2DB Convert complex cepstrum to dB power spectrum DB=(CC,NP,NC)
%
%  Inputs: cc(nf,n)     Complex ceptral coefficients excluding c(0), one frame per row
%          np           Size of output spectrum is np+1 [n]
%                       Alternatively, np can be a vector of output frequencies in the range 0 to 0.5
%          nc           Highest cepstral coefficient to use [np or, if np is a vector, n]
%                       Set nc=-1 to use n coefficients
%          c0(nf,1)     Cepstral coefficient cc(0) [0]
%
% Outputs: db(nf,np+2)  Power spectrum from DC to Nyquist in dB
%          f(1,np+2)    Normalized frequencies (0 to 0.5)
%
% The "complex cepstral coefficients", cc(n), are the inverse discrete-time Fourier transform
% of the log of the complex-valued spectrum. The cc(n) are real-valued and, for n<0, cc(n)=0.
% The "real cepstral coeffcients", rc(n), are the inverse discrete-time Fourier transform
% of the log of the magnitude spectrum; rc(0)=cc(0) and rc(n)=0.5*cc(n) for n~=0.
% For highest speed, choose np to be a power of 2.

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: v_lpccc2db.m 10865 2018-09-21 17:22:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,mc]=size(cc);
if nargin<2 || ~numel(np)
    if nargout
        np=mc;
    else
        np=128;
    end
end
k=10/log(10);
if nargin>=3 && numel(nc)==1 && nc==-1 nc=mc; end
if nargin<4 || ~numel(c0) c0=zeros(nf,1); end
if numel(np)>1 || np(1)<1
    if nargin<3 || ~numel(nc) nc=mc; end
    f=np(:)';
    if nc==mc
        db=k*(2*[c0 cc]*cos(2*pi*(0:mc)'*f));
    else
        db=k*(2*[c0 lpccc2cc(cc,nc)]*cos(2*pi*(0:nc)'*f));
    end
else
    if nargin<3 || ~numel(nc) nc=np; end
    if nc==mc
        db=k*(2*real(v_rfft([c0 cc].',2*np).'));
    else
        db=k*(2*real(v_rfft([c0 v_lpccc2cc(cc,nc)].',2*np).'));
    end
    f=linspace(0,0.5,np+1);
end
if ~nargout
    plot(f,db.');
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function varargout = svd( A, b )
%
%       [U, S, V] = svd(A);
%       [U, S, V] = svd(A, b);
%
%  Overloaded SVD method for psfMatrix objects.
%
%  Input:
%    A is a psfMatrix
%
%  Optional Input:
%    b is a blurred image.  Since A often uses a compact storage
%      scheme, this is sometimes needed to determine "real size"
%      of the matrix.
%
%  Output:
%    This depends on the type of blur, and boundary conditions.
%    * If A.boundary = 'periodic', then U and V are
%      transformMatrix objects, with U.transform V.transform = 'fft'
%      S is a column vector containing the eigenvalues of A.
%    * If A.boundary = 'reflexive', and A is symmetric, then
%      U and V are transformMatrix objects, with
%      U.transform = V.transform = 'dct'
%      S is a column vector containing the eigenvalues of A.
%    * In all other cases, a Kronecker product approximation of
%      A is first computed, and U and V are then kronMatrix objects.
%      S is a column vector containing the singular values of A
%      (they are not sorted, though).
%

%  J. Nagy  6/2/02
%  Modifications:
%    7/7/03 - J. Nagy, now allows to use a preliminary version
%                      of space variant SVD approximations.
%  22/03/07 - J. Nagy, this tries FFT, DCT and Kronecker product
%                      SVD bases to find which approximation
%                      is best.

switch A.type
  case 'invariant'
    P = A.psf;
    P1 = P.image;, c = P.center;
    PSF = P1{1};,  center = c{1};

    if (nargin == 2)
      PSF = padarray(PSF, size(b) - size(PSF), 'post');
    else
      b = PSF;  % This is used to define correct dimensions only.
    end

    switch A.boundary
    case 'periodic'
      [U, S, V] = fft_svd(PSF, center);
    case 'reflexive'
      if issymmetric(A)
         [U, S, V] = dct_svd(PSF, center);
      else
         [U, S, V] = svd_approx( kronApprox(A, b) );
      end
    case 'zero'
      [U, S, V] = svd_approx( kronApprox(A, b) );
    otherwise
      error('Invalid boundary condition')
    end

  case 'variant'
    [U, S, V] = svd_approx( kronApprox(A, b) );

  otherwise
    error('invalid matrix type')
end

if nargout == 1
  varargout{1} = S;
elseif nargout == 3
  varargout{1} = U;
  varargout{2} = S;
  varargout{3} = V;
else
  error('In correct number of output variables')
end

repo: IRtools


text:
function [R,s]=rotmat2vec(u,v)
%
% [R,s]=rotmat2vec(u,v)
%
% the rotation matrix from vector u to v, satisfying R*u*s=v
%
% author: Bruno Luong
% URL:http://www.mathworks.com/matlabcentral/newsreader/view_original/827969
%
% input:
%   u: a 3D vector in the source coordinate system;
%   v: a 3D vector in the target coordinate system;
%
% output:
%   R: a rotation matrix to transform normalized u to normalized v
%   s: a scaling factor, so that R*u(:)*s=v(:)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

s=norm(v(:))/norm(u(:));
u1=u(:)/norm(u(:));
v1=v(:)/norm(v(:));

k = cross(u1,v1);
if(~any(k)) % u and v are parallel
    R=eye(3);
    return;
end
% Rodrigues's formula:
costheta = dot(u1,v1);
R =[ 0 -k(3) k(2);
     k(3) 0 -k(1);
    -k(2) k(1) 0];
R = costheta*eye(3) + R + k*k'*(1-costheta)/sum(k.^2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Problem = UFget (matrix, UF_Index)
%UFGET loads a matrix from the UF Sparse Matrix Collection.
%
%   Problem = UFget(matrix) loads a matrix from the UF sparse matrix collection,
%   specified as either a number (1 to the # of matrices in the collection) or
%   as a string (the name of the matrix).  With no input parameters, index=UFget
%   returns an index of matrices in the collection.  A local copy of the matrix
%   is saved (be aware that as of March 2008 the entire collection is 9.5GB
%   in size).  If no input or output arguments are provided, the index is
%   printed.  With a 2nd parameter (Problem = UFget (matrix, index)), the index
%   file is not loaded.  This is faster if you are loading lots of matrices.
%
%   Examples:
%       index = UFget ;                     % loads index
%       index = UFget ('refresh') ;         % forces download of new index
%
%       Problem = UFget (6)                 % 4 ways of loading the same Problem
%       Problem = UFget ('HB/arc130')
%       Problem = UFget (6, index)
%       Problem = UFget ('HB/arc130', index)
%
%   See also UFgrep, UFweb, UFget_example, UFget_defaults, urlwrite.

%   Copyright 2008, Tim Davis, University of Florida.

%-------------------------------------------------------------------------------
% get the parameter settings
%-------------------------------------------------------------------------------

params = UFget_defaults ;
indexfile = sprintf ('%sUF_Index.mat', params.dir) ;
indexurl = sprintf ('%s/UF_Index.mat', params.url) ;

%-------------------------------------------------------------------------------
% get the index file (download a new one if necessary)
%-------------------------------------------------------------------------------

refresh = 0 ;
if nargin == 0
    % if the user passed in a zero or no argument at all, return the index file
    matrix = 0 ;
else
    % UFget ('refresh') downloads the latest index file from the web
    if (ischar (matrix))
        if (strcmp (matrix, 'refresh'))
            matrix = 0 ;
            refresh = 1 ;
        end
    end
end

if (~refresh)
    try
        % load the existing index file
        if (nargin < 2)
            load (indexfile) ;
        end
        % see if the index file is old; if so, download a fresh copy
        refresh = (UF_Index.DownloadTimeStamp + params.refresh < now) ;
    catch
        % oops, no index file, or a refresh is due.  download it.
        refresh = 1 ;
    end
end

err = '' ;      % to catch a download error, if any

if (refresh)
    % a new UF_Index.mat file to get access to new matrices (if any)
    try
        tmp = tempname ;                        % download to a temp file first
        old = sprintf ('%sUF_Index_old.mat', params.dir) ;
        urlwrite (indexurl, tmp) ;              % download the latest index file
        try
            movefile (indexfile, old, 'f') ;    % keep a backup of the old index
        catch
            % backup failed, continue anyway
        end
        movefile (tmp, indexfile, 'f') ;        % move the new index into place
    catch
        err = lasterr ;
    end
    load (indexfile) ;
    UF_Index.DownloadTimeStamp = now ;
    save (indexfile, 'UF_Index') ;
end

%-------------------------------------------------------------------------------
% return the index file if requested
%-------------------------------------------------------------------------------

if (matrix == 0)
    if (nargout == 0)
        % no output arguments have been passed, so print the index file
        fprintf ('\nUF sparse matrix collection index:  %s\n', ...
            UF_Index.LastRevisionDate) ;
        fprintf ('\nLegend:\n') ;
        fprintf ('(p,n)sym:  symmetry of the pattern and values\n') ;
        fprintf ('           (0 = unsymmetric, 1 = symmetric, - = not computed)\n') ;
        fprintf ('type:      real\n') ;
        fprintf ('           complex\n') ;
        fprintf ('           binary:  all entries are 0 or 1\n') ;
        nmat = length (UF_Index.nrows) ;
        for j = 1:nmat
            if (mod (j, 25) == 1)
                fprintf ('\n') ;
                fprintf ('ID   Group/Name                nrows-by-  ncols  nonzeros  (p,n)sym  type\n') ;
            end
            s = sprintf ('%s/%s', UF_Index.Group {j}, UF_Index.Name {j}) ;
            fprintf ('%4d %-30s %7d-by-%7d %9d ', ...
            j, s, UF_Index.nrows (j), UF_Index.ncols (j), UF_Index.nnz (j)) ;
            psym = UF_Index.pattern_symmetry (j) ;
            nsym = UF_Index.numerical_symmetry (j) ;
            if (psym < 0)
                fprintf ('  -  ') ;
            else
                fprintf (' %4.2f', psym) ;
            end
            if (nsym < 0)
                fprintf ('  -  ') ;
            else
                fprintf (' %4.2f', nsym) ;
            end
            if (UF_Index.isBinary (j))
                fprintf (' binary\n') ;
            elseif (~UF_Index.isReal (j))
                fprintf (' complex\n') ;
            else
                fprintf (' real\n') ;
            end
        end
    else
        Problem = UF_Index ;
    end

    if (~isempty (err))
        fprintf ('\nUFget: unable to download latest index; using old one.\n') ;
        disp (err) ;
    end
    return ;
end

%-------------------------------------------------------------------------------
% determine if the matrix parameter is a matrix index or name
%-------------------------------------------------------------------------------

[group matrix id] = UFget_lookup (matrix, UF_Index) ;

if (id == 0)
    error ('invalid matrix') ;
end

%-------------------------------------------------------------------------------
% download the matrix (if needed) and load it into MATLAB

matdir = sprintf ('%s%s%s%s.mat', params.dir, group) ;
matfile = sprintf ('%s%s%s.mat', matdir, filesep, matrix) ;
maturl = sprintf ('%s/%s/%s.mat', params.url, group, matrix) ;

if (~exist (matdir, 'dir'))
    mkdir (matdir) ;                        % create the Group directory
end

if (~exist (matfile, 'file'))
    fprintf ('downloading %s\n', maturl) ;
    fprintf ('to %s\n', matfile) ;
    tmp = tempname ;                        % download to a temp file first
    urlwrite (maturl, tmp) ;
    movefile (tmp, matfile, 'f') ;          % move the new matrix into place
end

load (matfile) ;


text:
function [filename] = imgFilename(iid)
% function [filename] = imgFileName(iid)
%
% Return image filename.
%
% INPUT
%	iid		Image ID.
%
% OUTPUT
%	filename	Image filename.
%
% David Martin <dmartin@eecs.berkeley.edu>
% January 2003

filename = fullfile(bsdsRoot,'images','test',sprintf('%d.jpg',iid));
if length(dir(filename))==1, return; end

filename = fullfile(bsdsRoot,'images','train',sprintf('%d.jpg',iid));
if length(dir(filename))==1, return; end

error(sprintf('Could not find image %d in %s/images/{train,test}.', ...
              iid,bsdsRoot));

author: Cloud-CV

repo: object-proposals


text:
%% Laplace Operator
%
% In this demo, we show how to use the OpenCV function |cv.Laplacian| to
% implement a discrete analog of the _Laplacian operator_.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d5/db5/tutorial_laplace_operator.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/Laplace_Demo.cpp>
%

%% Theory
%
% In a previous demo, we learned how to use the _Sobel Operator_. It was based
% on the fact that in the edge area, the pixel intensity shows a "jump" or a
% high variation of intensity. Getting the first derivative of the intensity,
% we observed that an edge is characterized by a maximum, as it can be seen in
% the figure:
%
% <<https://docs.opencv.org/3.2.0/Laplace_Operator_Tutorial_Theory_Previous.jpg>>
%
% And what happens if we take the second derivative?
%
% <<https://docs.opencv.org/3.2.0/Laplace_Operator_Tutorial_Theory_ddIntensity.jpg>>
%
% You can observe that the second derivative is zero! So, we can also use this
% criterion to attempt to detect edges in an image. However, note that zeros
% will not only appear in edges (they can actually appear in other meaningless
% locations); this can be solved by applying filtering where needed.
%
%% Laplacian Operator
%
% From the explanation above, we deduce that the second derivative can be used
% to _detect edges_. Since images are "2D", we would need to take the
% derivative in both dimensions. Here, the Laplacian operator comes handy.
%
% The _Laplacian operator_ is defined by:
%
% $$Laplace(f) = \frac{\partial^{2} f}{\partial x^{2}} +
%                \frac{\partial^{2} f}{\partial y^{2}}$$
%
% The Laplacian operator is implemented in OpenCV by the function
% |cv.Laplacian|. In fact, since the Laplacian uses the gradient of images, it
% calls internally the _Sobel_ operator to perform its computation.
%

%% Code
%
% The program:
%
% * Loads an image
% * Remove noise by applying a Gaussian blur and then convert the original
%   image to grayscale
% * Applies a Laplacian operator to the grayscale image and stores the output
%   image
% * Display the result in a window
%

%%
% load source image
src = cv.imread(fullfile(mexopencv.root(),'test','butterfly.jpg'), 'Color',true);

%%
% apply a Gaussian blur to reduce the noise
src = cv.GaussianBlur(src, 'KSize',[3 3]);

%%
% convert filtered image to grayscale
gray = cv.cvtColor(src, 'RGB2GRAY');

%%
% apply the Laplacian operator to the grayscale image
% (input is 8-bit, we set the output image depth to 16-bit to avoid overflow)
dst = cv.Laplacian(gray, 'KSize',3, 'DDepth','int16');

%%
% take absolute value and convert results back to 8-bit
dstabs = cv.convertScaleAbs(dst);

%%
% show result
imshow(dstabs)
title('Laplace Demo')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [A0, C0] = calc_dc_components(ai)

% Calculate DC components.
% A0 and C0 are bias coefficeis, corresponding to a frequency of zero.

    %% Maximum length of chain code
    k = size(ai, 2);

    %% Traversal time and distance
    t = calc_traversal_time(ai);
    s = calc_traversal_dist(ai);

    %% Basic period of the chain code
    T = t(k);

    %% DC Components: A0, C0
    sum_a0 = 0;
    sum_c0 = 0;

    for p = 1 : k

        delta_d = calc_traversal_dist(ai(p));
        delta_x = delta_d(:,1);

        delta_y = delta_d(:,2);
        delta_t = calc_traversal_time(ai(p));

        if (p > 2)
                zeta = s(p - 1, 1) - delta_x / delta_t * t(p - 1);
                delta = s(p - 1, 2) - delta_y / delta_t * t(p - 1);
        else
                zeta = 0;
                delta = 0;
        end

        if (p > 2)
            sum_a0 = sum_a0 + delta_x / (2 * delta_t) * ((t(p))^2 - (t(p - 1))^2) + zeta * (t(p) - t(p-1));
            sum_c0 = sum_c0 + delta_y / (2 * delta_t) * ((t(p))^2 - (t(p - 1))^2) + delta * (t(p) - t(p-1));
        else
            sum_a0 = sum_a0 + delta_x / (2 * delta_t) * (t(p))^2 + zeta * t(p);
            sum_c0 = sum_c0 + delta_y / (2 * delta_t) * (t(p))^2 + delta * t(p);
        end

    end

    %% Assign  to output
    A0 = sum_a0 / T;
    C0 = sum_c0 / T;
end


text:
function fb = acc_nav2body (acc_n, DCMnb_m)
% acc_nav2body: transforms accelerations from navigation frame to body
% frame.
%
% INPUT
%	acc_n: Nx3 matrix with [fn, fe, fd] accelerations in the navigation
%		frame.
%   DCMnb_m: Nx9 matrix with nav-to-body direct cosine matrices (DCM).
%       Each row of DCMnb_m contains the 9 elements of a particular DCMnb
%       matrix ordered as [a11 a21 a31 a12 a22 a32 a13 a23 a33].
%
% OUTPUT
%	fb: Nx3 matrix with [fx, fy, fz] simulated accelerations in the
%		body frame.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%   R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics}, vol. 17,
% issue 2, pp. 110-120, 2015. Eq. 10.
%
% Version: 002
% Date:    2020/11/03
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

fb = zeros(size(acc_n));

M = max(size(acc_n));

for k = 1:M

    dcmnb = reshape(DCMnb_m(k,:), 3, 3);
    fb(k,:) = ( dcmnb * ( acc_n(k,:)' ) )';
end

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo