text:
function [featureshape,trials,vectorize_trials] = utl_determine_featureshape(trials,shape,multitask)
% Uniformize the given trials and shape information.
% [FeatureShape,Trials,VectorizeTrials] = utl_determine_featureshape(Trials,Shape,Multitask)
%
% This function deals with the fact that trials can be represented either as [NxF] matrix of
% vectorized features (N=#observations, F=#features), or as [AxBxCx...xN] array of tensor-shaped
% features (A,B,C, ... = #elements in the respective dimension), which is a bit inconsistent for
% historical reasons.
%
% Takes trials in a variety of shapes, and an (optionally non-empty) desired shape vector and
% uniformizes them into a [#trials x #features] matrix of vectorized trials, the deduced feature
% shape, and a boolean flag of whether trials had to be vectorized. By default also handles a cell
% array of multiple trial arrays, which are assumed to stem from multiple compatibly shaped tasks.
%
% In:
%   Trials : A NxF matrix of observations with pre-vectorized features, or a [AxBxCx...xN] array
%            of tensor-shaped features. Can also be a cell array of multiple such arrays where all
%            dimensions except N must be identical.
%
%   Shape : A vector that is the desired output of size() applied to a single observation, that is,
%           a vector of sizes. If [], shape information will be deduced from Trials, otherwise
%           overrides the shape of Trials.
%
%   Multitask : If true, the output will be a cell array of trial matrices. Otherwise it will be
%               a matrix (and an error will be thrown if there was more than one task in the
%               input).
%
% Out:
%   FeatureShape : Final 1xD vector that describes the feature shape (same as size(one_trial)).
%
%   Trials : Vectorized version of input Trials (#trials x #features).
%
%   VectorizeTrials : whether trials had to be vectorized.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-02-04

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

if ~iscell(trials)
    trials = {trials}; end
featureshape = cell(1,length(trials));
vectorize_trials = cell(1,length(trials));
% for each task (each of which has multiple trials)...
for t=1:length(trials)
    if ndims(trials{t}) > 2 %#ok<ISMAT>
        featureshape{t} = size(trials{t}); featureshape{t} = featureshape{t}(1:end-1);
        if ~isempty(shape) && ~isequal(shape,featureshape{t})
            if prod(featureshape{t}) == prod(shape)
                warning('You are specifying a shape property but also multidimensional features of a different shape; using the explicit shape parameter.');
                featureshape{t} = shape;
            else
                error('You are specifying a shape property but also features with an incompatible number of elements. Please correct.');
            end
        end
        trials{t} = double(reshape(trials{t},[],size(trials{t},ndims(trials{t})))');
        vectorize_trials{t} = true;
    else
        if ~isempty(shape)
            if size(shape,1) > 1
                if all(all(bsxfun(@eq,shape(1,:),shape)))
                    featureshape{t} = [shape(1,1),shape(1,2),size(shape,1)];
                    if prod(featureshape{t}) == size(trials{t},2)
                        % the reason is that it is much more efficient to operate on a dense 3d array than a very sparse 2d array
                        warn_once('This method will by convention reshape block-diagonalized feature matrices with identical blocks into a 3d tensor. This warning will only come up once.');
                    else
                        error('Your shape parameter has a different number of features than your data.');
                    end
                else
                    % we don't implement block-diagonalization in here
                    error('This method does not handle implicitly block-diagonal features; please either reformulate in tensor form or pass a large sparse data matrix (pre-blockdiagonalized). Note that the tensor form is likely several times faster.');
                end
            elseif prod(shape) ~= size(trials{t},2)
                error('Your shape parameter has a different number of features than data.');
            else
                featureshape{t} = shape;
            end
        else
            featureshape{t} = [size(trials{t},2),1];
        end
        vectorize_trials{t} = false;
    end
end
vectorize_trials = unique([vectorize_trials{:}]);
if length(vectorize_trials)>1 || ~all(cellfun(@(x)isequal(x,featureshape{1}),featureshape))
    error('The number or shape of features must be the same for each task.'); end
featureshape = featureshape{1};
if nargin >=3 && ~multitask
    if length(trials) > 1
        error('Multi-task learning is disabled by multiple tasks were given.');
    else
        trials = trials{1};
    end
end

repo: REST


text:
function Ximg = preimage_rbf(Xtr,sig2,U,B,type,npcs,maxIts)
%
% function Ximg = preimage_rbf(Xtr,sig2,U,B,type,npcs,maxIts)
%   Reconstruction or denoising after kernel PCA with RBF kernels, i.e. to find the
%   approximate preimage (in the input space) of the corresponding feature space expansions
%
% Inputs
%   Xtr     : N by d matrix of the training data used to find the prinicipal components.
%   sig2    : parameter for the RBF kernel used, k(x,z)=exp(-norm(x-z)^2/sig2).
%   U       : the eigenvectors computed from the kernel PCA using RBF kernel with parameter sig2.
%   B       : for reconstruction, B is the compressed data,
%               i.e. the projections of the data on to the first n PCs;
%             for denoising, B is the Nt by d matrix of original noisy data;
%             if not specified, Xtr is denoised instead.
%   type    : 'reconstruct' or 'denoise'
%   npcs    : number of PCs used for approximation
%   maxIts  : maximum iterations allowed to update the preimage, 1000 by default.
%
% Outputs
%   Ximg    : the reconstructed or denoised data in the input space
%
% Usage e.g.
%   >> [lam,U] = kpca(Xtr,'RBF_kernel',sig2);
%   >> [lam, perm] = sort(-lam); lam = -lam; U = U(:,perm);
%   >> projections = kernel_matrix(Xtr,'RBF_kernel',sig2,Xtest)'*U;
%   >> Xr = preimage_rbf(Xtr,sig2,U,projections(:,1:npcs),'r'); % Reconstruction
%   >> Xd = preimage_rbf(Xtr,sig2,U(:,1:npcs),Xnoisy,'d');      % Denoising
%   >> Xdtr = preimage_rbf(Xtr,sig2,U(:,1:npcs));  % Denoising on the training data
%
% see also:
%    kpca, denoise_kpca, RBF_kernel
%
% Reference
%   Mika S., Schoelkopf B., Smola A., Muller K.-R., Scholz M., Ratsch G. (1999), ``Kernel
%   PCA and de-noising in feature spaces'', Advances in Neural Information Processing
%   Systems 11, 536-542, MIT Press.
%

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

MAXDX=1e-6; % Convergence criterion
[~, dim]=size(Xtr);
mX=mean(Xtr); sX=std(Xtr);

if nargin<4, B = Xtr; end
[Nt, dimB] = size(B);

if nargin<5
    if dimB==dim, type='denoise'; else type='reconstruct'; end
else
    if type(1)~='d'&type(1)~='r'| (type(1)=='d'&dimB~=dim),
        warning('Invalid type specified, default value is used!')
        if dimB==dim, type='denoise'; else type='reconstruct'; end
    end
end

if nargin<6
    if type(1)=='r'; npcs=dimB; else npcs=size(U,2); end
else
    if npcs>size(U,2) | (type(1)=='r' & npcs~=dimB),
        warning('Invalid number of PCs, default value is used!'),
        if type(1)=='r'; npcs=dimB; else npcs=size(U,2); end
    end
end

if nargin<7, maxIts=1000; end

U=U(:,1:npcs);

for n=1:Nt
    cont=1; t=0; ts=0;
    if type(1)=='r'
        % reconstuction
        rs = U*B(n,:)';
        x = zeros(1, dim); % set the initial value of approximate preimage for reconstruction
        k = RBF_kernel(x,Xtr,sig2);
    else
        % denoise
        x = B(n,:); % initial value of the approximate preimage for denoising is the noisy data
        k = RBF_kernel(x,Xtr,sig2);
        rs = U*(k'*U)';
    end

    % iteratively update the approximate preimage x
    %
    while cont,

        d=rs'*k; % the reconstruction error is (const-2d)

        if d==0;
            % choose a different starting value
%            [k,id]=min(k); x_new = Xtr(id,:);
            randn('state',cputime+ts); x_new=(randn(1,dim)+mX).*sX;
            fprintf('%5d> Starting value changed!(d=0) \n', ts);
        else
            % update approximate preimage with a linear combination of kpca training data
            x_new = sum(rs.*k*ones(1,dim).*Xtr)/d;
        end

        dx = norm(x_new - x);
        x = x_new;

        t=t+1; ts=ts+1;

        if dx<MAXDX;
            cont = 0;
         % fprintf('%5d> Converged! \n', ts);
        else
            if ts>=maxIts;
                cont=0;
                fprintf('%5d> Maximum iteration reached!\n', ts);
            elseif t>=500;
                % choose a different starting value
%                [k,id]=min(k); x = Xtr(id,:);
                randn('state',cputime+ts); x=(randn(1,dim)+mX).*sX;
                t=0;
            end
        end
        k = RBF_kernel(x,Xtr,sig2);

    end % while
    Ximg(n,:)=x_new;
end % for

author: peterhcharlton

repo: RRest


text:
function ind = ichoose(n,k)
%ICHOOSE Choose K indices from the set 1:N.
%   ICHOOSE(N,K) gives all combinations of choosing K indices from
%   the set 1:N without order and without repetitions.
%
%   Example:
%   	ind = ichoose(4,2);
%   	v = 'ABCD';
%   	v(ind) % ['AB';'AC';'AD';'BC';'BD';'CD']
%
%   See also NCHOOSEK

%   Author: Jonas Lundgren <splinefit@gmail.com> 2011

% Check input
if ~isscalar(n) || ~isfinite(n) || n < 0 || fix(n) < n
    error('ICHOOSE:N','N must be a non-negative integer.')
elseif ~isscalar(k) || ~isfinite(k) || k < 0 || fix(k) < k
    error('ICHOOSE:K','K must be a non-negative integer.')
elseif n < k
    ind = zeros(0,k);
    return
elseif k == 0
    ind = zeros(1,0);
    return
end

% Initiate indices
if n < 256
    ind = uint8(k:n)';
elseif n < 65536
    ind = uint16(k:n)';
else
    ind = (k:n)';
end

m = n-k+1;
b = 1:m;

% Generate indices
for j = 1:k-1
    a = b;
    b = cumsum(b);
    c = b(m) + 1 - b(1:m-1);
    if m > 1
        p = ones(b(m),1);
        p(c) = 1 - a(1:m-1);
        p = cumsum(p);
        ind = ind(p,:);
    end
    p = zeros(b(m),1);
    p(1) = k-j;
    p(c) = 1;
    p = cumsum(p);
    ind = [p,ind]; %#ok
end

author: Sable

repo: mcbench-benchmarks


text:
function [C,h,Ph,F] = spm_ar_reml(YY,X,m,N)
% ReML estimation of covariance components from y*y'
% FORMAT [C,h,Ph,F] = spm_ar_reml(YY,X,m,N);
%
% YY  - (m x m) sample covariance matrix Y*Y'  {Y = (m x N) data matrix}
% X   - (m x p) design matrix
% m   - (1) order of AR(m) model
% N   - number of samples
%
%
% C   - (m x m) estimated errors = h(1)*Q{1} + h(2)*Q{2} + ...
% h   - (q x 1) ReML hyperparameters h: normalised AR coeficients
% Ph  - (q x q) conditional precision of h (unnormalised)
%
% F   - [-ve] free energy F = log evidence = p(Y|X,Q) = ReML objective
%
% Performs a Fisher-Scoring ascent on F to find ReML variance parameter
% estimates.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Karl Friston
% $Id: spm_ar_reml.m 5219 2013-01-29 17:07:07Z spm $

% assume a single sample if not specified
%--------------------------------------------------------------------------
try
    N;
catch
    N  = 1;
end

% assume AR(1) if not specified
%--------------------------------------------------------------------------
try
    m;
catch
    m  = 1;
end

% ortho-normalise X
%--------------------------------------------------------------------------
if isempty(X)
    X = sparse(length(YY),1);
else
    X = orth(full(X));
end

% initialise h
%--------------------------------------------------------------------------
m     = m + 1;
n     = length(YY);
dh    = zeros(m,1);
dFdh  = zeros(m,1);
dFdhh = zeros(m,m);
L     = zeros(m,m);

% initialise and specify hyperpriors
%--------------------------------------------------------------------------
hE  = sparse(1,1,1,m,1);
hP  = speye(m,m)/exp(32);
h   = hE;

% initialise precision components
%--------------------------------------------------------------------------
for i = 1:m
    Q{i} = spdiags(ones(n,2),[(1 - i) (i - 1)],n,n);
end

% scale data
%--------------------------------------------------------------------------
Ys    = norm(YY,1)/N;
YY    = YY/Ys;

% ReML (EM/VB)
%--------------------------------------------------------------------------
for k = 1:64

    % compute current estimate of covariance
    %----------------------------------------------------------------------
    iC    = sparse(n,n);
    for i = 1:m
            iC = iC + Q{i}*h(i);
    end
    C     = inv(iC);

    % E-step: conditional covariance cov(B|y) {Cq}
    %======================================================================
    Cq    = pinv(X'*iC*X);

    % M-step: ReML estimate of hyperparameters
    %======================================================================

    % Gradient dF/dh (first derivatives)
    %----------------------------------------------------------------------
    P     = C - X*Cq*X';
    U     = (iC*YY/N*iC - iC)*P;
    for i = 1:m

        % dF/dh = -trace(dF/diC*iC*Q{i}*iC)
        %------------------------------------------------------------------
        QP{i}     = Q{i}*P;
        dFdh(i)   = -trace(QP{i}*U)*N/2;

    end

    % Expected curvature E{dF/dhh} (second derivatives)
    %----------------------------------------------------------------------
    for i = 1:m
        for j = i:m

            % dF/dhh = -trace{P*Q{i}*P*Q{j}}
            %--------------------------------------------------------------
            dFdhh(i,j) = -trace(QP{i}*QP{j})*N/2;
            dFdhh(j,i) =  dFdhh(i,j);

        end
    end

    % add hyperpriors
    %----------------------------------------------------------------------
    e     = h     - hE;
    dFdh  = dFdh  - hP*e;
    dFdhh = dFdhh - hP;

    % update regulariser
    %----------------------------------------------------------------------
    if ~rem(k,8)
       L  = speye(m,m)*norm(dFdhh)/128;
    end

    % Fisher scoring: update dh = -inv(ddF/dhh)*dF/dh
    %----------------------------------------------------------------------
    Ph    = -dFdhh;
    dh    = -pinv(dFdhh - L)*dFdh;

    % preclude numerical overflow
    %----------------------------------------------------------------------
    h     = h + dh;

    % Convergence (1% change in log-evidence)
    %======================================================================
    dF    = dFdh'*dh;
    fprintf('%-30s: %i %30s%e\n','  ReML Iteration',k,'...',full(dF));
    if dF < 1e-1, break, end

end

% log evidence = ln p(y|X,Q) = ReML objective = F = trace(R'*iC*R*YY)/2 ...
%--------------------------------------------------------------------------
if nargout > 3
    R = P*iC;
    F = - trace(R*YY*R')/2 ...
        - e'*hP*e/2 ...
        - N*n*log(2*pi)/2 ...
        - N*spm_logdet(C)/2 ...
        + N*spm_logdet(Cq)/2 ...
        -   spm_logdet(Ph)/2 ...
        +   spm_logdet(hP)/2;
end

% rescale (NB - Q{1) = 2*I
%--------------------------------------------------------------------------
C    =  C*Ys;
h    = -h(2:m)/(h(1)*2);

author: spm


text:
function r = loggmpdf(X, model)
% Compute log pdf of a Gaussian mixture model.
% Written by Mo Chen (sth4nth@gmail.com).
mu = model.mu;
Sigma = model.Sigma;
w = model.weight;

n = size(X,2);
k = size(mu,2);
logRho = zeros(k,n);

for i = 1:k
    logRho(i,:) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
end
r = logsumexp(bsxfun(@plus,logRho,log(w)'),1);

function y = loggausspdf(X, mu, Sigma)
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= chol(Sigma);
if p ~= 0
    error('ERROR: Sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U)));   % normalization constant
y = -(c+q)/2;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


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

    enumeration
        Minimize('Minimize')
        Maximize('Maximize')
    end

    properties
        name(1,:) char
    end

    methods
        function obj = ObjFcnDirectionTypeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('ObjFcnDirectionTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ObjFcnDirectionTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ObjFcnDirectionTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% Some parameters to set - make sure that your code works at image borders!
close all;
patchSize = 5;
sigma = 20; % standard deviation (different for each image!)
h = 0.55; %decay parameter
windowSize = 15; % --> it takes 99 seconds to run

%NOTE: for each image, please also read its CORRESPONDING 'clean' or
%reference image. We will need this later to do some analysis
%NOTE2: the noise level is different for each image (it is 20, 10, and 5 as
%indicated in the image file names)
%NOTE3 : all the functions to denoise the image work both with grayscale
%images and with RGB images

%Load the image to denoise and the refernce image
imageNoisy = imread('images/alleyNoisy_sigma20.png');
imageReference = imread('images/alleyReference.png');

% converting the image in grayscale is much faster roughly 3 times faster
imageNoisy = rgb2gray(imageNoisy);
imageReference = rgb2gray(imageReference);

% we can resize the image so the code can run even faster
%imageNoisy = imresize(imageNoisy, 0.5);
%imageReference = imresize(imageReference, 0.5);

% Implement the non-local means functions and compare the performances

%############################## SIGMA = 20 ######################
% tic
% filtered = nonLocalMeansWithoutIntegral(imageNoisy, sigma, h, patchSize, windowSize, false);
% toc
%
% figure('name', 'NL-Means Denoised Image with Naive Approach -- SIGMA = 20 h = 0.55');
% imshow(filtered);

tic
filtered = nonLocalMeans(imageNoisy, sigma, h, patchSize, windowSize, false);
toc

figure('name', 'NL-Means Denoised Image with Integral Images Approach -- SIGMA = 20 h = 1');
imshow(filtered);

h = 0.15;

% tic
% filtered = nonLocalMeansWithoutIntegral(imageNoisy, sigma, h, patchSize, windowSize, false);
% toc
%
% figure('name', 'NL-Means Denoised Image with Naive Approach -- SIGMA = 20 h = 0.25');
% imshow(filtered);

tic
filtered = nonLocalMeans(imageNoisy, sigma, h, patchSize, windowSize, false);
toc

figure('name', 'NL-Means Denoised Image with Integral Images Approach -- SIGMA = 20 h = 0.25');
imshow(filtered);

%############################### SIGMA = 5 ####################

imageNoisy = imread('images/townNoisy_sigma5.png');
imageReference = imread('images/townReference.png');

% converting the image in grayscale is much faster roughly 3 times faster
imageNoisy = rgb2gray(imageNoisy);
imageReference = rgb2gray(imageReference);

sigma = 5; % standard deviation (different for each image!)
h = 0.55; %decay parameter

% tic
% filtered = nonLocalMeansWithoutIntegral(imageNoisy, sigma, h, patchSize, windowSize, false);
% toc
%
% figure('name', 'NL-Means Denoised Image with Naive Approach -- SIGMA = 5 h = 0.55');
% imshow(filtered);

tic
filtered = nonLocalMeans(imageNoisy, sigma, h, patchSize, windowSize, false);
toc

figure('name', 'NL-Means Denoised Image with Integral Images Approach -- SIGMA = 5 h = 0.55');
imshow(filtered);

h = 0.25;

% tic
% filtered = nonLocalMeansWithoutIntegral(imageNoisy, sigma, h, patchSize, windowSize, false);
% toc
%
% figure('name', 'NL-Means Denoised Image with Naive Approach -- SIGMA = 5 h = 0.25');
% imshow(filtered);

tic
filtered = nonLocalMeans(imageNoisy, sigma, h, patchSize, windowSize, false);
toc

figure('name', 'NL-Means Denoised Image with Integral Images Approach -- SIGMA = 5 h = 0.25');
imshow(filtered);

%% Let's show your results!

%Show the denoised image
% figure('name', 'NL-Means Denoised Image');
% imshow(filtered);

%Show difference image
diff_image = abs(imageReference - filtered);
figure('name', 'Difference Image');
imshow(diff_image ./ max(max((diff_image))));

%Print some statistics ((Peak) Signal-To-Noise Ratio)
disp('For Noisy Input');
[peakSNR, SNR] = psnr(imageNoisy, imageReference);
disp(['SNR: ', num2str(SNR, 10), '; PSNR: ', num2str(peakSNR, 10)]);

disp('For Denoised Result');
[peakSNR, SNR] = psnr(filtered, imageReference);
disp(['SNR: ', num2str(SNR, 10), '; PSNR: ', num2str(peakSNR, 10)]);

%Feel free (if you like only :)) to use some other metrics (Root
%Mean-Square Error (RMSE), Structural Similarity Index (SSI) etc.)

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [x,zo]=estnoiseg_orig(yf,tz,pp)
%ESTNOISEG - estimate MMSE noise spectrum [x,zo]=(yf,tz,pp)
%
% Usage:    ninc=round(0.016*fs);   % frame increment [fs=sample frequency]
%           ovf=2;                  % overlap factor
%           f=rfft(enframe(s,hanning(ovf*ninc,'periodic'),ninc),ovf*ninc,2);
%           f=f.*conj(f);           % convert to power spectrum
%           x=estnoiseg(f,ninc/fs); % estimate the noise power spectrum
%
% Inputs:
%   yf      input power spectra (one row per frame)
%   tz      frame increment in seconds
%           Alternatively, the input state from a previous call (see below)
%   pp      algorithm parameters [optional]
%
% Outputs:
%   x       estimated noise power spectra (one row per frame)
%   zo      output state
%
% The algorithm parameters are defined in reference [1] from which equation
% numbers are given in parentheses. They are as follows:
%
%        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 estnoiseg in chunks of arbitrary size. Thus the following are equivalent:
%
%                   (a) dp=estnoiseg(yp(1:300),tinc);
%
%                   (b) [dp(1:100,:),z]=estnoiseg(yp(1:100,:),tinc);
%                       [dp(101:200,:),z]=estnoiseg(yp(101:200,:),z);
%                       [dp(201:300,:),z]=estnoiseg(yp(201:300,:),z);

% This is intended to be a precise implementation of [1] for a frame rate of 62.5 Hz.
% Time constants are adjusted for other frame rates.
%
% Refs:
%    [1] 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

%	   Copyright (C) Mike Brookes 2012
%      Version: $Id: estnoiseg.m 3387 2013-08-23 12:32:47Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nr,nrf]=size(yf);          % number of frames and freq bins
x=zeros(nr,nrf);            % initialize output arrays
if isempty(yf) && isstruct(tz)             % no real data
    zo=tz;                  % just keep the same state
else
    if isstruct(tz)         % take parameters from a previous call
        nrcum=tz.nrcum;     % cumulative number of frames
        xt=tz.xt;           % smoothed power spectrum
        pslp=tz.pslp;       % correction factor (9)
        tinc=tz.tinc;       % frame increment
        qq=tz.qq;           % parameter structure
    else
        tinc = tz;          % second argument is frame increment
        nrcum=0;            % no frames so far
        % default algorithm constants
        qq.tax=0.0717;      % noise output smoothing time constant = -tinc/log(0.8) (8)
        qq.tap=0.152;       % speech prob smoothing time constant = -tinc/log(0.9) (23)
        qq.psthr=0.99;      % threshold for smoothed speech probability [0.99] (24)
        qq.pnsaf=0.01;      % noise probability safety value [0.01] (24)
        qq.pspri=0.5;       % prior speech probability [0.5] (18)
        qq.asnr=15;         % active SNR in dB [15] (18)
        qq.psini=0.5;       % initial speech probability [0.5] (23)
        qq.tavini=0.064;        % assumed speech absent time at start [64 ms]

        if nargin>=3 && ~isempty(pp)  % update fields from pp input
            qqn=fieldnames(qq);
            for i=1:length(qqn)
                if isfield(pp,qqn{i})
                    qq.(qqn{i})=pp.(qqn{i});
                end
            end
        end
        pslp=repmat(qq.psini,1,nrf); % initialize smoothed speech presence prob
        xt=[];                       % initialize just in case the first call has no data
    end

    % unpack parameters needed within the loop

    psthr=qq.psthr;     % threshold for smoothed speech probability [0.99] (24)
    pnsaf=qq.pnsaf;     % noise probability safety value [0.01] (24)

    % derived algorithm constants

    ax=exp(-tinc/qq.tax); % noise output smoothing factor = 0.8 (8)
    axc=1-ax;
    ap=exp(-tinc/qq.tap); % noise output smoothing factor = 0.9 (23)
    apc=1-ap;
    xih1=10^(qq.asnr/10); % speech-present SNR
    xih1r=1/(1+xih1)-1;
    pfac=(1/qq.pspri-1)*(1+xih1); % p(noise)/p(speech) (18)

    if nrcum==0 && nr>0       % initialize values for first frame
        xt=qq.psini*mean(yf(1:max(1,min(nr,round(1+qq.tavini/tinc))),:),1);  % initial noise estimate
    end

    % loop for each frame
    for t=1:nr
        yft=yf(t,:);        % noisy speech power spectrum
        ph1y=(1+pfac*exp(xih1r*yft./xt)).^(-1); % a-posteriori speech presence prob (18)
        pslp=ap*pslp+apc*ph1y; % smoothed speech presence prob (23)
        ph1y=min(ph1y,1-pnsaf*(pslp>psthr)); % limit ph1y (24)
        xtr=(1-ph1y).*yft+ph1y.*xt; % estimated raw noise spectrum (22)
        xt=ax*xt+axc*xtr;  % smooth the noise estimate (8)
        x(t,:)=xt;  % save the noise estimate
    end
    if nargout>1    % we need to store the state for next time
        zo.nrcum=nrcum+nr;      % number of frames so far
        zo.xt=xt;          % smoothed power spectrum
        zo.pslp=pslp;               % correction factor (9)
        zo.tinc=tinc;     % must be the last one
        zo.qq=qq;
    end
    if ~nargout
        clf;
        subplot(212);
        plot((1:nr)*tinc,10*log10([sum(yf,2) sum(x,2)]))
        ylabel('Frame Energy (dB)');
        xlabel(sprintf('Time (s)   [%d ms frame incr]',round(tinc*1000)));
        axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
        subplot(211);
        plot(1:nrf,10*log10([sum(yf,1)'/nr sum(x,1)'/nr]))
        ylabel('Power (dB)');
        xlabel('Frequency bin');
        axisenlarge([-1 -1.05]);
        legend('input','noise','Location','Best');
    end
end

author: LimingShi


text:
function Pset=str2poly(strcell,varnames)
% POLY2STR - convert multivariate polynomials from strings to arrays
% usage: Pset=poly2str(strcell,varnames)
%
% INPUTS: strcell, varnames
%  strcell is a cell array of strings, one per polynomial
%    all numbers must be whole or decimal form (no scientific notation)
%    brackets and double negatives are not allowed.
%    All spaces in the strings are ignored.
%  varnames (optional) is a cell array of strings which must include all
%    variable names that occur in strcell [default: {'x1','x2',...}]
%
% OUTPUTS: Pset
%  Pset is a cell array of polynomial coefficients where row j of Pset{i}
%    is [c,k1,k2,...] representing the monomial c*x1^k1*x2^k2*...
%
% EXAMPLE:
%  str2poly({'x1*x2^2'})    % returns {[1,1,2]}
%  str2poly({'-x1*x2^2'})   % returns {[-1,1,2]}
%  str2poly({'-1*x1*x2^2'}) % returns {[-1,1,2]}
%  str2poly({'- 1 * x 1 * x 2 ^ 2'}) % returns {[-1,1,2]}
%
% SEE ALSO:
%  groebner, poly2str

% Author: Ben Petschel 19/6/2009
%
% Change history:
%  19/6/2009 - first release
%  22/6/2009 - fixed bug in handling "-x1" etc.
%  17/7/2009 - allow whitespace in strings
%  20/3/2010 - changed array representation to rectangular instead of n-d

if nargin<2,
  varnames = {}; % will fill in varnames later, as needed
end;
Pset = cell(size(strcell));

for i=1:numel(strcell),
  str = strcell{i};

  % ignore all spaces
  remain = str;
  remain = strrep(remain,' ','');

  % first insert "+" before every "-", for easier term separation later
  str = '';
  while ~isempty(remain),
    if remain(1)=='-',
      if isempty(str),
        str = '-';
      else
        str = [str,'+-'];
      end;
    end;
    [tok,remain]=strtok(remain,'-');
    str = [str,tok];
  end;

  % now separate into terms
  remain = str;
  terms = [];
  while ~isempty(remain),
    coeff = 1;
    exps = [];
    [term,remain]=strtok(remain,'+');
    % process term, separating into products
    rem2 = term;
    while ~isempty(rem2),
      [part,rem2]=strtok(rem2,'*');
      if ~isnan(str2double(part)),
        % part is a number, multiply it by the coefficient
        coeff = coeff * str2double(part);
      else
        % part is var^n or var or -var^n or -var
        if part(1)=='-',
          coeff = -coeff;
          part=part(2:end);
          if isempty(part),
            error('isolated minus sign');
          end;
        end;
        ind = find(part=='^');
        if isempty(ind),
          var = part;
          expo = 1;
        else
          ind = ind(1);
          var = part(1:ind-1);
          expo = str2double(part(ind+1:end));
          if isnan(expo),
            % exponent is not valid
            error('invalid exponent in part %s',part);
          end;
        end;
        if isempty(varnames),
          % see if var is of form 'x1' etc
          if (length(var)>1) && (var(1)=='x') && ~isnan(str2double(var(2:end))),
            % variable name is ok
            k = str2double(var(2:end));
          else
            error('variable name not of form "x1","x2",...');
          end;
        else
          ind = strmatch(var,varnames);
          if isempty(ind),
            error('variable name not in list');
          end;
          k = ind(1);
        end;
        if length(exps)<k,
          exps = [exps,zeros(1,k-length(exps))];
        end;
        exps(k) = exps(k)+expo;
      end;
    end; % while ~isempty(rem2),
    % term is processed, now add it to the list of terms
    if size(terms,2)<(length(exps)+1),
      terms = [terms,zeros(size(terms,1),length(exps)+1-size(terms,2))];
    elseif length(exps)<size(terms,2)-1,
      exps = [exps,zeros(1,size(terms,2)-1-length(exps))];
    end;
    terms = [terms;coeff,exps];
  end; % while ~isempty(remain),

  Pset{i}=terms;
end; % for i=1:numel(strcell),

end % main function str2poly(...)


text:
function num_cores( desiredCores,flag_parpool )
% for implicit parallel processing
%
% (c) Marc Fischer, Thomas Kuestner
% ---------------------------------------------------------------------

availableCores = maxNumCompThreads;
if availableCores <= desiredCores
    maxNumCompThreads(availableCores);
else
    maxNumCompThreads(desiredCores);
end;

%for explicit parallel processing:
if flag_parpool
    c = parcluster();
    p = gcp('nocreate');
    if isempty(p)
        p = parpool(c,desiredCores-1); % -1 for work on local pc
    end;
end;

end

author: thomaskuestner


text:
clear all;
%%%====== Settings ======%%%
model = 'Deep SR-ITM'; % Deep SR-ITM (ICCV'19) or Multi-purpose CNN (ACCV'18)
yuv_format = '420'; % YUV file format
SDR_file = './data/test/testset_SDR.yuv'; % input .yuv file
wid = 1920; % width
hei = 1080; % height
num_fr = 28; % number of frames in the YUV file
scale = 2; % scale factor for SR
pred_file = sprintf('./pred/pred_x%d.yuv', scale); % result .yuv file
gpuDevice(1); % GPU ID
%%%======================%%%
addpath('utils');
disp(['Testing for scale ', num2str(scale), '...'])
% initialize
fclose(fopen(pred_file, 'w'));
[fwidth,fheight] = yuv_factor(yuv_format);

% load net
disp('Loading net...')
if strcmp(model, 'Deep SR-ITM')
    netstruct = load(sprintf('./net/x%d.mat', scale));
elseif strcmp(model, 'Multi-purpose CNN')
    netstruct = load(sprintf('./net/Multi-purpose_CNN_x%d.mat', scale));
end
net = dagnn.DagNN.loadobj(netstruct.net);
move(net,'gpu');
net.mode = 'test' ;
pred_index = net.getVarIndex('pred');
net.conserveMemory = true;

% test
disp('Testing starts...')
for fr = 1:num_fr
    % read frames
    SDR_YUV = uint8(load_yuv(SDR_file, fr, hei, wid, fheight, fwidth, 'SDR'));
    % normalize
    SDR_YUV = single(SDR_YUV)/255;
    % change type
    SDR_YUV = gpuArray(SDR_YUV);
    % prediction
    net.eval({'input', SDR_YUV});
    pred = gather(net.vars(pred_index).value);
    pred = min(max(pred, 0), 1);
    % save yuv file
    save_yuv(pred*1023, pred_file, hei, wid, fheight, fwidth, 'HDR');
    disp('Frame saved!')
end
disp('Done!')

author: sooyekim

repo: Deep-SR-ITM


text:
function [mod_res_apx exp_res_apx]=exp_approx2(dboun, exp_par1, exp_par2,show_res)
nexp=size(exp_par1,1);
exp_res_apx=zeros(nexp,3);
mod_res_apx=zeros(nexp,3);

for i=1:nexp
    dind=(dboun(i,1):dboun(i,2));

    Y=(exp_par1(i,1).^dind*exp_par1(i,2)-exp_par2(i,1).^dind*exp_par2(i,2))';
    p1=exp_par1(i,1);
    p2=exp_par2(i,1);
    sr_a=p1/(1-p1^2)/(p1-p2)/(1-p1*p2);
    sr_b=p2/(1-p2^2)/(p2-p1)/(1-p1*p2);
    H=[sr_a*p1.^dind'+sr_b*p2.^dind'];

    N_est=inv(H'*H)*H'*Y;

    exp_res_apx(i,:)=[p1,p2,sqrt(N_est)];
    mod_res_apx(i,:)=[p1+p2,-p1*p2,sqrt(N_est)];
end

if (show_res==1)
    for i=1:nexp
        dind=0:dboun(i,2);
        val_a=exp_par1(i,1).^dind*exp_par1(i,2);
        val_b=exp_par2(i,1).^dind*exp_par2(i,2);
        p1=exp_res_apx(i,1);
        p2=exp_res_apx(i,2);
        N=exp_res_apx(i,3)^2;
        sr_a=N*p1/(1-p1^2)/(p1-p2)/(1-p1*p2);
        sr_b=N*p2/(1-p2^2)/(p2-p1)/(1-p1*p2);
        val_c=sr_a*p1.^dind'+sr_b*p2.^dind';
        figure;
        plot(dind,val_a-val_b);grid;
        hold on
        plot(dind,val_c,'r');
    end
end

repo: ekfmonoslam


text:
function pde = Darcydata3
%% DARCYDATA3 Jump coefficients
%
% Use for grid with size h=1/4

pde = struct('f', @f,'g_D',@g_D,'g_N',@g_N);

% p = round(5*rand(16,1));
p = [2  3  2  3  3  4  3  5  1  2  1  0  2  1  2  4]';
pde.K = repmat(10.^(-p),2,1); % repeat for each triangle in the small square

    function s = f(pt)
    x = pt(:,1); y = pt(:,2);
    s = sin(2*pi*x).*cos(2*pi*y);
    end
    function s = g_D(pt)
    s = 0;
    end
    function s = g_N(pt,vargin)
    s = 0;
    end
end

repo: ifem


text:
function laguerre_polynomial_test08 ( p, e )

%*****************************************************************************80
%
%% LAGUERRE_POLYNOMIAL_TEST08 tests L_POWER_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polynomial
%    factors.
%
%    Input, integer E, the exponent of X.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLYNOMIAL_TEST08\n' );
  fprintf ( 1, '  Compute a power product table for L(n,x):\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tij = integral ( 0 <= x < +oo ) x^e L(i,x) L(j,x) exp(-x) dx\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum degree P = %d\n', p );
  fprintf ( 1, '  Exponent of X, E = %d\n', e );

  table = l_power_product ( p, e );

  r8mat_print ( p + 1, p + 1, table, '  Power product table:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function behaviorAnnotator( fName, aName, tName )
% Caltech Behavior Annotator.
%
% Can only be used to annotated seq files. To annotate other types of video
% files, they need to be converted to seq files first. For example, to
% convert an avi file into an seq file, use the following command:
%  seqIo([seqName],'frImgs',struct('codec','jpg'),'aviName',[aviName])
% Where [seqName] denotes the name of the seq you want to save to and
% [aviName] is the avi file. See seqIo and seqIo>frimgs for more details.
%
% Use arrow keys to play the video. Controls are as follows:
% (1) If stopped: [L] play backward, [R] play forward
% (2) If playing forward: [L] halve speed, [R] double speed, [U/D] stop
% (3) If playing backward: [L] double speed, [R] halve speed,[U/D] stop
% You can explicitly enter a frame or use the slider to jump in the video.
%
% USAGE
%  behaviorAnnotator( [fName], [aName], [tName] )
%
% INPUTS
%  fName    - optional seq file to load at start
%  aName    - optional annotation file to load or import at start
%  tName    - optional tracking or detection file to load at start
%
% OUTPUTS
%
% EXAMPLE
%  behaviorAnnotator
%
% See also behaviorData, seqIo
%
% Piotr's Image&Video Toolbox      Version 2.60
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% handles to gui objects / other globals
[hFig,menu,pTop,pMid,pBot,dispApi,A,trk] = deal([]);

% initialize all
makeLayout();
menuApi = menuMakeApi();
dispApi = dispMakeApi();
annApi  = annMakeApi();
menuApi.vidClose();

% open vid, annotation and tracking results if given
if(nargin>=1 && ~isempty(fName)), menuApi.vidOpen(fName); end
if(nargin>=2 && ~isempty(aName)), menuApi.annOpen(aName); end
if(nargin>=3 && ~isempty(tName)), menuApi.trkOpen(tName); end

  function makeLayout()
    % common properties
    name = 'Caltech Behavior Annotator';
    bg='BackgroundColor'; fg='ForegroundColor';
    fs='FontSize'; ha='HorizontalAlignment';
    units = {'Units','pixels'}; st='String'; ps='Position';

    % initial figures size / pos
    set(0,'Units','pixels');  ss = get(0,'ScreenSize');
    if( ss(3)<800 || ss(4)<600 ); error('screen too small'); end;
    figPos = [(ss(3)-600)/2 (ss(4)-500)/2 600 500];

    % create main figure
    figPrp = {'NumberTitle','off', 'Toolbar','auto', 'MenuBar','none', 'Color','k'};
    hFig = figure(figPrp{:},'Visible','off', ps,figPos, 'Name',name );
    set(hFig,'DeleteFcn',@exitProg,'ResizeFcn',@figResized);

    % mid panel
    pMid.hAx=axes(units{:},'Parent',hFig); set(pMid.hAx,'XTick',[],'YTick',[]);
    pMid.hSl=uicontrol(hFig,'Style','slider','Min',0,'Max',1,bg,'k'); imshow(0);
    pMid.hBar=uicontrol(hFig,'Style','pushbutton',units{:},'Enable','inactive');
    pMid.hAdd=uicontextmenu('Parent',hFig); pMid.hAdds=[]; pMid.hTxt=[];

    % top panel
    pnlProp = [units {bg,[.1 .1 .1],'BorderType','none'}];
    txtPrp = {'Style','text',bg,[.1 .1 .1],fg,'w',ha};
    edtPrp = {'Style','edit',bg,[.1 .1 .1],fg,'w',ha};
    pTop.h = uipanel(pnlProp{:},'Parent',hFig);
    pTop.hBehLbl=uicontrol(pTop.h,txtPrp{:},'Left',st,'behavior:');
    pTop.hBh=uicontrol(pTop.h,units{:},'Style','popupmenu',fg,'w',bg,'k',st,{''});
    pTop.hFrmLbl=uicontrol(pTop.h,txtPrp{:},'Left',st,'frame:');
    pTop.hFrmInd=uicontrol(pTop.h,edtPrp{:},'Right');
    pTop.hTmLbl=uicontrol(pTop.h,txtPrp{:},'Left',st,'time:');
    pTop.hTmVal=uicontrol(pTop.h,txtPrp{:},'Left');
    pTop.hFrmNum=uicontrol(pTop.h,txtPrp{:},'Left');
    pTop.hPlyLbl=uicontrol(pTop.h,txtPrp{:},'Left',st,'speed:');
    pTop.hPlySpd=uicontrol(pTop.h,txtPrp{:},'Left',st,'0x');

    % bottom panel
    pBot.h=uipanel(pnlProp{:},'Parent',hFig); uiPr=[pBot.h,units,'Style'];
    uiButPr=[uiPr,'pushbutton', bg,[.9 .9 .9], fs,11, 'FontWeight','bold'];
    pBot.hStrm=uicontrol(uiPr{:},'popupmenu',fg,'w',bg,'k',st,{''});
    pBot.hBh=uicontrol(uiPr{:},'popupmenu',fg,'w',bg,'k',st,{''});
    pBot.hArrL=uicontrol(uiButPr{:},fg,'k',st,'<<');
    pBot.hSl=uicontrol(uiPr{:},'slider',bg,'k');
    pBot.hBar=uicontrol(uiPr{:},'pushbutton','Enable','inactive');
    pBot.hArrR=uicontrol(uiButPr{:},fg,'k',st,'>>');
    pBot.hSetL=uicontrol(uiPr{:},'pushbutton',fg,'k',bg,[.9 .9 .9],fs,8);
    pBot.hSetR=uicontrol(uiPr{:},'pushbutton',fg,'k',bg,[.9 .9 .9],fs,8);
    pBot.hZm=uicontrol(uiPr{:},'popupmenu',fg,'w',bg,'k',...
      st,int2str2([30 60 120 250*2.^(0:5)]),'Value',5);
    pBot.hDel=uicontrol(uiButPr{:},fg,[.5 0 0],st,'X');

    % icons for set region button
    I=ones(20,18); I(9:15,10:16)=triu(ones(7));
    I(14:15,4:13)=0; I=repmat(I,[1 1 3]);
    set(pBot.hSetL,'CData',I(:,end:-1:1,:));
    set(pBot.hSetR,'CData',I);

    % set the keyPressFcn for all focusable components (except popupmenus)
    set( [hFig, pMid.hSl pBot.hArrL pBot.hSl pBot.hArrR ...
      pBot.hSetL pBot.hSetR pBot.hDel], 'keyPressFcn',@keyPress );

    % create menus
    menu.hVid = uimenu(hFig,'Label','Video');
    menu.hVidOpn = uimenu(menu.hVid,'Label','Open');
    menu.hVidsOpn = uimenu(menu.hVid,'Label','Open dual');
    menu.hVidExp = uimenu(menu.hVid,'Label','Export');
    menu.hVidCls = uimenu(menu.hVid,'Label','Close');
    menu.hVidInf = uimenu(menu.hVid,'Label','Info');
    menu.hVidAud = uimenu(menu.hVid,'Label','Audio');
    menu.hAnn = uimenu(hFig,'Label','Annotation');
    menu.hAnnNew = uimenu(menu.hAnn,'Label','New');
    menu.hAnnOpn = uimenu(menu.hAnn,'Label','Open');
    menu.hAnnCls = uimenu(menu.hAnn,'Label','Close');
    menu.hAnnSav = uimenu(menu.hAnn,'Label','Save');
    menu.hAnnCnf = uimenu(menu.hAnn,'Label','Config');
    menu.hAnnMrg = uimenu(menu.hAnn,'Label','Merge');
    menu.hTrk = uimenu(hFig,'Label','Tracking');
    menu.hTrkOpn = uimenu(menu.hTrk,'Label','Open');
    menu.hTrkCls = uimenu(menu.hTrk,'Label','Close');

    % set hFig to visible upon completion
    set(hFig,'Visible','on'); drawnow;

    function keyPress( h, evnt )
      char=int8(evnt.Character); if(isempty(char)), char=0; end;
      % arrow keys control either video, main slider, or bottom slider
      if( char>=28 && char<=31 ) % L/R/U/D = 28/29/30/31
        if( h==hFig )
          if(char>=30), flag=0; else flag=double(char-28)*2-1; end
          dispApi.setSpeedCb(flag);
        elseif( h~=pBot.hSl && h~=pMid.hSl )
          flag = mod(double(char),2)*2-1;
          annApi.setFrame(flag);
        end
        return;
      end
      % all other keys require an annotation to be loaded
      if(isempty(A)), return; end; bp=[];
      % '-'/'=' control jumps, ','/'.' control moveLeft/Right
      if(char=='-'), bp='prevBeh'; end
      if(char=='='), bp='nextBeh'; end
      if(char==','), bp='moveLeft'; end
      if(char=='.'), bp='moveRight'; end
      if(~isempty(bp)), annApi.buttonPress(bp); return; end
      % make bindings for behavior creation
      [member,type]=ismember( char, int8(A.getKeys()) );
      if(member), annApi.insertBh(type); return; end;
    end

    function figResized( h, evnt ) %#ok<INUSD>
      % aspect ratio of video
      if(isempty(dispApi)), info=[]; else info=dispApi.getInfo(); end
      if(isempty(info)), ar=4/3; else ar=info.width/info.height; end
      % enforce minimum size (min width=500+pad*2)
      pos=get(hFig,ps); pad=8; htBot=26; htSl=20; htTop=20;
      mWd=500+pad*2; mHt=500/ar+htBot+htSl+htTop+pad*2;
      persistent posPrv;
      if(pos(3)<mWd || pos(4)<mHt && ~isempty(posPrv))
        set(hFig,ps,[posPrv(1:2) mWd mHt]); figResized(); return;
      end; posPrv=pos;
      % overall layout
      wd=pos(3)-2*pad; ht=pos(4)-2*pad-htSl-htTop-htBot;
      wd=min(wd,ht*ar); ht=min(ht,wd/ar); x=(pos(3)-wd)/2; y=pad;
      set(pBot.h,ps,[x y wd htBot]); y=y+htBot;
      set(pMid.hSl,ps,[x y wd htSl]); y=y+htSl;
      if(~isempty(pMid.hTxt)), set(pMid.hTxt,ps,[wd/2 ht]); end
      set(pMid.hAx,ps,[x y wd ht]); y=y+ht;
      set(pTop.h,ps,[x y wd htTop]);
      % position stuff in top panel
      x=10;
      set(pTop.hBehLbl,ps,[x 3 50 14]); x=x+50;
      set(pTop.hBh,ps,[x 7 140 14]); x=wd-300;
      set(pTop.hPlyLbl,ps,[x 3 40 14]); x=x+40;
      set(pTop.hPlySpd,ps,[x 3 40 14]); x=x+40;
      set(pTop.hTmLbl,ps,[x 3 30 14]); x=x+30;
      set(pTop.hTmVal,ps,[x 3 45 14]); x=x+50;
      set(pTop.hFrmLbl,ps,[x 3 35 14]); x=x+35;
      set(pTop.hFrmInd,ps,[x 3 50 16]); x=x+50;
      set(pTop.hFrmNum,ps,[x 3 60 14]);
      % position stuff in bottom panel
      x=5; wd0=wd-325;
      set(pBot.hStrm,ps,[x 5 35 20]); x=x+40;
      set(pBot.hBh,  ps,[x 5 85 20]); x=x+95;
      set(pBot.hArrL,ps,[x 4 25 20]); x=x+25;
      set(pBot.hSl, ps,[x 4 wd0 20]); x=x+wd0;
      set(pBot.hArrR,ps,[x 4 25 20]); x=x+30;
      set(pBot.hSetL,ps,[x 4 20 20]); x=x+20;
      set(pBot.hSetR,ps,[x 4 20 20]); x=x+25;
      set(pBot.hZm,  ps,[x 5 50 20]); x=x+60;
      set(pBot.hDel, ps,[x 4 20 20]);
      % update display
      if(~isempty(dispApi)); dispApi.requestUpdate(); end;
    end

    function exitProg( h, evnt ) %#ok<INUSD>
      menuApi.vidClose();
    end
  end

  function api = annMakeApi()
    % create api
    clrs=[]; fM=0;
    api = struct( 'updateDisp',@updateDisp, 'updateAnn',@updateAnn, ...
      'setFrame',@setFrame, 'buttonPress',@buttonPress, ...
      'insertBh',@insertBh,'setCenter',@setCenter,'dispTrack',@dispTrack );
    set( pBot.hStrm, 'callback', @(h,evnt) buttonPress('setStrm'));
    set( pBot.hBh,   'callback', @(h,evnt) buttonPress('setType'));
    set( pBot.hSetL, 'callback', @(h,evnt) buttonPress('moveLeft'));
    set( pBot.hSetR, 'callback', @(h,evnt) buttonPress('moveRight'));
    set( pBot.hDel,  'callback', @(h,evnt) buttonPress('delete'));
    set( pBot.hArrL, 'callback', @(h,evnt) buttonPress('prevBeh'));
    set( pBot.hArrR, 'callback', @(h,evnt) buttonPress('nextBeh'));
    set( pBot.hSl,   'callback', @(h,evnt) setFrame(0));
    set( pBot.hZm,   'callback', @(h,evnt) buttonPress('setZoom'));
    set( pTop.hBh,   'callback', @(h,evnt) buttonPress('gotoBeh'));

    function updateAnn()
      % new annotation loaded or annotation closed
      isAnn=~isempty(A); if(isAnn), en='on'; else en='off'; end;
      hAll=[pTop.hBh pBot.hBh pBot.hSetL pBot.hArrL pBot.hSl ...
        pBot.hArrR pBot.hSetR pBot.hZm pBot.hDel];
      set(hAll,'Enable',en); set([pBot.hBar pMid.hBar],'Visible',en);
      if(~isAnn), strs={''}; else strs=A.getNames(); end
      set(pBot.hBh,'Value',1,'String',strs);
      set(pTop.hBh,'Value',1,'String',{''});
      if(~isAnn), strs={''}; else strs=int2str2(1:A.nStrm()); end
      set(pBot.hStrm,'Value',1,'String',strs);
      if(isAnn), clrs=[.3 .3 .3; uniqueColors(ceil((A.k()-1)/6),6)]; end
      if(~isempty(pMid.hTxt)), set(pMid.hTxt,'String',''); end
      % update context menu for inserting behaviors
      if(isAnn), ns=A.getNames(); ks=A.getKeys(); k=A.k(); else k=0; end
      delete(pMid.hAdds); pMid.hAdds=zeros(1,k);
      ls=cell(1,k); for t=1:k, ls{t}=[ns{t} ' (' ks(t) ')']; end
      for t=1:k, pMid.hAdds(t)=uimenu(pMid.hAdd,'Label',ls{t}); end
      for t=1:k, set(pMid.hAdds(t),'callback',@(h,e) insertBh(t)); end
      % finally update display
      updateDisp();
    end

    function updateDisp()
      if(isempty(A)), return; end
      % get current frame / annotation info
      f=dispApi.getFrame(); id=A.getId(f); type=A.getType(id); n=A.n();
      nFrame=A.nFrame(); bs=A.getBnds(); ts=A.getTypes(); ns=A.getNames();
      % set bounds for zoomed slider
      h=pBot.hZm; s=get(h,'String'); v=get(h,'Value'); w=str2double(s{v});
      fL=max(0,fM-w/2); fR=min(nFrame,fL+w); fL=max(0,fR-w);
      % update standard GUI controls
      clr={'ForegroundColor',clrs(type,:)};
      set(pBot.hBh,'Value',type,clr{:}); ss=ns(ts);
      for i=1:n, ss{i}=sprintf('%i-%i %s',bs(i)+1,bs(i+1),ss{i}); end
      set(pTop.hBh,'String',ss,'Value',id,clr{:});
      if( nFrame==1 ), set(pBot.hSl,'Enable','off'); else
        set(pBot.hSl,'Min',fL,'Max',fR-1,'Value',f,'Enable','on');
        s=1/(fR-fL-1); set(pBot.hSl,'SliderStep',[s s]);
      end
      % display text label containing info about every stream
      strm0=get(pBot.hStrm,'Value'); nSt=A.nStrm(); s=cell(1,nSt);
      for i=1:nSt, A.setStrm(i); j=A.getId(f);
        s{i}=A.getName(j); if(i==strm0 && nSt>1), s{i}=['[' s{i} ']']; end
        s{i}=['\color[rgb]{' num2str(clrs(A.getType(j),:)) '}' s{i} ' / '];
      end
      s=[s{:}]; s=s(1:end-3); s(s=='_')='-'; A.setStrm(strm0);
      set(pMid.hTxt,'String',s);
      % update hBar for zoomed slider
      idL=A.getId(fL); idR=A.getId(fR-1);
      bs1=[fL bs(idL+1:idR) fR]-fL; ts1=ts(idL:idR);
      colBar( pBot.hBar, pBot.hSl, bs1, ts1, 4 );
      % update hBar for main slider
      colBar( pMid.hBar, pMid.hSl, bs, ts, 6, fL, fR );
      % finally backup annotation
      menuApi.annBackup();

      function colBar( hBar, hSl, bs, ts, h, fL, fR )
        % set position of hBar, adjust according to width of slider (8)
        p=get(hSl,'Position'); w0=16; w1=(p(3)-2*w0)/(bs(end)-bs(1));
        p(1)=p(1)+w0; p(2)=p(2)+p(4)-h; p(3)=p(3)-2*w0; p(4)=h;
        if(w1<8), p(1)=p(1)+4; p(3)=p(3)+round(w1)-8; end
        set(hBar,'Position',p); w=ceil(p(3)); h=ceil(p(4));
        % get hbar image
        nFrame=bs(end); bs=round(bs/nFrame*w); I=zeros(w,1);
        for i1=1:length(ts), I(bs(i1)+1:bs(i1+1),:)=ts(i1); end
        I=permute(clrs(I,:),[3 1 2]); I=I(ones(1,h),:,:);
        if(nargin==7)
          fL=round(fL/nFrame*w); fR=max(1,round(fR/nFrame*w));
          I(:,[fL+1 fR],:)=1; I([1 end],fL+1:fR,:)=1;
        end
        set(hBar,'CData',I);
      end
    end

    function setFrame(flag)
      assert(~isempty(A)); rng=get(pBot.hSl,{'Min','Max'});
      f = min(max(round(get(pBot.hSl,'Value')+flag),rng{1}),rng{2});
      setCenter = (f==rng{1} && f>0) || (f==rng{2} && f<A.nFrame()-1);
      set(pBot.hSl,'Value',f); dispApi.setFrame(f,0,setCenter);
    end

    function buttonPress( str )
      assert(~isempty(A));
      f=dispApi.getFrame(); id=A.getId(f);
      switch str
        case 'setStrm'
          strm = get(pBot.hStrm,'Value');
          A.setStrm(strm); dispApi.requestUpdate();
        case 'setType'
          type = get(pBot.hBh,'Value');
          A.setType(id,type); dispApi.requestUpdate();
        case 'moveRight'
          if(id==A.n()), return; end %id==A.getId(get(pBot.hSl,'Max'))
          A.move(id+1,f+1); dispApi.requestUpdate();
        case 'moveLeft'
          if(id==1), return; end %id==A.getId(get(pBot.hSl,'Min'))
          A.move(id,f); dispApi.requestUpdate();
        case 'delete'
          A.setType(id,1); %A.delete( id );
          dispApi.requestUpdate();
        case 'prevBeh'
          f1=A.getStart(id); if(id>1 && f==f1), f1=A.getStart(id-1); end
          dispApi.setFrame(f1,0);
        case 'nextBeh'
          f1=A.getEnd(id); if(id<A.n()), f1=f1+1; end
          dispApi.setFrame(f1,0);
        case 'setZoom'
          dispApi.setFrame(f,0)
        case 'gotoBeh'
          f1=A.getStart(get(pTop.hBh,'Value'));
          dispApi.setFrame(f1,0);
        otherwise
          assert(false);
      end
    end

    function insertBh( type )
      if(isempty(A)), return; end
      A.add(type,dispApi.getFrame());
      dispApi.requestUpdate();
    end

    function setCenter( frame ), fM=round(frame); end

    function hs = dispTrack( frame )
      if(isempty(trk)), hs=[]; return; end; hs0=[]; hs1=[]; hs2=[];
      cols=[0 1 0; 1 0 0; 0 0 1; 1 1 0; 0 1 1; 1 0 1];
      if(isfield(trk,'bbs')) % display detections
        bbs=trk.bbs{frame+1}; if(isempty(bbs)), bbs=zeros(0,6); end
        bbs=bbs(bbs(:,5)>0,:); n=size(bbs,1);
        hs0=zeros(1,n); rPrp={'LineWidth',2,'LineStyle','-','EdgeColor'};
        if(size(bbs,2)==6), is=mod(bbs(:,6)-1,6)+1; else is=ones(1,n); end
        for b=1:n, hs0(b)=rectangle('Position',bbs(b,1:4),rPrp{:},...
            cols(is(b),:)); end
      end
      if(isfield(trk,'Y')) % display (x,y) tracking results
        if(iscell(trk.Y)), Y=trk.Y{frame+1}; else Y=trk.Y(:,:,frame+1); end
        if(isempty(Y)), Y=zeros(0,2); end; m=size(Y,1); hs1=zeros(1,m);
        Y=max(Y,1); inf=dispApi.getInfo(); r=min(inf.width,inf.height)/70;
        Y(:,2)=min(Y(:,2),inf.height); Y(:,1)=min(Y(:,1),inf.width);
        rPrp={'LineWidth',2,'LineStyle','-','Curvature',[1 1],...
          'EdgeColor','w','FaceColor'};
        %D=pdist2(Y,Y)+diag(nan(m,1)); close=min(D)<125^2;
        %for i=1:m, if(close(i)), cols(i,:)=max(.7,cols(i,:)); end; end
        for i=1:m, hs1(i)=rectangle(rPrp{:},cols(i,:),...
            'Position',[Y(i,:)-r 2*r 2*r]); end
      end
      if(isfield(trk,'E')) % display ellipse tracking results
        if(iscell(trk.E)), E=trk.E{frame+1}; else E=trk.E(:,:,frame+1); end
        if(isempty(E)), E=zeros(0,5); end; m=size(E,1); hs2=zeros(3,m);
        hold on; for i=1:m, [hs2(1,i),hs2(2,i),hs2(3,i)]=plotEllipse(...
            E(i,1),E(i,2),E(i,3),E(i,4),E(i,5),cols(i,:),[],2); end
        hold off; hs2=hs2(:)';
      end
      hs = [hs0 hs1 hs2];
    end
  end

  function api = dispMakeApi()
    % create api
    [sr, audio, info, nFrame, speed, curInd, hs, ...
      hImg, needUpdate, prevTime, looping ]=deal([]);
    api = struct( 'setVid',@setVid, 'setAud',@setAud, ...
      'setSpeedCb',@setSpeedCb, 'getInfo',@getInfo, ...
      'getFrame',@getFrame, 'setFrame',@setFrame, ...
      'requestUpdate',@requestUpdate, 'exportVid',@exportVid );
    set(pMid.hSl,    'Callback',@(h,evnt) setFrameCb(0));
    set(pTop.hFrmInd,'Callback',@(h,evnt) setFrameCb(1));

    function setVid( sr1 )
      % reset local variables
      if(isstruct(sr)), sr=sr.close(); end
      if(~isempty(hs)), delete(hs); hs=[]; end
      [sr, audio, info, nFrame, speed, curInd, hs, ...
        hImg, needUpdate, prevTime, looping ]=deal([]);
      sr=sr1; nFrame=0; looping=0; speed=-1; setFrame( 0, 0 );
      % update GUI
      if(~isstruct(sr)), cla(pMid.hAx); pMid.hTxt=[]; else
        info=sr.getinfo(); nFrame=info.numFrames;
        sr.seek(0); s=max(1e-6,1/(nFrame-1)); ss={'SliderStep',[s,s]};
        if(nFrame>1), set(pMid.hSl,'Max',nFrame-1,ss{:}); end
        hImg = imshow( zeros(info.height,info.width,'uint8') );
        set(hImg,'UIContextMenu',pMid.hAdd);
        pMid.hTxt=text(0,0,'','FontSize',25,'Units','pixels',...
          'HorizontalAlignment','center', 'VerticalAlignment','top');
        %fprintf('fps of video = %f\n', info.fps); % temp display
      end
      set(pMid.hSl,'Value',0); v=(nFrame>1)+1;
      en={'off','on'}; set(pMid.hSl,'Enable',en{v});
      en={'inactive','on'}; set(pTop.hFrmInd,'Enable',en{v});
      set(pTop.hFrmInd,'String','0'); set(pTop.hTmVal,'String','0:00');
      set(pTop.hFrmNum,'String',[' / ' int2str(nFrame)]);
      % update display
      feval(get(hFig,'ResizeFcn'));
      requestUpdate();
    end

    function exportVid( nm )
      % prompt for range of frames to export
      prompt={'Start frame:','End frame:','Frame Skip','Quality (0-100)'};
      wi=getInfo(); dfs={'1', num2str(wi.numFrames),'1','80'};
      rng=str2double(inputdlg(prompt,'Select Export Range',1,dfs));
      if(any(isnan(rng)) || rng(1)>rng(2)), error('invalid range'); end
      f0=max(1,rng(1)); f1=min(rng(2),wi.numFrames); skip=max(1,rng(3));
      wi.codec='jpg'; wi.quality=rng(4);
      try %#ok<ALIGN> % export video
        for f=f0:skip:f1
          setFrame(f,0); I=getframe(pMid.hAx); I=I.cdata;
          c=2; I=I(1+c:end-c,1+c:end-c,:); h=size(I,1); w=size(I,2);
          if(f==f0), wi.height=h; wi.width=w; sw=seqIo(nm,'w',wi); end
          sw.addframe(I);
        end; sw.close();
      catch err, sw.close(); throw(err); end
    end

    function setAud(y,fs,nb)
      a = abs(1 - (length(y)/fs) / (nFrame/info.fps));
      if(a>.01), error('Audio/video mismatch.'); end
      audio.fPlay=audioplayer(y,fs,nb);
      audio.bPlay=audioplayer(flipud(y),fs,nb);
      audio.fs=fs; audio.ln=length(y); setSpeed(speed);
      requestUpdate();
    end

    function dispLoop()
      if(looping), return; end; looping=1;
      while( 1 )
        % exit if appropriate, or if vid not loaded do nothing
        if(~isstruct(sr) || ~isstruct(info)), looping=0; return; end

        % stop playing video if at begin/end
        if((speed>0&&curInd==nFrame-1) || (speed<0&&curInd==0))
          setSpeed(0); needUpdate=1;
        end

        % increment/decrement curInd appropriately
        if( speed~=0 )
          t=clock(); eTime=etime(t,prevTime);
          del = speed * max(10,info.fps) * min(.1,eTime);
          if( speed>0 ), del=min(del, nFrame-curInd-1 ); end
          if( speed<0 ), del=max(del, -curInd ); end
          setFrame(curInd+del, speed); prevTime=t; needUpdate=1;
        end

        % update display if necessary
        if(~needUpdate), looping=0; return; else
          sr.seek( round(curInd) ); I=sr.getframe();
          if(~isempty(hs)), delete(hs); hs=[]; end
          hs=annApi.dispTrack(round(curInd));
          assert(~isempty(I)); set(hImg,'CData',I);
          set(pMid.hSl,'Value',curInd); c=round(curInd/info.fps);
          set(pTop.hFrmInd,'String',int2str(round(curInd+1)));
          c0=floor(c/3600); c=mod(c,3600); c1=floor(c/60); c2=mod(c,60);
          if(~c0), set(pTop.hTmVal,'String',sprintf('%i:%02i',c1,c2)); else
            set(pTop.hTmVal,'String',sprintf('%i:%02i:%02i',c0,c1,c2)); end
          annApi.updateDisp(); needUpdate=false; drawnow();
        end
      end
    end

    function setSpeedCb( flag )
      if(~isstruct(sr)),return; end
      if( flag==0 )
        setSpeed( 0 );
      elseif( speed==0 )
        setSpeed( flag ); prevTime=clock();
      elseif( sign(speed)==flag && abs(speed)<2^12 )
        setSpeed( speed*2 );
      elseif( sign(speed)~=flag && abs(speed)>1/8 )
        setSpeed( speed/2 );
      elseif( sign(speed)~=flag )
        setSpeed( 0 );
      end
      requestUpdate();
    end

    function setSpeed( speed1 )
      if((speed1>0&&curInd==nFrame-1)||(speed1<0&&curInd==0)),speed1=0;end
      speed=speed1; p=abs(speed); if(speed<0), ss='-'; else ss=''; end
      if(p<1 && p~=0), s=['1/' int2str(1/p)]; else s=int2str(p); end
      set(pTop.hPlySpd,'String',[ss s 'x']);
      if(~isempty(audio))
        stop(audio.fPlay); stop(audio.bPlay); if(speed==0), return; end;
        st=curInd/(nFrame-1); if(speed<0), st=1-st; end; st=st*audio.ln+1;
        if(speed>0), plr=audio.fPlay; else plr=audio.bPlay; end;
        set(plr,'SampleRate',audio.fs*p); play(plr,st);
      end
    end

    function setFrameCb(flag)
      if( flag==0 )
        f=round(get(pMid.hSl,'Value')); set(pMid.hSl,'Value',f);
      elseif(flag==1)
        f=str2double(get(pTop.hFrmInd,'String'));
        if(isnan(f)), requestUpdate(); return; else f=f-1; end
      end
      setFrame(f,0);
    end

    function setFrame( curInd1, speed1, setCenter )
      curInd=max(0,min(curInd1,nFrame-1));
      if( speed~=speed1 ), setSpeed(speed1); end
      if(nargin<3 || setCenter), annApi.setCenter(curInd); end
      requestUpdate();
    end

    function requestUpdate(), needUpdate=true; dispLoop(); end

    function info1 = getInfo(), info1=info; end

    function curInd1 = getFrame(), curInd1=round(curInd); end
  end

  function api = menuMakeApi()
    % create api
    [fVid, fAnn, lastSave]=deal([]);
    api = struct('vidClose',@vidClose, 'annClose',@annClose, ...
      'vidOpen',@vidOpen, 'audOpen',@audOpen, 'trkOpen',@trkOpen, ...
      'annOpen',@annOpen, 'annBackup',@annBackup );
    set(menu.hVidOpn,'Callback',@(h,envt) vidOpen(1) );
    set(menu.hVidsOpn,'Callback',@(h,envt) vidOpen(2) );
    set(menu.hVidExp,'Callback',@(h,envt) vidExport() );
    set(menu.hVidCls,'Callback',@(h,envt) vidClose() );
    set(menu.hVidInf,'Callback',@(h,envt) vidInfo() );
    set(menu.hVidAud,'Callback',@(h,envt) audOpen() );
    set(menu.hAnnNew,'Callback',@(h,envt) annNew(0) );
    set(menu.hAnnOpn,'Callback',@(h,envt) annOpen(0) );
    set(menu.hAnnCls,'Callback',@(h,envt) annClose() );
    set(menu.hAnnSav,'Callback',@(h,envt) annSave() );
    set(menu.hAnnCnf,'Callback',@(h,envt) annNew(1) );
    set(menu.hAnnMrg,'Callback',@(h,envt) annOpen(1) );
    set(menu.hTrkOpn,'Callback',@(h,envt) trkOpen(0) );
    set(menu.hTrkCls,'Callback',@(h,envt) trkOpen([]) );

    function updateMenus()
      m=menu; if(isempty(fVid)), en='off'; else en='on'; end
      set([m.hVidExp m.hVidCls m.hVidInf m.hVidAud m.hAnnNew m.hAnnOpn ...
        m.hTrkOpn m.hTrkCls],'Enable',en);
      if(isempty(A)), en='off'; else en='on'; end
      set([m.hAnnSav m.hAnnCls m.hAnnCnf m.hAnnMrg],'Enable',en);
      nm='Caltech Behavior Annotator';
      if(~isempty(fVid)), [~,nm1]=fileparts(fVid); nm=[nm ' - ' nm1]; end
      set(hFig,'Name',nm); annApi.updateAnn(); dispApi.requestUpdate();
    end

    function vidClose()
      if(~isempty(A)), annClose(); end; trkOpen([]);
      fVid=[]; dispApi.setVid([]); updateMenus();
    end

    function vidOpen( flag )
      if(isempty(fVid)), d='.'; else d=fileparts(fVid); end
      if(all(ischar(flag)))
        [d,f]=fileparts(flag); if(isempty(d)), d='.'; end;
        d=[d '/']; f=[f '.seq']; flag=1;
      elseif( flag==1 )
        [f,d]=uigetfile('*.seq','Select video',[d '/*.seq']);
      elseif( flag==2 )
        [f,d]=uigetfile('*.seq','Select first video',[d '/*.seq']);
        [f2,d2]=uigetfile('*.seq','Select second video',[d '/*.seq']);
        if( f2==0 ), return; end
      end
      if( f==0 ), return; end; vidClose(); fVid=[d f];
      try
        if( flag==1 ), sr=seqIo(fVid,'r'); else
          sr=seqIo({fVid,[d2 f2]},'rdual'); end
        dispApi.setVid(sr); updateMenus();
      catch er
        errordlg(['Failed to load: ' fVid '. ' er.message],'Error');
        fVid=[]; return;
      end
    end

    function vidInfo()
      info=dispApi.getInfo(); txt=evalc('disp(info);'); %#ok<NASGU>
      while(txt(end)==10), txt=txt(1:end-1); end; txt=[char([10 10]) txt];
      pos=get(hFig,'Position'); pos(1:2)=pos(1:2)+100; pos(3:4)=[400 400];
      h=figure('NumberTitle','off', 'Toolbar','auto', ...
        'Color','k', 'MenuBar','none', 'Visible','on', ...
        'Name','', 'Resize','off','Position',pos); pos(1:2)=0;
      uicontrol(h,'Style','text','FontSize',10,'BackgroundColor','w',...
        'HorizontalAlignment','left','Position',pos,'String',txt,...
        'FontName','FixedWidth');
    end

    function vidExport()
      assert(~isempty(fVid)); fNm=[fVid(1:end-4) '-copy.seq'];
      [f,d] = uiputfile('*.seq','Select export file',fNm);
      if(f==0), return; end; f=[d f];
      try dispApi.exportVid(f); catch er
        errordlg(['Failed to export: ' f '. ' er.message],'Error'); end
    end

    function audOpen( fAud )
      if( nargin==0 ), [f,d]=uigetfile('*.wav','Select audio',...
          [fVid(1:end-3) 'wav']); if(f==0), return; end; fAud=[d f]; end
      try
        [y,fs,nb]=wavread(fAud); dispApi.setAud(y,fs,nb); %#ok<REMFF1>
      catch er
        errordlg(['Failed to load: ' fAud '. ' er.message],'Error');
      end
    end

    function annClose()
      assert(~isempty(A)); qstr='Save Current Annotation?';
      button = questdlg(qstr,'Save','yes','no','yes');
      if(strcmp(button,'yes')); annSave(); end
      [fAnn,A,lastSave]=deal([]); updateMenus();
    end

    function annOpen( flag )
      assert(~isempty(fVid)); A1=A; if(~isempty(A)), annClose(); end; e='';
      if( all(ischar(flag)) ), f=flag; flag=0;
        [d,f,e]=fileparts(f); if(isempty(d)), d='.'; end; d=[d '/'];
        if(isempty(e) && exist([d f '.txt'],'file')), e='.txt'; end
        if(isempty(e) && exist([d f '.bAnn'],'file')), e='.bAnn'; end
      else
        if(isempty(fAnn)), fAnn=[fVid(1:end-3) 'txt']; end
        [f,d]=uigetfile('*.txt;*.bAnn','Select Annotation',fAnn);
      end
      if( f==0 ), return; end; fAnn=[d f e];
      try
        if( flag==1 )
          A=A1; assert(~isempty(A)); A.merge(fAnn);
        elseif( flag==0 )
          A=behaviorData('load',fAnn);
          info=dispApi.getInfo(); nFrame=info.numFrames;
          if(A.nFrame()~=nFrame), error('Annotation/video mismatch.'); end
        end
        updateMenus();
      catch er
        errordlg(['Failed to load: ' fAnn '. ' er.message],'Error');
        [fAnn,A,lastSave]=deal([]); return;
      end
    end

    function annSave()
      assert(~isempty(fVid) && ~isempty(A) && ~isempty(fAnn));
      [f,d] = uiputfile('*.txt;*.bAnn','Select annotation',fAnn);
      if( f==0 ), return; end; fAnn=[d f]; A.save(fAnn);
    end

    function annNew( update )
      assert(~isempty(fVid)); A1=A; if(~isempty(A)), annClose(); end
      fNm=[fileparts(fVid) '/config.txt'];
      [f,d] = uigetfile('*.txt','Select config file',fNm);
      if( f==0 ), return; end; f=[d f]; fAnn=[fVid(1:end-3) 'txt'];
      try
        if( update )
          A=A1; assert(~isempty(A)); A.recreate(f); updateMenus();
        else
          info=dispApi.getInfo(); nFrame=info.numFrames;
          A=behaviorData('create',f,nFrame); updateMenus();
        end
      catch er
        errordlg(er.message, 'File Error'); A=[];
      end
    end

    function annBackup()
      if( isempty(lastSave) || etime(clock,lastSave)>60 )
        [d,f]=fileparts(fAnn); f=[d '/' f '-backup.txt'];
        assert(~isempty(A)); A.save(f); lastSave=clock();
      end
    end

    function trkOpen( flag )
      if(isempty(flag)), trk=[]; dispApi.requestUpdate(); return; end
      if( all(ischar(flag)) ), f=flag; else fNm=[fVid(1:end-3) 'mat'];
        [f,d]=uigetfile('*.mat','Load Tracking',fNm); f=[d f];
      end
      assert(~isempty(fVid)); if(f==0), return; end
      try
        tmp=load(f); assert(any(isfield(tmp,{'Y','bbs','E'})));
        if(isfield(tmp,'Y')), trk.Y=tmp.Y; T=size(trk.Y,3);
        elseif(isfield(tmp,'bbs')), trk.bbs=tmp.bbs; T=length(trk.bbs);
        elseif(isfield(tmp,'E')), trk.E=tmp.E; T=size(trk.E,3);
        end
        info=dispApi.getInfo(); nFrame=info.numFrames;
        if(T~=nFrame), error('Tracking/video mismatch.'); end
      catch er
        trk=[]; errordlg(['Failed to load: ' f '. ' er.message],'Error');
      end
      dispApi.requestUpdate();
    end

  end
end

author: zouchuhang

repo: LayoutNet


text:
% ------------------------------------------------------------------------
%        Get ROI reference-space locations and region annotations
% ------------------------------------------------------------------------

%% SET FILE LOCATIONS OF TRANSFORM AND ROIS

% file location of transform and transformed image
transform_location =         'C:\Drive\Histology\for tutorial - sample data\SS096_done\processed\transformations\SS096_1_1_transform_data';
transformed_slice_location = 'C:\Drive\Histology\for tutorial - sample data\SS096_done\processed\transformations\SS096_1_1_transformed.tif';

% file location of ROIs (image / array of the same size as the reference ie 800 x 1140)
roi_location = 'C:\Drive\Histology\for tutorial - sample data\SS096_1_1_ROIs.tif';

% directory of reference atlas files
annotation_volume_location = 'C:\Drive\Histology\for tutorial\annotation_volume_10um_by_index.npy';
structure_tree_location = 'C:\Drive\Histology\for tutorial\structure_tree_safe_2017.csv';

% plane used to view ROI ('coronal' -- most common, 'sagittal', 'transverse')
plane = 'coronal';

% Synthetic ROIs for testing
% rois = zeros(800,1140,'uint8');
% rois(250:300, 600:610) = 200; rois(480:500, 200:210) = 200;
% imwrite(rois,roi_location)
%
% Using a set of x and y coordinates from the ImageJ function Analyze Particles to generate an ROI image
% roi_array = zeros(800,1140,'uint8');
% roi_array_values = csvread('C:\ROI_files\cfos_cells.csv', 1, 5);
% y = roi_array_values(:, 1);
% x = roi_array_values(:, 2);
%
% for i = 1:length(roi_array_values)-1
%     roi_array(x(i),y(i)) = 200;
% end

%% LOAD THE DATA

% load the transformed slice image
transformed_slice_image = imread(transformed_slice_location);

% load the transform from the transform file
transform_data = load(transform_location);
transform_data = transform_data.save_transform;

% get the actual transformation from slice to atlas
slice_to_atlas_transform = transform_data.transform;

% get the position within the atlas data of the transformed slice
slice_num = transform_data.allen_location{1};
slice_angle = transform_data.allen_location{2};

% load the rois
rois = imread(roi_location);

% if the rois come from a transformed roi image of non-contiguous roi
% pixels (e.g. an ROI pixel for each neuron), then run this line to ensure
% a one-to-one mapping between ROIs in the original and transformed images:
% rois = uint8(imregionalmax(rois));

% load the reference brain annotations
if ~exist('av','var') || ~exist('st','var')
    disp('loading reference atlas...')
    av = readNPY(annotation_volume_location);
    st = loadStructureTree(structure_tree_location);
end

% select the plane for the viewer
if strcmp(plane,'coronal')
    av_plot = av;
elseif strcmp(plane,'sagittal')
    av_plot = permute(av,[3 2 1]);
elseif strcmp(plane,'transverse')
    av_plot = permute(av,[2 3 1]);
end

%% GET REFERENCE-SPACE LOCATIONS AND REGION ANNOTATIONS FOR EACH ROI

% I will do this for every *pixel* in the roi image with nonzero value,
% but this code can be modified, e.g. to do it by clusters of pixels

% show the transformed ROI, together with the transformed image
figure; imshow(imfuse(rois, transformed_slice_image));
title('transformed slice image, fused with ROIs')

% make sure the rois are in a properly size image
assert(size(rois,1)==800&size(rois,2)==1140&size(rois,3)==1,'roi image is not the right size');

% initialize array of locations (AP, DV, ML relative to bregma) in reference space
% and the correponding region annotations
roi_location = zeros(sum(rois(:)>0),3);
roi_annotation = cell(sum(rois(:)>0),3);

% get location and annotation for every roi pixel
[pixels_row, pixels_column] = find(rois>0);

% generate other necessary values
bregma = allenCCFbregma(); % bregma position in reference data space

atlas_resolution = 0.010; % mm
ref_size = size(squeeze(av_plot(1,:,:)));
offset_map = get_offset_map(slice_angle, ref_size);

% loop through every pixel to get ROI locations and region annotations
for pixel = 1:length(pixels_row)

    % get the offset from the AP value at the centre of the slice, due to
    % off-from-coronal angling
    offset = offset_map(pixels_row(pixel),pixels_column(pixel));

    % use this and the slice number to get the AP, DV, and ML coordinates
    if strcmp(plane,'coronal')
        ap = -(slice_num-bregma(1)+offset)*atlas_resolution;
        dv = (pixels_row(pixel)-bregma(2))*atlas_resolution;
        ml = (pixels_column(pixel)-bregma(3))*atlas_resolution;
    elseif strcmp(plane,'sagittal')
        ml = -(slice_num-bregma(3)+offset)*atlas_resolution;
        dv = (pixels_row(pixel)-bregma(2))*atlas_resolution;
        ap = -(pixels_column(pixel)-bregma(1))*atlas_resolution;
    elseif strcmp(plane,'transverse')
        dv = -(slice_num-bregma(2)+offset)*atlas_resolution;
        ml = (pixels_row(pixel)-bregma(3))*atlas_resolution;
        ap = -(pixels_column(pixel)-bregma(1))*atlas_resolution;
    end

    roi_location(pixel,:) = [ap dv ml];

    % finally, find the annotation, name, and acronym of the current ROI pixel
    ann = av_plot(slice_num+offset,pixels_row(pixel),pixels_column(pixel));
    name = st.safe_name{ann};
    acr = st.acronym{ann};

    roi_annotation{pixel,1} = ann;
    roi_annotation{pixel,2} = name;
    roi_annotation{pixel,3} = acr;

end

 roi_table = table(roi_annotation(:,2),roi_annotation(:,3), ...
                    roi_location(:,1),roi_location(:,2),roi_location(:,3), roi_annotation(:,1), ...
     'VariableNames', {'name', 'acronym', 'AP_location', 'DV_location', 'ML_location', 'avIndex'});

 disp(roi_table(1:10,:))

% now, use roi_locations and roi_annotations for your further analyses

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function pass = test_periodic(pref)
% Test 'periodic' syntax for linear ODEs.

if ( nargin == 0 )
    pref = cheboppref();
end
tol = 1e4*pref.bvpTol;

%% A simple ODE:

% Chebop:
N = chebop(@(u) diff(u, 2) - u);

% Periodic RHS:
rhs = chebfun(@(x) sin(pi*x));

% Apply periodic BCs manually:
N.bc = @(u) [u(-1) - u(1) ; feval(diff(u),-1) - feval(diff(u), 1)];
u = mldivide(N, rhs, pref);

% Apply via 'periodic':
N.bc = 'periodic';
v = mldivide(N, rhs, pref);

% Compare results:
pass(1) = norm(u - v) < tol;

%% A periodic piecewise system:

d = [-pi, 0, pi];
A = chebop(d);
A.op = @(x, u, v) [u-diff(v) ; diff(u,2)+v];
x = chebfun('x',d);
f = [chebfun(0, d) ; cos(x)];
A.bc = 'periodic';
uv = mldivide(A, f, pref);

trueSoln = [cos(x+3*pi/4)/sqrt(2) ; cos(x+pi/4)/sqrt(2)];
pass(2) = norm(uv - trueSoln) < tol;

%% Eigenvalue problem:

d = [-pi, 0, pi];
A = chebop(d);
A.op = @(x, u, v) [u-diff(v) ; diff(u,2)+v];
A.bc = 'periodic';
B = chebop(d);
B.op = @(x, u, v) [v + u ; diff(v)];

[V, D] = eigs(A, B, 5, 0, pref);
e = diag(D);

% Sort the eigenvalues to ensure things will work on all machines. Eigs() does
% not appear to return the eigenvalues in a consistent way for different
% machines. We expect five pair of eigenvalues to appear, check whether they are
% all there
% `sort(x)` sorts complex values by abs() and then by angle(). In order have
% consistent sorting, we move all the eigenvalues up into the first quadrant
% before sorting them.
[ignored, idx] = sort(real(e));
e = e(idx);

e12 = e(1:2);
e35 = e(3:5);
[ignored, idx] = sort(imag(e12));
e12 = e12(idx);
[ignored, idx] = sort(imag(e35));
e35 = e35(idx);

e = [e12; e35];

pass(3) = norm(real(e) - [0 0 1 1 1].', inf) + ...
    norm(imag(e) - [-1 1 -1 0 1].', inf) < tol;
pass(4) = norm(V{1}(pi) - V{1}(pi), inf) + norm(V{2}(pi) - V{2}(pi), inf) < tol;

%% Test the TRIGCOLLOC class. FIRST ORDER AND CONSTANT COEFFICIENTS:
%  u' + u = cos(x), on [0 2*pi].

% Set domain, operator L, and rhs f.
dom = [0 2*pi];
L = chebop(@(u) diff(u) + u, dom);
f = chebfun(@(x) cos(x), dom);

% Solve with TRIGCOLLOC.
L.bc = 'periodic';
u = L \ f;

% Compare with exact solution.
exact = chebfun(@(x) 1/2*cos(x) + 1/2*sin(x), dom, 'periodic');
pass(5) = norm(u - exact, inf) < tol;
pass(6) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test the TRIGCOLLOC class. FIRST ORDER AND VARIABLES COEFFICIENTS:
%  u' + (1+cos(x))u = cos(2x), on [-2*pi 2*pi].

% Set domain, c, and rhs f.
dom = [-2*pi 2*pi];
f = chebfun(@(x) cos(2*x), dom);

% Set chebop L. We construct the variable coefficient inside the chebop.
L = chebop(@(x, u) diff(u) + (1 + cos(x)).*u, dom);
L.bc = 'periodic';

% Solve with TRIGCOLLOC.
u = L \ f;

pass(7) = norm(L*u - f) < tol;
pass(8) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(9) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test the TRIGCOLLOC class. SECOND ORDER AND CONSTANT COEFFICIENTS:
%  u'' + 10u' + 5u = cos(x), on [-2*pi 2*pi].

% Set domain, constants coefficients a and b, operator L,
% and rhs f.
dom = [-2*pi 2*pi];
a = 10;
b = 5;
L = chebop(@(u) diff(u, 2) + a*diff(u) + b*u, dom);
L.bc = 'periodic';
f = chebfun(@(x) cos(x), dom);

% Solve with TRIGCOLLOC.
u = L \ f;

% Compare with exact solution.
exact = chebfun(@(x) 1/29*cos(x) + 5/58*sin(x), dom, 'periodic');
pass(10) = norm(u - exact, inf) < tol;
pass(11) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test the TRIGCOLLOC class. SECOND ORDER AND VARIABLE COEFFICIENTS:
%  (2+cos(4x))u'' + sin(cos(2x))u' + exp(cos(x))u = cos(x), on [-pi pi].

% Set domain, variable coefficients a, b and c, and rhs f.
dom = [-pi pi];
a = chebfun(@(x) 2 + cos(4*x), dom);
b = chebfun(@(x) sin(cos(2*x)), dom, 'periodic');
c = chebfun(@(x) exp(cos(x)), dom);
f = chebfun(@(x) cos(x), dom);

% Set chebop. The variale coefficients have been constructed outside the
% chebop, some with 'periodic', some without it.
L = chebop(@(u) a.*diff(u, 2) + b.*diff(u) + c.*u, dom);
L.bc = 'periodic';

% Solve with TRIGCOLLOC.
u = L \ f;

pass(12) = norm(L*u - f) < tol;
pass(13) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(14) = abs(feval(diff(u), dom(1)) - feval(diff(u), dom(2))) < tol;
pass(15) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test the TRIGCOLLOC class. THIRD ORDER AND VARIABLE COEFFICIENTS:
%  (2+cos(x))u''' + sin(cos(2x))u'' + exp(cos(x))u' + sin(x)u = cos(x),
%  on [-pi pi].

% Set domain, variable coefficients a, b, c and d, and rhs f.
dom = [-pi pi];
a = chebfun(@(x) 2 + cos(x), dom);
b = chebfun(@(x) sin(cos(2*x)), dom);
c = chebfun(@(x) exp(cos(x)), dom);
d = chebfun(@(x) sin(x), dom);
f = chebfun(@(x) cos(x), dom);

% Set chebop.
L = chebop(@(u) a.*diff(u, 3) + b.*diff(u, 2) + c.*diff(u) + d.*u, dom);
L.bc = 'periodic';

% Solve with TRIGCOLLOC.
u = L \ f;

pass(16) = norm(L*u - f) < tol;
pass(17) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(18) = abs(feval(diff(u), dom(1)) - feval(diff(u), dom(2))) < tol;
pass(19) = abs(feval(diff(u, 2), dom(1)) - feval(diff(u, 2), dom(2))) < tol;
pass(20) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test the TRIGCOLLOC class. FOURTH ORDER AND VARIABLE COEFFICIENTS:
%  (2+cos(x))u'''' + sin(cos(2x))u''' + exp(cos(x))u'' + ...
%  sin(x)u' + sin(2*x)u = cos(10*x), on [-pi pi].

% Set domain, variable coefficients aa, bb, cc, dd and ee, and rhs f.
dom = [-pi pi];
a = chebfun(@(x) 2 + cos(x), dom);
b = chebfun(@(x) sin(cos(2*x)), dom);
c = chebfun(@(x) exp(cos(x)), dom);
d = chebfun(@(x) sin(x), dom);
e = chebfun(@(x) sin(2*x), dom);
f = chebfun(@(x) cos(10*x), dom);

% Set chebop.
L = chebop(@(u) a.*diff(u, 4) + b.*diff(u, 3) + c.*diff(u, 2) + ...
    d.*diff(u) + e.*u, dom);
L.bc = 'periodic';

% Solve with TRIGCOLLOC.
u = L \ f;

pass(21) = norm(L*u - f) < tol;
pass(22) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(23) = abs(feval(diff(u), dom(1)) - feval(diff(u), dom(2))) < tol;
pass(24) = abs(feval(diff(u, 2), dom(1)) - feval(diff(u, 2), dom(2))) < tol;
pass(25) = abs(feval(diff(u, 3), dom(1)) - feval(diff(u, 3), dom(2))) < tol;
pass(26) = isequal(get(u.funs{1}, 'tech'), @trigtech);

%% Test breakpoint introduced by the domain.
%  u' + u = cos(x), on [0 pi 2*pi].

dom = [0 pi 2*pi];
L = chebop(@(u) diff(u) + u, dom);
f = chebfun(@(x) cos(x), dom);
L.bc = 'periodic';
u = L \ f;

pass(27) = norm(L*u - f) < tol;
pass(28) = abs(u(dom(1)) - u(dom(end))) < tol;
pass(29) = isequal(get(u.funs{1}, 'tech'), @chebtech2);

%% Test breakpoint introduced by a coefficient.
%  u'' + abs(x)u = 1, on [-1 1].

dom = [-1 1];
L = chebop(@(x,u) diff(u,2) + abs(x).*u, dom);
L.bc = 'periodic';
u = L \ 1;

pass(30) = norm(L*u - 1) < tol;
pass(31) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(32) = isequal(get(u.funs{1}, 'tech'), @chebtech2);

%% Test the TRIGSPEC class. FIRST ORDER AND VARIABLES COEFFICIENTS:
%  u' + (1+cos(x))u = cos(2x), on [-2*pi 2*pi].

% Set domain, c, and rhs f.
dom = [0 2*pi];
f = chebfun(@(x) cos(2*x), dom);

% Set chebop L. We construct the variable coefficient inside the chebop.
L = chebop(@(x, u) diff(u) + (1 + cos(x)).*u, dom);
L.bc = 'periodic';

% Solve with TRIGSPEC.
pref.discretization = @trigspec;
u = solvebvp(L, f, pref);

pass(33) = norm(L*u - f) < tol;
pass(34) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(35) = isequal(get(u.funs{1}, 'tech'), @trigtech);
pass(36) = isreal(u);

%% Test the TRIGSPEC class. SECOND ORDER AND VARIABLE COEFFICIENTS:
%  (2+cos(4x))u'' + sin(cos(2x))u' + exp(cos(x))u = cos(x), on [-pi pi].

% Set domain, variable coefficients a, b and c, and rhs f.
dom = [-pi pi];
a = chebfun(@(x) 2 + cos(4*x), dom);
b = chebfun(@(x) sin(cos(2*x)), dom, 'periodic');
c = chebfun(@(x) exp(cos(x)), dom);
f = chebfun(@(x) cos(x), dom);

% Set chebop. The variale coefficients have been constructed outside the
% chebop, some with 'periodic', some without it.
L = chebop(@(u) a.*diff(u, 2) + b.*diff(u) + c.*u, dom);
L.bc = 'periodic';

% Solve with TRIGSPEC.
pref.discretization = 'coeffs';
u = solvebvp(L, f, pref);

pass(37) = norm(L*u - f) < tol;
pass(38) = abs(u(dom(1)) - u(dom(2))) < tol;
pass(39) = abs(feval(diff(u), dom(1)) - feval(diff(u), dom(2))) < tol;
pass(40) = isequal(get(u.funs{1}, 'tech'), @trigtech);
pass(41) = isreal(u);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function train_CNN_model(varargin)

setup('useGPU', true);

tic

% Load the data
root_loc = 'D:\Datasets\detection_validation';
location = [root_loc, '/prep_data/'];

faceCheckersLoc_train = dir([location 'face_validator_train_*']);
faceCheckersLoc_test = dir([location 'face_validator_test_*']);

corrs_all = [];
rmses_all = [];

faceCheckers = struct;

% As we will be training a classifier that will act as a regressor,
% binarize it
num_bins = 30;

for i=1:numel(faceCheckersLoc_train)

    load([location faceCheckersLoc_train(i).name]);

    % set a max value to the error
    errors(errors > 3) = 3;
    errors_train = quantizeContinuous(errors, 0, 3, num_bins);
    examples_train = single(examples);
    clear examples

    mean_ex = mean(examples_train);
    std_ex = std(examples_train);
    std_ex = std_ex / 256;
    examples_train = bsxfun(@times, bsxfun(@minus, examples_train, mean_ex), 1./std_ex);
    num_examples_train = size(examples_train,1);

    % Add rows and columns untill we have 60 x 60 image

    % keep adding rows
    while(size(mask,1) < 60)
        mask = cat(1, mask, false(1, size(mask,2)));
        triX = cat(1, triX, -ones(1, size(mask,2)));
    end

    % keep adding cols
    while(size(mask,2) < 60)
        mask = cat(2, mask, false(size(mask,1),1));
        triX = cat(2, triX, -ones(size(triX,1),1));
    end

    examples_r = single(zeros(size(mask, 1), size(mask, 2), num_examples_train));

    img_curr = zeros(size(mask));
    for e=1:num_examples_train

        img_curr(mask) = examples_train(e,:);
        examples_r(:, :, e) = img_curr;

    end

    imdb.images.data = examples_r;
    clear examples_r
    imdb.images.label = errors_train';
    imdb.images.id = 1:numel(errors_train);

    % Split data for training and validation (20%)
    imdb.images.set = ones(1, numel(errors_train));
    imdb.images.set(round(4*end/5):end) = 2;

    % Visualize some of the data
    figure(10) ; clf ; colormap gray ;
    subplot(1,2,1) ;
    pos_samples = squeeze(imdb.images.data(:,:,imdb.images.label==1));
    vl_imarraysc(pos_samples(:,:,1:20)) ;
    axis image off ;
    title('Positive training data') ;

    subplot(1,2,2) ;
    neg_samples = squeeze(imdb.images.data(:,:,imdb.images.label>5));
    vl_imarraysc(neg_samples(:,:,1:20)) ;
    axis image off ;
    title('Negative training data') ;

    net = initializeFaceCNN(num_bins) ;

    trainOpts.batchSize = 100 ;
    trainOpts.numEpochs = 20;
    trainOpts.continue = true ;
    trainOpts.gpus = [1] ;
    trainOpts.learningRate = 0.001 ;
    trainOpts.expDir = ['trained/intermediate/face_validator_' num2str(i)];
    trainOpts.errorFunction = 'regression';
    trainOpts = vl_argparse(trainOpts, varargin);

    % Call training function in MatConvNet
    [net,info] = cnn_train_reg(net, imdb, @getBatch, trainOpts) ;

    % Move the CNN back to the CPU if it was trained on the GPU
    if numel(trainOpts.gpus) > 0
      net = vl_simplenn_move(net, 'cpu') ;
    end

    % Save the result for later use
    net.layers(end) = [] ;
    net.imageMean = mean_ex ;
    net.imageStd = std_ex;

%     save('data/face-experiment/facecnn.mat', '-struct', 'net') ;

    % Evaluate on the test data, also evaluate using correlation and RMSE
    load([location faceCheckersLoc_test(i).name]);

    % set a max value to the error
    errors(errors > 3) = 3;
    errors_test = errors;
    errors_q = quantizeContinuous(errors, 0, 3, num_bins);
    errors_test_q = errors_q;
    examples_test = single(examples);
    examples_test = bsxfun(@times, bsxfun(@minus, examples_test, net.imageMean), 1./net.imageStd);
    num_examples_test = size(examples_test,1);

    % keep adding rows
    while(size(mask,1) < 60)
        mask = cat(1, mask, false(1, size(mask,2)));
        triX = cat(1, triX, -ones(1, size(mask,2)));
    end

    % keep adding cols
    while(size(mask,2) < 60)
        mask = cat(2, mask, false(size(mask,1),1));
        triX = cat(2, triX, -ones(size(triX,1),1));
    end

    examples_r = single(zeros(size(mask, 1), size(mask, 2), num_examples_test));

    img_curr = zeros(size(mask));
    for e=1:num_examples_test

        img_curr(mask) = examples_test(e,:);
        examples_r(:, :, e) = img_curr;

    end
    examples_test = zeros(size(examples_r,1), size(examples_r,2), 1, size(examples_r,3));
    examples_test(:) = single(examples_r(:));

    ids = 1:99:num_examples_test;

    errors_all_test = zeros(num_examples_test,1);
    for k=1:numel(ids)-1
        examples_test_sm = single(examples_test(:,:,:,ids(k):ids(k+1)-1));
        res = vl_simplenn(net, examples_test_sm, [], []);
        res = gather(res(end).x);
        [~,res] = sort(res, 3, 'descend') ;
        res = squeeze(res);
        res = res(1,:);

        errors_test_rec = unQuantizeContinuous(res, 0, 3, num_bins)';
        errors_all_test(ids(k):ids(k+1)-1) = errors_test_rec;
    end
    errors_all_test = errors_all_test(1:ids(end));
    corrs_test = corr(errors_all_test, errors_test(1:numel(errors_all_test)));
    rmse_test = sqrt(mean((errors_all_test-errors_test(1:numel(errors_all_test))).^2));

    corrs_all = cat(1, corrs_all, corrs_test);
    rmses_all = cat(1, rmses_all, rmse_test);

    net.corrs = corrs_test;
    net.rmse = rmse_test;
    save([trainOpts.expDir, '/facecnn.mat'], '-struct', 'net');

    % The CNN
    faceCheckers(i).cnn = net;

    % The orientation
    faceCheckers(i).centres = centres;

    % The info for preprocessing
    faceCheckers(i).mask = mask;
    faceCheckers(i).nPix = nPix;
    faceCheckers(i).minX = minX;
    faceCheckers(i).minY = minY;

    faceCheckers(i).destination = shape(:,1:2);
    tri = load([location faceCheckersLoc_train(i).name], 'triangulation');
    faceCheckers(i).triangulation = tri.triangulation;
    faceCheckers(i).triX = triX;
    faceCheckers(i).mask = mask;
    faceCheckers(i).alphas = alphas;
    faceCheckers(i).betas = betas;

    faceCheckers(i).mean_ex = mean_ex;
    faceCheckers(i).std_ex = std_ex;

    WriteOutFaceCheckersCNNbinary("trained/validator_cnn.txt", faceCheckers);
end

save('trained/faceCheckers.mat', 'faceCheckers', 'corrs_all', 'rmses_all');

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

author: TadasBaltrusaitis

repo: OpenFace


text:
%CREATE TEMPLATES
%Letter
A=imread('letters_numbers\A.bmp');B=imread('letters_numbers\B.bmp');
C=imread('letters_numbers\C.bmp');D=imread('letters_numbers\D.bmp');
E=imread('letters_numbers\E.bmp');F=imread('letters_numbers\F.bmp');
G=imread('letters_numbers\G.bmp');H=imread('letters_numbers\H.bmp');
I=imread('letters_numbers\I.bmp');J=imread('letters_numbers\J.bmp');
K=imread('letters_numbers\K.bmp');L=imread('letters_numbers\L.bmp');
M=imread('letters_numbers\M.bmp');N=imread('letters_numbers\N.bmp');
O=imread('letters_numbers\O.bmp');P=imread('letters_numbers\P.bmp');
Q=imread('letters_numbers\Q.bmp');R=imread('letters_numbers\R.bmp');
S=imread('letters_numbers\S.bmp');T=imread('letters_numbers\T.bmp');
U=imread('letters_numbers\U.bmp');V=imread('letters_numbers\V.bmp');
W=imread('letters_numbers\W.bmp');X=imread('letters_numbers\X.bmp');
Y=imread('letters_numbers\Y.bmp');Z=imread('letters_numbers\Z.bmp');
%Number
one=imread('letters_numbers\1.bmp');  two=imread('letters_numbers\2.bmp');
three=imread('letters_numbers\3.bmp');four=imread('letters_numbers\4.bmp');
five=imread('letters_numbers\5.bmp'); six=imread('letters_numbers\6.bmp');
seven=imread('letters_numbers\7.bmp');eight=imread('letters_numbers\8.bmp');
nine=imread('letters_numbers\9.bmp'); zero=imread('letters_numbers\0.bmp');
%*-*-*-*-*-*-*-*-*-*-*-
letter=[A B C D E F G H I J K L M...
    N O P Q R S T U V W X Y Z];
number=[one two three four five...
    six seven eight nine zero];
character=[letter number];
templates=mat2cell(character,42,[24 24 24 24 24 24 24 ...
    24 24 24 24 24 24 24 ...
    24 24 24 24 24 24 24 ...
    24 24 24 24 24 24 24 ...
    24 24 24 24 24 24 24 24]);
save ('templates','templates')
clear all


text:
classdef EDNARMOEA < ALGORITHM
% <multi/many> <real/integer> <expensive>
% Efficient dropout neural network based AR-MOEA
% delta --- 0.05 --- Threshold of judging the diversity
% wmax  ---   20 --- Number of generations before updating Kriging models
% Ke    ---    3 --- Number of the solutions to be revaluated in each iteration

%------------------------------- Reference --------------------------------
% D. Guo, X. Wang, K. Gao, Y. Jin, J. Ding, and T. Chai. Evolutionary
% optimization of high-dimensional multiobjective and many-objective
% expensive problems assisted by a dropout neural network. IEEE
% Transactions on Systems, Man, and Cybernetics: Systems, 2022, 52(4):
% 2084-2097.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            assert(~isempty(ver('nnet')),'The execution of EDN-ARMOEA requires the Deep Learning Toolbox.');

            %% Parameter setting
            [delta,wmax,Ke] = Algorithm.ParameterSet(0.05,20,3);
            W     = UniformPoint(Problem.N,Problem.M);
            NI    = 11*Problem.D-1;
            P     = UniformPoint(NI,Problem.D,'Latin');
            A     = Problem.Evaluation(repmat(Problem.upper-Problem.lower,NI,1).*P+repmat(Problem.lower,NI,1));
            tr_x  = A.decs;
            tr_y  = A.objs;
            [tr_xx,ps] = mapminmax(tr_x');tr_xx=tr_xx';
            [tr_yy,qs] = mapminmax(tr_y');tr_yy=tr_yy';
            Params.ps  = ps;Params.qs=qs;
            RatioOld   = [];

            %% Train the model
            [net, Params] = trainmodel(tr_xx, tr_yy, Params);

            while Algorithm.NotTerminated(A)
                %% Update the model
                net=updatemodel(tr_xx, tr_yy, Params, net);

                %% Generate the sampling points and random population
                popsize=Problem.N;
                PopDec=repmat(Problem.upper-Problem.lower,popsize,1).*rand(popsize,Problem.D)+repmat(Problem.lower,popsize,1);
                [PopObj, PopMSE]=Estimate(PopDec, net, Params, Problem.M);
                [Archive,RefPoint,Range, Ratio] = UpdateRefPoint(PopObj,W,[]);
                if isempty(RatioOld)
                    RatioOld=Ratio;
                end

                %% Start the interations
                w=1;
                while w < wmax
                    MatingPool = MatingSelection(PopObj,RefPoint,Range);
                    OffspringDec  = OperatorGA(Problem,PopDec(MatingPool,:),{1,20,1,20});
                    [OffspringObj, OffspringMSE] = Estimate(OffspringDec, net, Params, Problem.M);
                    [Archive,RefPoint,Range, Ratio] = UpdateRefPoint([Archive;OffspringObj],W,Range);
                    MediatePopDec=[PopDec;OffspringDec];
                    MediatePopObj=[PopObj;OffspringObj];
                    MediatePopMSE=[PopMSE;OffspringMSE];
                    [Index,Range]       = EnvironmentalSelection(MediatePopObj,RefPoint,Range,popsize);
                    PopDec=MediatePopDec(Index,:);
                    PopObj=MediatePopObj(Index,:);
                    PopMSE=MediatePopMSE(Index,:);
                    w=w+1;
                end
                flag=RatioOld-Ratio<delta;
                PopNew=IndividualSelect(PopDec, PopObj, PopMSE, Ke, flag);
                RatioOld=Ratio;
                New = Problem.Evaluation(PopNew);
                A   = [A,New];
                [tr_x, tr_y]=SelectTrainData(A, 11*Problem.D-1, length(New));
                [tr_xx,ps]=mapminmax(tr_x');tr_xx=tr_xx';
                [tr_yy,qs]=mapminmax(tr_y');tr_yy=tr_yy';
                Params.ps=ps;Params.qs=qs;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function test_suite = test_moxunit_util_strjoin
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_strjoin_basics
    aeq=@(a,varargin)assertEqual(moxunit_util_strjoin(varargin{:}),a);

    aeq('a b cc',{'a','b','cc'});
    aeq('a>#<b>#<cc',{'a','b','cc'}, '>#<');
    aeq(sprintf('a\tb\tcc'),{'a','b','cc'}, '\t');
    aeq('a\b\cc',{'a','b','cc'}, '\\');
    aeq('a+b=cc',{'a','b','cc'}, {'+','='});

function test_strjoin_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        moxunit_util_strjoin(varargin{:}),'');

    % first argument must be cell with strings
    aet('a');
    aet(struct)
    aet(2);
    aet([]);
    aet({2});

    % second argument must be string or cell with strings
    aet({'a','b','c'},2);
    aet({'a','b','c'},struct);
    aet({'a','b','c'},[]);
    aet({'a','b','c'},[]);
    aet({'a','b','c'},{'x'});
    aet({'a','b','c'},{1,2});

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function y=polyValGen(p,x,direction)
%%POLYVALGEN Return the value of a scalar polynomial evaluated at the point
%       or points given in x. The polynomial can be specified in the format
%       y=p(1)*x^N+p(2)*x^(N-1)+...+p(N)*x+p(N+1)
%       or in the format
%       y=p(N+1)*x^N+p(N)*x^(N-1)+...+p(2)*x+p(1)
%       as specified by the direction option. This function is similar to
%       polyval, but allows one to choose between the formats without
%       calling flip to reverse the order of the elements in p.
%
%INPUTS: p An NX1 or 1XN set of real or complex polynomial coefficients;
%          wN>=1.
%        x A real or complex scalar point or a matrix of points at which
%          the scalar polynomial should be evaluated.
% direction An optional parameter specifying which of the two polynomial
%          formats is used. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) p(1) is
%            the coefficient of the highest order x term.
%          1 p(1) is the coefficient of the lowest order x term.
%
%OUTPUTS: y The value or values of the polynomial in p evaluated at the
%           point(s) in x. This is the same size as x, unless x is an empty
%           matrix, in which case y is 0.
%
%The function is evaluated using Horner's rule, which is described in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Horner's Rule." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/HornersRule.html
%
%July 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(direction))
   direction=0;
end

numP=length(p);

y=zeros(size(x));
if(direction==0)
    if(numP>0)
        y(:)=p(1);
    end

    for k=2:numP
        y=x.*y+p(k);
    end
else
    if(numP>0)
        y(:)=p(numP);
    end

    for k=(numP-1):-1:1
        y=x.*y+p(k);
    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:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: poisspdf_of_lambda.m
% Used when you want to integrate a poisson pdf with
% respect to lambda.
%
% 022300 tdr created
% 030400 tdr replaced global variable approach to getting x and n to binopdf.
% 031600 tdr converted from binopdf_of_p.m to poisson dist.
% 032200 tdr replaced for loop with vector argument.

function y = poisspdf_of_lambda(lambda, x)

if nargin ~= 2,
    error('Requires two input arguments (lambda,x)');
end

y = poisspdf_0ok(x,lambda);