text:
%% flandmark_compilemex.m
% Compilation script for flandmark mex-files.
% Tested with OpenCV version 2.3 (linux, )
%
% 09-12-12, Michal Uricar

clc;
clearvars; close all;

% Change this to your version of opencv
OpenCVmajor = 2;
OpenCVminor = 4;

% Setup paths accpording to system platform
if (isunix&&~ismac)
    fprintf('Compiling mex-files on Unix platform...\n');

    include = ' -I/usr/include/opencv/ -I/usr/include/ ';
    libpath = '/usr/lib/';
    libflandmark = ' -L../bin/libflandmark/ -lflandmark_static ';

    if (OpenCVminor < 3)
        cxcore = dir([libpath 'libcxcore.so.' num2str(OpenCVmajor) '.' num2str(OpenCVminor)]);
        cxcore = cxcore.name;
        cv = dir([libpath 'libcv.so.' num2str(OpenCVmajor) '.' num2str(OpenCVminor)]);
        cv = cv.name;
        cvaux = dir([libpath 'libcvaux.so.' num2str(OpenCVmajor) '.' num2str(OpenCVminor)]);
        cvaux = cvaux.name;

        cvlibs = [libpath cxcore ' ' libpath cv ' ' libpath cvaux ' '];
    else
        opencv_core = dir([libpath 'libopencv_core.so.' num2str(OpenCVmajor) '.' num2str(OpenCVminor)]);
        opencv_core = opencv_core.name;
        opencv_imgproc = dir([libpath 'libopencv_imgproc.so.' num2str(OpenCVmajor) '.' num2str(OpenCVminor)]);
        opencv_imgproc = opencv_imgproc.name;

        cvlibs = [libpath opencv_core ' ' libpath opencv_imgproc ' '];
    end;

    % compile mexfiles
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_detector']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_base_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_detector_base']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_write_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_write_model']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_load_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_load_model']);
    eval(['mex -O -largeArrayDims ./mex/lbp2sparse_mex.cpp -I../libflandmark/' libflandmark ' -o ../matlab_toolbox/mex/lbp2sparse']);
    eval(['mex -O -largeArrayDims ./mex/argmax_mex.cpp -I../libflandmark/ ' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/argmax_mex']);
    eval(['mex -O -largeArrayDims ./mex/getNormalizedFrame_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/getNormalizedFrame']);

    fprintf('Compilation finished.\n');
end;

if (ispc)
    fprintf('Compiling mex-files on Windows platform...\n');

    include = [' -Ic:\OpenCV' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '\build\include\opencv\', ...
               ' -Ic:\OpenCV' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '\build\include\' ];
    libpath = ['c:\OpenCV' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '\build\x64\vc10\lib\'];
    libflandmark = ' -L../bin7/libflandmark/Release -lflandmark_static ';

    if (OpenCVminor < 3)
        cxcore = dir([libpath 'cxcore' num2str(OpenCVmajor) num2str(OpenCVminor) '*.lib']);
        cxcore = cxcore.name;
        cv = dir([libpath 'cv' num2str(OpenCVmajor) num2str(OpenCVminor) '*.lib']);
        cv = cv.name;
        cvaux = dir([libpath 'cvaux' num2str(OpenCVmajor) num2str(OpenCVminor) '*.lib']);
        cvaux = cvaux.name;

        cvlibs = [libpath cxcore ' ' libpath cv ' ' libpath cvaux ' '];
    else
        opencv_core = dir([libpath 'opencv_core' num2str(OpenCVmajor) num2str(OpenCVminor) '*.lib']);
        opencv_core = opencv_core.name;
        opencv_imgproc = dir([libpath 'opencv_imgproc' num2str(OpenCVmajor) num2str(OpenCVminor) '*.lib']);
        opencv_imgproc = opencv_imgproc.name;

        cvlibs = [libpath opencv_core ' ' libpath opencv_imgproc ' '];
    end;

    % compile mexfiles
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/flandmark_detector']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_base_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/flandmark_detector_base']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_write_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/flandmark_write_model']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_load_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/flandmark_load_model']);
    eval(['mex -O -largeArrayDims ./mex/lbp2sparse_mex.cpp -I../libflandmark/' libflandmark ' -output ../matlab_toolbox/mex/lbp2sparse']);
    eval(['mex -O -largeArrayDims ./mex/argmax_mex.cpp -I../libflandmark/ ' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/argmax_mex']);
    eval(['mex -O -largeArrayDims ./mex/getNormalizedFrame_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/getNormalizedFrame']);

    fprintf('Compilation finished.\n');
end;

if (ismac)
    fprintf('Compiling mex-files on Mac platform...\n');

    include = ' -I/opt/local/include/opencv/ -I/opt/local/include/ -I/usr/local/include/opencv/ -I/usr/local/include/';
    if(numel(dir('/usr/local/lib/libopencv*')))
        libpath = '/usr/local/lib/';
    else
        libpath = '/opt/local/lib/';
    end
    libflandmark = ' -L../bin/libflandmark/ -lflandmark_static ';

    % files = dir([libpath 'libopencv*.dylib']);
    % e.g. libopencv_core.2.4.2.dylib

    if (OpenCVminor < 3)
        cxcore = dir([libpath 'libcxcore.' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '.dylib']);
        cxcore = cxcore.name;
        cv = dir([libpath 'libcv.' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '.dylib']);
        cv = cv.name;
        cvaux = dir([libpath 'libcvaux.' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '.dylib']);
        cvaux = cvaux.name;

        cvlibs = [libpath cxcore ' ' libpath cv ' ' libpath cvaux ' '];
    else
        opencv_core = dir([libpath 'libopencv_core.' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '.dylib']);
        opencv_core = opencv_core.name;
        opencv_imgproc = dir([libpath 'libopencv_imgproc.' num2str(OpenCVmajor) '.' num2str(OpenCVminor) '.dylib']);
        opencv_imgproc = opencv_imgproc.name;

        cvlibs = [libpath opencv_core ' ' libpath opencv_imgproc ' '];
    end;

    % compile mexfiles
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_detector']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_detector_base_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_detector_base']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_write_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_write_model']);
    eval(['mex -O -largeArrayDims ./mex/flandmark_load_model_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/flandmark_load_model']);
    eval(['mex -O -largeArrayDims ./mex/lbp2sparse_mex.cpp -I../libflandmark/' libflandmark ' -o ../matlab_toolbox/mex/lbp2sparse']);
    eval(['mex -O -largeArrayDims ./mex/argmax_mex.cpp -I../libflandmark/ ' include libflandmark cvlibs ' -output ../matlab_toolbox/mex/argmax_mex']);
    eval(['mex -O -largeArrayDims ./mex/getNormalizedFrame_mex.cpp -I../libflandmark/' include libflandmark cvlibs ' -o ../matlab_toolbox/mex/getNormalizedFrame']);

    fprintf('Compilation finished.\n');
end;

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function W=RegCsp(EEGdata,gnd,genSs,genMs,beta,gamma)
% RCSP: Regularized Common Spatial Pattern
%
% %[Prototype]%
% function W=RegCsp(EEGdata,gnd,genSs,genMs,beta,gamma)
%
% %[Author Notes]%
% Author: Haiping LU
% Email : hplu@ieee.org   or   eehplu@gmail.com
% Release date: March 20, 2012 (Version 1.0)
% Please email me if you have any problem, question or suggestion
%
% %[Algorithm]%:
% This function implements the Regularized Common Spatial Pattern
% (R-CSP) algorithm presented in the follwing paper:
%    Haiping Lu, How-Lung Eng, Cuntai Guan, K.N. Plataniotis, and A.N. Venetsanopoulos,
%    "Regularized Common Spatial Pattern With Aggregation for EEG Classification in Small-Sample Setting",
%    IEEE Trans. on Biomedical Engineering, Vol. 57, No. 12, Pages
%    2936-2946, Dec. 2010.
% Please reference this paper when reporting work done using this code.
%
% The following is an earlier conference version
%    Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos,
%    "Regularized Common Spatial Patterns with Generic Learning for EEG Signal Classification",
%    in Proceedings of the 31st Annual International Conference of the
%    IEEE Engineering in Medicine and Biology Society (EMBC), Sep., 2009.
%
% %[Syntax]%: W=RegCsp(EEGdata,gnd,genSs,genMs,beta,gamma)
%
% %[Inputs]%:
%    EEGdata: the input EEG data for training, a 3D array with size [numT,numCh,nTrl]
%       numT is the number of samples in each channle (T in the paper)
%       numCh is the number of channels (N in the paper)
%       nTrl is the number of trails for training (2*M in the paper)
%
%    gnd: the ground truth class labels for the nTrl trials, size nTrl x 1
%
%    genSs: size numCh x numCh x 2, the sum of generic covariance matrices
%       from generic training trials for two classes, equation (6) in the paper
%
%    genMs: size 2 x 1, the number of generic training trials for two classes
%
%    beta,gamma: the reguarlizatio parameters in the paper, setting both to
%       zero gives the conventional CSP
%
% %[Outputs]%:
%    W: \hat{W} in the paper, the projection matrix, we often use only the
%         first 6 columns from W
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %[Notes]%:
% A. Developed using Matlab R2006a
% B. Revision history:
%       Version 1.0 released on March 20, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%
[numT,numCh,nTrl] = size(EEGdata);%Input data, see documentation
Ns=zeros(2,1);
Sigma2=cell(2,1);%Two average spatial covariance matrices, one for each class
%%%%%%%%%%%%%%%%%%%
for i=1:2%for each class
    Idxs=find(gnd==i);
    EEG=EEGdata(:,:,Idxs);
    Ns(i)=length(Idxs);
    C=zeros(numCh,numCh,Ns(i));%Sample covariance matrix, equation (1) in the paper
    for trial=1:Ns(i)
        E=EEG(:,:,trial)';
        tmpC = (E*E');
        C(:,:,trial) = tmpC./trace(tmpC);%normalization
    end
    Csum=sum(C,3);
    %Reguarlization, see equations (3) and (4) in the paper
    Sigma1=((1-beta)*Csum+beta*genSs(:,:,i))/((1-beta)*Ns(i)+beta*genMs(i));
    Sigma2{i}=(1-gamma)*Sigma1+gamma*trace(Sigma1)*eye(numCh)/numCh;
end

%Equation (7) in the paper
SigmaComps=Sigma2{1}+Sigma2{2};
[Ucomps,lmds] = eig(SigmaComps);
[lmds,Idxs] = sort(diag(lmds),'descend');
Ucomps = Ucomps(:,Idxs);

%Note equations (8) to (12) in the conference version is now condensed in
%one equation (8) in the journal version
%Equation (8) in the CONFERENCE paper
P=sqrt(inv(diag(lmds)))*Ucomps';

%Equation (9) in CONFERENCE the paper
Sgm1=P*Sigma2{1}*P';

%Equation (11) in the CONFERENCE paper
[B,D] = eig(Sgm1);
[D,Idxs] = sort(diag(D),'descend');
B = B(:,Idxs);

%Equation (12) in the CONFERENCE paper
%Equation (8) in the JOURNAL paper
W=(B'*P);
%Normalize the projrection matrix
for i=1:length(Idxs), W(i,:)=W(i,:)./norm(W(i,:)); end

%Sort columns, take first and last columns first, etc
W0=W;
W=zeros(size(W));
i=0;
for d=1:numCh
    if (mod(d,2)==0)
        W(d,:)=W0(numCh-i,:);
        i=i+1;
    else
        W(d,:)=W0(1+i,:);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [a]=minus(b,c)
%A=B-C
%   [A]=MINUS(B,C) Subtraction of two TT-matrices, B and C
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
a=b+(-1.0)*c;
return
end

author: oseledets

repo: TT-Toolbox


text:
% erpregoutfunc() - sub function of erpregout() used to regress
%                   out the ERP from the data
%
% Usage:
%   totdiff = erpregout(fact, data, erp);
%
% Inputs:
%   fact    - factor
%   data    - [float] 1-D data (time points).
%   erp     - [float] 1-D data (time points).
%
% Outputs:
%   totdif  - residual difference
%
% Author: Arnaud Delorme, Salk, SCCN, UCSD, CA, April 29, 2004

%123456789012345678901234567890123456789012345678901234567890123456789012

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

function totdiff = erpregoutfunc(fact, data, erp);

    totdiff = mean(abs(data - fact*erp));

author: PatternRecognition

repo: OpenBMI


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [yc,dy,para] = getTrafoFromVelocityRK4(vc,yc,varargin)
%
% compute transformation yc by integrating velocity field in time using
% a fourth-order Runge-Kutta method. Here, vc is assumed to be stationary.
% For instationary velocities, use getTrafoFromInstatinaryVelocityRK4.m
%
% For more details see Sec. 3 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% Input:
%
%  vc   - discrete velocity field (nodal, cell-centered, or staggered)
%  yc   - particle positions
%
% Additional REQUIRED Input (provided through varargin)
%
%  omega        - spatial domain (required!)
%  m            - number of cells in each direction (required!)
%
% Optional Input (provided through varargin)
%
%  doDerivative - compute derivative w.r.t. velocity
%  tspan        - time interval (default: [0 1])
%  N            - number of time discretization points (nodal) (default: 5)
%  storeInter   - store intermediate transformation (e.g., for visualization)
%
% Output:
%
%  yc           - end point of characteristics
%  dy           - derivative w.r.t. vc
%  para         - info such as CFL
%
% =========================================================================
function [yc,dy,para] = getTrafoFromVelocityRK4(vc,yc,varargin)

if nargin==0
    runMinimalExample
    return;
end
doDerivative = (nargout>1);
omega = [];
m     = [];
tspan = [0 1];
N     = 5;
storeInter = false;
for k=1:2:length(varargin)     % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;
if isempty(omega) || isempty(m)
    error('%s - omega and m must be provided through varargin or trafo(''set'',''omega'',omega,''m'',m')
end
dt   = (tspan(2)-tspan(1))/(N-1); % note dt is negatave when going backwards in time
dy = [];
if doDerivative
    dy = sparse(numel(yc),numel(vc));
end
% compute the CFL number (not actually required for stability in Lagrangian
% methods, but an indicator of how many voxels the particles may move in
% one time step.
h    = (omega(2:2:end)-omega(1:2:end))./m;
% vt   =  sqrt(sum(reshape(vc.^2,[],dim),2));
% CFL  = max(vt)*dt/min(h);
CFL = 0.0;
para = struct('CFL',CFL,'dt',dt,'N',N,'omega',omega,'m',m,'h',h);
if storeInter
    para.YC = zeros(numel(yc),N);
    para.YC(:,1) = yc;
end
for k=1:N-1
    [vi,dvidy] = linearInterGrid(vc,omega,m,yc,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yc);
        dyi = Ty + dvidy*dy;
        dytemp = dyi;
    end
    ytemp = vi;
    yi    = yc + .5*dt*vi;

    [vi,dvidy] = linearInterGrid(vc,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        dyi = Ty + dvidy*(dy+.5*dt*dyi);
        dytemp = dytemp + 2*dyi;
    end
    ytemp = ytemp + 2*vi;
    yi    = yc + .5*dt*vi;

    [vi,dvidy] = linearInterGrid(vc,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        dyi = Ty + dvidy*(dy+.5*dt*dyi);
        dytemp = dytemp + 2*dyi;
    end
    ytemp = ytemp + 2*vi;
    yi    = yc + dt*vi;
    [vi,dvidy] = linearInterGrid(vc,omega,m,yi,'doDerivative',doDerivative);
    if doDerivative
        Ty = getLinearInterGridMatrix(omega,m,yi);
        dyi = Ty + dvidy*(dy+dt*dyi);
        dytemp = dytemp + dyi;
        dy  = dy + (dt/6)*dytemp;
    end
    ytemp = ytemp + vi;
    yc    = yc + (dt/6)*ytemp;
    if storeInter, para.YC(:,k+1) = yc; end
end

function runMinimalExample

omegaV = [-1 1 -1 1];
omega = .1*omegaV;
m     = [24 24];
tspan     = [2 3];
N     = 2;
yc    = getNodalGrid(omega,m)+.0002;

% test cell-centered grid
regularizer('set','regularizer','mbCurvature','alpha',1)
xc      = reshape(getCellCenteredGrid(omegaV,m),[],2);
vc      = [.3*sign(xc(:,1)).*xc(:,1).^2 sin(pi*xc(:,2))];
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))

% test nodal grid
regularizer('set','regularizer','mbElasticNodal','alpha',1)
xc      = reshape(getNodalGrid(omegaV,m),[],2);
vc      = [.3*sign(xc(:,1)).*xc(:,1).^2 sin(pi*xc(:,2))];
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))

% test staggered grid
regularizer('set','regularizer','mbElastic','alpha',1)
vc  = grid2grid(vc(:),m,'nodal','staggered');
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))

omegaV = [-1 1 -1 1 -1 1];
omega = .1*omegaV;
m     = [16 16 8];
tspan = [4 3];
N     = 10;
yc    = getNodalGrid(omega,m)+.0002;

regularizer('set','regularizer','mbCurvature','alpha',1)
xc      = reshape(getCellCenteredGrid(omegaV,m),[],3);
vc      = [.3*sign(xc(:,1)).*xc(:,1).^2 sin(pi*xc(:,2)) xc(:,3)];
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
% [yc,dy] = fctn(vc(:));
checkDerivative(fctn,0*vc(:))

% test nodal grid
regularizer('set','regularizer','mbElasticNodal','alpha',1)
xc      = reshape(getNodalGrid(omegaV,m),[],3);
vc      = [.3*sign(xc(:,1)).*xc(:,1).^2 sin(pi*xc(:,2)) xc(:,3)];
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))

% test staggered grid
regularizer('set','regularizer','mbElastic','alpha',1)
vc  = grid2grid(vc(:),m,'nodal','staggered');
fctn    = @(vc) getTrafoFromVelocityRK4(vc,yc,'omega',omegaV,'m',m,'T',tspan,'N',N);
checkDerivative(fctn,vc(:))


text:
%% (Internal) Display an enumeration of options to a string
%
%     This function format a string for fprintf like functions in order to
%     present a message srtMessage followed by an enumeration of string
%     elements included in the cell str cOptionEnum
%
%   str_aux = disp_option_enumeration(srtMessage, cOptionEnum)
%
% Arguments:
%
%      + srtMessage: A title message previous to the enumeration
%
%      + cOptionEnum: cellstring with the items to enumerate
%
% Output:
%
%      + str_aux: resulting string
%
% Example:
%
%
% See also disp_string_framed, disp_string_title
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function str_aux = disp_option_enumeration(srtMessage, cOptionEnum)

%  Escape problematic strings, such as file paths
cOptionEnum = cellfun(@(a)(strrep(a, '\', '\\')), cOptionEnum, 'UniformOutput', false );

str_aux = [ repmat(' + ', length(cOptionEnum), 1) char(cOptionEnum) repmat('\n', length(cOptionEnum), 1 ) ];
str_aux = [ srtMessage '\n' rowvec(str_aux') ];

author: marianux

repo: ecg-kit

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


text:
function s=get_name(a,detailed)
s=[get_name(a.algorithm)];
if length(a.child)==1,
    s=[s '(' get_name(a.child{1}) ')'];
else
   s=[s '(' get_name(a.child{1}) '...)'];
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [obj] = update(obj)
    % Updates the NLP problems before load it to NLP solver
    %
    % @note Overload the superclass (NonlinearProgram) for additional
    % updates

    % initialize the type of the variables
    obj.VariableArray = NlpVariable.empty();
    obj.CostArray  = NlpFunction.empty();
    obj.ConstrArray = NlpFunction.empty();
    phase = obj.Phase;
    n_phase = length(phase);

    index_s = 1;
    % register variables
    for i=1:n_phase

        % register variables
        obj = regVariable(obj, phase(i).OptVarTable);

        index_f = sum([obj.VariableArray.Dimension]);
        obj.PhaseVarIndices(i,:)= [index_s, index_f];
        index_s = index_f+1;

        % register constraints
        obj = regConstraint(obj, phase(i).ConstrTable);

        % register cost functions
        obj = regObjective(obj, phase(i).CostTable);

    end

    obj = update@NonlinearProgram(obj);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function solution = gnc_category_registration(problem,SDP,path,varargin)
%% Solve outlier-robust category registration using GNC
%% Heng Yang, July 06, 2021

params = inputParser;
params.CaseSensitive = false;
params.addParameter('lambda',0.1, @(x) isscalar(x));
params.addParameter('refine',false, @(x) islogical(x));
params.addParameter('refineversion','v2', @(x) ischar(x));
params.parse(varargin{:});
lambda  = params.Results.lambda;
refine  = params.Results.refine;
refineversion = params.Results.refineversion;

N       = problem.N;
barc2   = 1.0;
shapes  = problem.shapes;
scene   = problem.scene;
noiseBoundSq = problem.noiseBoundSq;
cBound  = problem.cBound;
tBound  = problem.translationBound;
t0      = tic;
fprintf('\n\n===================================== GNC-TLS =======================================\n')

allPoints = 1:N;
weights   = ones(N,1);
stopTh    = 1e-20;
maxSteps  = 1e2;
divFactor = 2;
itr       = 0;
pre_f_cost= 1e6;
cost_diff = 1e6;
breakyes  = 0;

fprintf('epsilon: %3.2e, div: %1.1f, maxiters: %3.0e.\n',stopTh,divFactor,maxSteps);
fprintf('-------------------------------------------------------------------------------------\n')
fprintf(' itr |    obj        delobj  |    mu    |   sumw   |  sumout  |   maxres   |   gap   |\n')

while (true)
    if max(abs(weights)) < 1e-6
        msg      = 'GNC encounters numerical issues, the solution is likely to be wrong.';
        breakyes = 3;
    end
    if itr == maxSteps
        breakyes = 2;
        msg      = 'Maximum iterations reached.';
    end
    if cost_diff < stopTh
        breakyes = 1;
        msg      = sprintf('GNC converged %3.2e < %3.2e.',cost_diff,stopTh);
    end
    if breakyes > 0
        fprintf('%s\n',msg);
        break
    end

    problem.weights     = weights;
    [R,t,c,ofinfo]      = outlier_free_category_registration(problem,path,'lambda',lambda);
    shape               = combine_shapes(shapes,c);
    residuals           = sum( (scene - R*shape - t).^2, 1) / noiseBoundSq;
    f_cost              = sum(min(residuals,barc2)) + lambda * (c'*c);

    maxResidual = max(residuals);
    if itr < 1
        mu = max(1 / ( 5 * maxResidual / barc2 - 1 ),1e-6); % make sure mu is positive
    end

    % update weights in closed-form
    th1 = (mu+1)/mu * barc2;
    th2 = (mu)/(mu+1) * barc2; % th1 > th2
    for i = 1:N
        if residuals(i) - th1 >= 0
            weights(i) = 0;
        elseif residuals(i) - th2 <= 0
            weights(i) = 1;
        else
            weights(i) = sqrt( barc2*mu*(mu+1)/residuals(i) ) - mu;
%             assert(weights(i)>= 0 && weights(i) <=1, 'weights calculation wrong!');
        end
    end
    cost_diff   = abs(f_cost - pre_f_cost);
    fprintf('%4d | %3.4e %3.4e | %3.2e | %3.2e | %3.2e | %3.4e | %1.1e |\n',...
            itr,f_cost,cost_diff,mu,sum(weights),N-sum(weights),maxResidual,ofinfo.gap);

    % increase mu and compute cost difference
    mu          = mu * divFactor;
    itr         = itr + 1;
    pre_f_cost  = f_cost;
end
R_est = R;
t_est = t;
c_est = c;

theta_est               = zeros(N,1);
theta_est(weights>0.5)  = 1;
theta_est(weights<0.5)  = -1;

% if max(max(-c_est,0)) > 0
%     c_est(c_est <= 0)       = 1e-3;
%     addpath(genpath(path.manoptpath))
%     rrPar.blk               = SDP.blk;
%     rrPar.translationBound  = problem.translationBound;
%     rrPar.cBound            = problem.cBound;
%     rrPar.N                 = problem.N;
%     rrPar.K                 = problem.K;
%
%     [R_est,t_est]  = invert_transformation(R_est,t_est);
%     [~,~,out]      = nlp_catreg_v2(SDP.C,rrPar,R_est,t_est,c_est,theta_est);
%     [out.R,out.t]  = invert_transformation(out.R,out.t);
%     rmpath(genpath(path.manoptpath))
%
%     R_est          = out.R;
%     t_est          = out.t;
%     c_est          = out.c;
%     theta_est      = out.theta;
%
% end

% make sure c_est is norm bounded
nncflag = 0;
if sum(c_est <= 0) > 0
    nncflag           = 1;
    fprintf('GNC c has %d entries below 0.\n',sum(c_est <= 0));
    c_est(c_est <= 0) = 0;
end
bdcflag = 0;
if norm(c_est) >= cBound
    bdcflag           = 1;
    fprintf('GNC c has norm %g > %g.\n',norm(c_est),cBound);
    c_est = c_est / norm(c_est) * cBound;
end
% make sure t_est is norm bounded
bdtflag = 0;
if norm(t_est) >= tBound
    bdtflag     = 1;
    fprintf('GNC t has norm %g > %g.\n',norm(t_est),tBound);
    t_est = t_est / norm(t_est) * tBound;
end

% re-evaluate the cost
shape               = combine_shapes(shapes,c_est);
residuals           = sum( (scene - R_est*shape - t_est).^2, 1) / noiseBoundSq;
f_est               = sum(min(residuals,barc2)) + lambda * (c_est'*c_est);

% if refine
%     addpath(genpath(path.manoptpath))
%     rrPar.blk               = SDP.blk;
%     rrPar.translationBound  = problem.translationBound;
%     rrPar.cBound            = problem.cBound;
%     rrPar.N                 = problem.N;
%     rrPar.K                 = problem.K;
%     switch refineversion
%     case 'v1'
%         [~,fopt,out] = nlp_catreg(SDP.C,rrPar,R_est,t_est,c_est,theta_est);
%     case 'v2'
%         [R_est,t_est] = invert_transformation(R_est,t_est);
%         [~,fopt,out]  = nlp_catreg_v2(SDP.C,rrPar,R_est,t_est,c_est,theta_est);
%         [out.R,out.t] = invert_transformation(out.R,out.t);
%     otherwise
%         error('Unknown refine version.')
%     end
%     if fopt < f_est
%         fprintf('        MANOPT cost %3.4e < GNC cost %3.4e.\n',fopt,f_est);
%         R_est   = out.R;
%         t_est   = out.t;
%         c_est   = out.c;
%         theta_est = out.theta;
%         f_est   = fopt;
%     end
%     rmpath(genpath(path.manoptpath))
% end

time_gnc = toc(t0);

solution.type           = 'GNC-TLS';
solution.weights        = weights;
solution.theta_est      = theta_est;
solution.R_est          = R_est;
solution.t_est          = t_est;
solution.c_est          = c_est;
solution.itr            = itr;
solution.divFactor      = divFactor;
solution.time           = time_gnc;
solution.f_est          = f_est;
solution.residuals      = residuals;
solution.detectedOutliers = allPoints(theta_est<0);
solution.nncflag        = nncflag;
solution.bdcflag        = bdcflag;
solution.bdtflag        = bdtflag;

% print some info
fprintf('f_est = %g, divFactor=%g, itr=%d, time_gnc=%g[s].\n',f_est,divFactor,solution.itr,time_gnc);
fprintf('=====================================================================================\n\n\n')
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function [mod_ref_data, mod_deg_data]= input_filter( ref_data, ref_Nsamples, ...
    deg_data, deg_Nsamples)

mod_ref_data= DC_block( ref_data, ref_Nsamples);
mod_deg_data= DC_block( deg_data, deg_Nsamples);

mod_ref_data= apply_filters( mod_ref_data, ref_Nsamples);
mod_deg_data= apply_filters( mod_deg_data, deg_Nsamples);

author: vipchengrui

repo: traditional-speech-enhancement


text:
function test (varargin)
% TODO: in stack information, replace pointers to temporary files with
%       pointes to test section in the mfile being tested.
% TODO: add test case for get_stack_trace
% TODO: add test case for functions other than m-files that have
%       an accompanying homonymous m-file containing the test blocks.
% TODO: support subfunctions. (Nested functions will never be supported
%!      because they cannot stand alone.)
% TODO: is seems that a test block, at the end of a file, without a blank
%       line after it, is not supported.

    %%%%%%%%%%
    switch nargin
    case 0
        error('test:notEnoughInputs', ...
        'Not enough input arguments.');
    case 1
        mfunc = varargin{1};
        param = '';

        mfile = which(mfunc);
        if isempty(mfile)
            error('test:UndefinedFunction', ...
            'Function %s not found in path.', mfunc);
        end
        [mdir, mfunc, ext] = fileparts(mfile);
        if ~strcmp(ext, '.m')
            mfile = fullfile(mdir, [mfunc '.m']);
            if ~exist(mfile, 'file')
                error('test:FileNotFound',...
                'File not found: %s', mfile);
            end
        end

        %[tdir, tfunc] = fileparts(tempname);
        %tfile = fullfile(tdir, [tfunc '.m']);
        % under Vista, it's significantly faster to addpath
        % an empty directory rather than the standard temporary dir.
        tdir = tempname;
        tfunc = 'testit';
        tfile = fullfile(tdir, [tfunc '.m']);

        pdir = tdir;
        create_pdir = true;

        test_type = 'func';
    case 2
        mfunc = varargin{1};
        param = varargin{2};

        eval(sprintf('%s;', param));
        mfunc2 = sprintf('%s(%s)', mfunc, param);
        mfile = which(mfunc2);
        if isempty(mfile)
            error('test:UndefinedFunction',...
            'Function %s not found in path.', mfunc2);
        end
        [mdir, mfunc, ext] = fileparts(mfile);
        if ~strcmp(ext, '.m')
            mfile = fullfile(mdir, [mfunc '.m']);
            if ~exist(mfile, 'file')
                error('test:FileNotFound',...
                'File not found: %s', mfile);
            end
        end

        %[tdir, tfunc] = fileparts(tempname);
        %tfile = fullfile(tdir, [tfunc '.m']);
        % under Vista, it's significantly faster to addpath
        % an empty directory rather than the standard temporary dir.
        tdir = tempname;
        tfunc = 'testit';
        tfile = fullfile(tdir, [tfunc '.m']);

        pdir = tdir;
        create_pdir = true;

        test_type = 'method';
    case 3
        mfunc = varargin{1};
        param = char(varargin{2});  % so that printf works when param = []
        mdir  = varargin{3};

        if ~isempty(param) && strcmp(mdir, 'private')
            mdir = [fileparts(which(param)) '\private'];
        end
        if ~(   strend([filesep 'private' filesep], mdir) ...
             || strend([filesep 'private']        , mdir) )
            error('test:DirEndInPrivate',...
            ['Third parameter %s should end in '...
            '"%cprivate".'], mdir, filesep);
        end

        mfile = fullfile(mdir, [mfunc '.m']);
        if ~exist(mfile, 'file')
            error('test:FileNotFound',...
            'File not found: %s', mfile);
        end

        tdir = dirup(mdir);
        create_pdir = false;
        while true
            [temp, tfunc] = fileparts(tempname); %#ok<ASGLU>
            tfile = fullfile(tdir, [tfunc '.m']);
            if ~exist(tfile, 'file'),  break;  end
        end

        if isempty(param)
            pdir = tdir;
            test_type = 'priv_func';
        else
            if ~strend([filesep '@' param], tdir)
                error('test:DirEndInClass',...
                ['Third parameter %s should end in '...
                '"%c@%s%cprivate".'], mdir, filesep, param, filesep);
            end
            pdir = dirup(tdir);
            test_type = 'priv_method';
        end
    otherwise
        error('test:tooManyInputs',...
        'Too many input parameters.');
    end
    %mdir, tdir, pdir  % DEBUG

    %%%%%%%%%%
    %mfile  % DEBUG
    %type(mfile)  % DEBUG
    [mfid, msg] = fopen (mfile, 'rt');
    if (mfid == -1)
        error('test:FileOpening', msg);
    end

    % extract blocks (test, error, shared) in m-file
    block_start = [];
    is_in_block = false;
    i = 1;  % # of %! lines
    j = 1;  % # of test or error blocks
    mline = fgets(mfid);
    while ~isequal(mline, -1)  % while not end of file
        if (   strncmp(mline, '%!test', 6) ...
            || strncmp(mline, '%!error', 7) ...
            || strncmp(mline, '%!shared', 8) )
            % start block
            is_in_block = true;
            block_start(j) = i;
            if strncmp(mline, '%!test', 6)
                block_type{j} = 'test';
            elseif strncmp(mline, '%!error', 7)
                block_type{j} = 'error';
            elseif strncmp(mline, '%!shared', 8)
                block_type{j} = 'shared';
            end
        elseif is_in_block && strncmp(mline, '%!', 2)
            % read block contents
            tline{i} = mline(4:end);
            i = i + 1;
        elseif is_in_block
            % end block
            block_end(j) = i-1;
            j = j + 1;
            is_in_block = false;
        end
        mline = fgets(mfid);
    end
    if is_in_block
        % end latest block
        block_end(j) = i-1;
    end
    fclose(mfid);

    if isempty(block_start) || all(strcmp(block_type, 'shared'))
        error('test:noBlock',...
        'Found no test or error section in file %s.', mfile);
    end
    myassert (length(block_start), length(block_end));
    myassert (length(block_start), length(block_type));

    %%%%%%%%%%

    %% create test file
    old_path = path;
    if create_pdir,  mkdir(pdir);  end
    addpath(pdir);
    %addpath(pdir, '-frozen');

    %% write and run each block
    for j=1:length(block_start)
        if strcmp(block_type{j}, 'shared'),  continue;  end

        [tfid, msg] = fopen(tfile, 'wt');
            if (tfid == -1),  error('test:FileOpening', msg);  end

        if strcmp(test_type, 'priv_method')
            fprintf(tfid, 'function %s (obj)\n', tfunc);
        else
            fprintf(tfid, 'function %s\n', tfunc);
        end

        found_mfunc_call = false;

        %% Write optional, previous, shared block.
        k = find (strcmp({block_type{1:j}}, 'shared'), 1, 'last');
        if ~isempty(k)
            myassert (strcmp(block_type{k}, 'shared'));
            fprintf(tfid, '%% shared\n');
            for i=block_start(k):block_end(k)
                fprintf(tfid, '%s', tline{i});
                found_mfunc_call = found_mfunc_call ...
                                   || is_function_called (tline{i}, mfunc);
            end
            fprintf(tfid, '\n');
        end

        %% Is it an error block?
        if strcmp(block_type{j}, 'error')
            fprintf(tfid, 'try\n');
            % In an error block, error is the expected behaviour.
            % We do we following:
            %     try,
            %         % INSERT USER CODE HERE
            %     catch,
            %         return;
            %     end
            %     error('test:error:fail', '');
            % If the user code generates an error, it is catched and
            % the test runs OK. If the user code does _not_ generate
            % an error(contrary to what was expected), then the test fails.
        end

        for i=block_start(j):block_end(j)
            fprintf(tfid, '%s', tline{i});
            found_mfunc_call = found_mfunc_call ...
                               || is_function_called (tline{i}, mfunc);
        end

        if (~found_mfunc_call)
            warning('test:noFuncCall', ...
                'Found no call to function %s itself in its %s block #%d.',...
                mfunc, block_type{j}, j);
        end

        if strcmp(block_type{j}, 'error')
            fprintf(tfid, 'catch\n');
            fprintf(tfid, '    return;\n');
            fprintf(tfid, 'end\n');
            fprintf(tfid, 'clear s;\n');
            fprintf(tfid, 's.identifier=''test:error:fail'';\n');
            fprintf(tfid, 's.message=''No error found in error section.'';\n');
            fprintf(tfid, 'error(s);\n');
        end

        fprintf(tfid, 'end\n');
        fclose(tfid);
        %type(tfile); % DEBUG

        % reload function definition:
        clear(tfile);
        if strcmp(test_type, 'priv_method')
            % we have added a new class methods; reload the class:
            clear(which(param));
        else
            clear(which(tfunc));
        end

        %% run test file
        try
            if strcmp(test_type, 'priv_method')
                eval(sprintf('%s(%s);', tfunc, param));
            else
                %which(tfunc),  type(which(tfunc))  % DEBUG
                feval(tfunc);
                %run([tdir filesep tfunc]);
            end
        catch %#ok<CTCH>
            s = lasterror; %#ok<LERR>
            idx = strcmp({s.stack.name}, mfunc);
            if strcmp(test_type, 'priv_func') && any(idx)
                s.stack(idx).file = mfile;  % it needs full path
            end
            type(which(tfunc));
            s.message = [s.message, get_stack_trace(s)];
            %if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
                delete(tfile);
            %end
            path(old_path);
            if create_pdir,  rmdir(pdir);  end
            rethrow(s);
        end
    end

    %% remove test file
    if exist(tfile, 'file'),  delete(tfile);  end
    path(old_path);
    if create_pdir,  rmdir(pdir);  end
end

%!#test
%! % test()
%! myassert (false)

%!test
%! % test
%! myassert(true);
%! lasterr('', 'test:testing');

%!test
%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:testing');
%! lasterr('', '');

%!error
%! % test
%! myassert(false);

%!shared
%! var = pi;

%!test
%! % test
%! myassert (var == pi);
%! var = 2*pi;

%!test
%! % test
%! myassert (var == pi);

%!shared
%! var = 3*pi;

%!test
%! % test
%! myassert (var == 3*pi);

%!shared  % to clear previous %!shared

%!error
%! test;

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:notEnoughInputs');

%!error
%! test(tempname);

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:UndefinedFunction');

%!error
%! test(1, 2, 3, 4);

%! % test
%! s = lasterror;
%! myassert(s.identifier, 'test:tooManyInputs');

%!test
%! % It's not an error to have a single, empty test section.
%! func = tempname;
%! file = [func '.m'];
%!
%! fid = fopen(file, 'wt');
%! fprintf(fid, '%%!test\n');
%! fprintf(fid, '%%! %% %s\n', func);
%! fclose(fid);
%!
%! old_path = path;
%! addpath(fileparts(file));
%!
%! lasterr ('', '');
%! try,  test(func);  end
%! s = lasterror;
%! %type (file);  % DEBUG
%! delete (file);
%! path(old_path);
%! myassert (isempty(s.identifier));

%!test
%! % Issue an error if there is shared section but not test or error sections.
%! func = tempname;
%! file = [func '.m'];
%!
%! fid = fopen(file, 'wt');
%! fprintf(fid, '%%!shared\n');
%! fprintf(fid, '%%! \n');
%! fclose(fid);
%!
%! old_path = path;
%! addpath(fileparts(file));
%!
%! lasterr ('', '');
%! try,  test(func);  end
%! s = lasterror;
%! %type (file);  % DEBUG
%! delete (file);
%! path(old_path);
%! myassert (s.identifier, 'test:noBlock');

%!shared
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         afunc.m
%! % Then test() creates:
%! % dir/
%! %     atest.m
%!
%! class_code = {...
%!     'function answer = aclass (value)'
%!     '    if (nargin < 1),  value = [];  end'
%!     '    answer.value = value;'
%!     '    answer = class(answer, ''aclass'');'
%!     'end'
%! };
%! func_code = {...
%!     'function answer = afunc (obj)'
%!     '    answer = obj.value;'
%!     'end'
%!     ''
%!     '%!shared'
%!     '%! value = pi;'
%!     '%! obj = aclass(value);'
%!     ''
%!     '%!error'
%!     '%! % afunc'
%!     '%! myassert(false);'
%!     ''
%!     '%!test'
%!     '%! myassert(afunc(obj), value);'
%! };
%! %class_code, func_code  % DEBUG
%!
%! dir = tempname;
%! class_dir = fullfile(dir, '@aclass');
%! temp = mkdir(class_dir);
%!     myassert(temp ~= 0);
%!
%! class_file = fullfile(class_dir, 'aclass.m');
%! [class_fid, msg] = fopen (class_file, 'wt');
%!     if (class_fid == -1),  error(msg);  end
%! fprintf(class_fid, '%s\n', class_code{:});
%! temp = fclose(class_fid);
%!     myassert (temp ~= -1);
%!
%! func_dir = class_dir;
%! func_file = fullfile(func_dir, 'afunc.m');
%! func_fid = fopen (func_file, 'wt');
%!     myassert (func_fid ~= -1);
%! fprintf(func_fid, '%s\n', func_code{:});
%! temp = fclose(func_fid);
%!     myassert (temp ~= -1);
%!
%! %type(class_file);  type(func_file);  % DEBUG
%! %pause  % DEBUG
%!
%! old_path = path;
%! addpath (dir);  % NOT addpath (class_dir);

%!test
%! % Test a class method:
%! % test('afunc', 'aclass');
%!
%! try,
%!     test('afunc', 'aclass');
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!test
%! % Test a class private method:
%! % test('afunc', 'aclass', 'adir');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         private/
%! %             afunc.m
%! % Then test() creates:
%! % dir/
%! %     @aclass/
%! %         atest.m
%!
%! func_file_old = func_file;
%! func_dir = fullfile(class_dir, 'private');
%! temp = mkdir(func_dir);
%!     myassert(temp ~= 0);
%! func_file = fullfile(func_dir, 'afunc.m');
%! [status, msg, msgid] = movefile (func_file_old, func_file);
%!     if (status ~= 1),  error(msgid, msg);  end
%!
%! try,
%!     %keyboard  % DEBUG
%!     test('afunc', 'aclass', func_dir);
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!test
%! % Test a class private method, letting test() guess adir:
%! % test('afunc', 'aclass', 'private');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     @aclass/
%! %         aclass.m  % class constructor
%! %         private/
%! %             afunc.m
%! % Then test() creates:
%! % dir/
%! %     @aclass/
%! %         atest.m
%!
%! func_file_old = func_file;
%! func_dir = fullfile(class_dir, 'private');
%! temp = mkdir(func_dir);
%!     myassert(temp ~= 0);
%! func_file = fullfile(func_dir, 'afunc.m');
%! [status, msg, msgid] = movefile (func_file_old, func_file);
%!     if (status ~= 1),  error(msgid, msg);  end
%!
%! try,
%!     %keyboard  % DEBUG
%!     test('afunc', 'aclass', 'private');
%! catch,
%!     s = lasterror;
%!     if ~strcmp(s.identifier, 'MATLAB:UndefinedFunction')
%!         path (old_path);
%!         [status, msg, msgid] = rmdir(dir, 's');
%!             if (status == 0),  error(msgid, msg);  end
%!     end
%!     rethrow (s);
%! end
%!
%! path (old_path);
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!shared  % to clear previous %!shared

%!test
%! % Test a private function:
%! % test('afunc', [], 'adir');
%! %
%! % Here we generate the following directory-tree:
%! % dir/
%! %     private/
%! %         afunc.m
%! % Then test() creates:
%! % dir/
%! %     atest.m
%!
%! func_code = {...
%!     'function out = afunc (in)'
%!     '    if isequal(in, 0),  error(''afunc:any'', '''');  end'
%!     '    out = 1/in;'
%!     'end'
%!     ''
%!     '%!error'
%!     '%! afunc(0)';
%!     ''
%!     '%!test'
%!     '%! myassert(afunc(pi), 1/pi);'
%! };
%! %func_code  % DEBUG
%!
%! dir = tempname;
%! temp = mkdir(dir);
%!     myassert(temp ~= 0);
%! temp = mkdir(dir, 'private');
%!     myassert(temp ~= 0);
%!
%! func_dir = fullfile(dir, 'private');
%! func_file = fullfile(func_dir, 'afunc.m');
%! [func_fid, msg] = fopen (func_file, 'wt');
%!     if (func_fid == -1),  error(msg);  end
%! fprintf(func_fid, '%s\n', func_code{:});
%! temp = fclose(func_fid);
%!     myassert (temp ~= -1);
%! %type(func_file);  % DEBUG
%!
%! old_path = path;
%! addpath (dir);  % NOT addpath (func_dir);
%! test('afunc', [], func_dir);
%! path (old_path);
%!
%! [status, msg, msgid] = rmdir(dir, 's');
%!     if (status == 0),  error(msgid, msg);  end

%!#test
%! % Check if test files are being deleted,
%! % even in case of error.
%!
%! code = {...
%!     'function answer = func'
%!     '    persistent temp'
%!     '    if isempty(temp)'
%!     '        temp = which(get_caller_name);'
%!     '    end'
%!     '    answer = temp;'
%!     'end'
%! };
%! test1 = {...
%!     '%!test'
%!     '%! func;'
%!     '%! myassert(true);'
%! };
%! test2 = {...
%!     '%!test'
%!     '%! func;'
%!     '%! myassert(false);  % THROW ERROR'
%! };
%!
%! dir = tempname;
%! mkdir(dir);
%! file = fullfile(dir, 'func.m');
%!
%! old_path = path;
%! addpath(dir);
%!
%! for i=1:2
%!     switch i
%!     case 1
%!         testi = test1;
%!     case 2
%!         testi = test2;
%!     end
%!
%!     fid = fopen(file, 'wt');
%!     fprintf(fid, '%s\n', code{:});
%!     fprintf(fid, '%s\n', testi{:});
%!     fclose(fid);
%!
%!     clear('func');
%!     try,  test('func');  end
%!
%!     tfile = func;
%!     try,
%!         switch i
%!         case 1
%!             myassert (~exist(tfile, 'file'));
%!         case 2
%!             myassert (~exist(tfile, 'file'));
%!         end
%!     catch,
%!         s = lasterror;
%!         path(old_path);
%!         rmdir(dir, 's');
%!         rethrow (s);
%!     end
%! end
%! path(old_path);
%! rmdir(dir, 's');

%!test
%! % test ()
%! lastwarn('', '');  % clear last warning
%! warning('off', 'test:noFuncCall');

%!test
%! % generate warning

%!test
%! % test ()
%! [msg, msgid] = lastwarn;
%! myassert (msgid, 'test:noFuncCall');  % check warning
%! warning('on', 'test:noFuncCall');

%!test
%! % test ()
%! lastwarn('', '');  % clear last warning

%!shared
%! % test ()  % avoid test()'s warning

%!test
%! [msg, msgid] = lastwarn;
%! myassert (msgid, '');  % check warning is still clear


text:
function [J] = kJ2J(kJ)
% Convert energy or work from kilojoules to joules.
% Chad A. Greene 2012
J = kJ*1000;


text:
%% R_squared
% Below is a demonstration of the features of the |R_squared| function

%%
clear; close all; clc;

%% Syntax
% |[R_sq]=R_squared(y,yf);|

%% Description
% This function calculates the coefficient of determination (R-Squared) for
% the (e.g. measurement) data |y| and the data |yf| (e.g. a model fit).
% NaN entries in either inputs are ignored.

%% Examples
%

%%
% Plot settings
lineWidth=3;
markerSize=50;
fontSize=25;

%%
% Example data
t=linspace(0,2*pi,25);
y=sin(t)+0.1*randn(size(t));

%%
% Example fit
yf=sin(t);

%%
% Compute R-squared
[R_sq]=R_squared(y,yf)

%%
% Visualize

cFigure; hold on;
title(['R^2=',sprintf('%.4f',R_sq)],'Interpreter','Tex');
plot(t,y,'k.','MarkerSize',markerSize);
plot(t,yf,'g-','LineWidth',lineWidth);
set(gca,'FontSize',fontSize);
axis tight; axis square; box on; grid on;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% GPTOOLBOXES Load in the relevant toolboxes for GP.

importLatest('netlab');
importLatest('mocap');
importTool('ndlutil');
importLatest('prior');
importLatest('mltools');
importLatest('optimi');
importLatest('datasets');
importLatest('kern');
importLatest('noise'); % Only needed for C++ load ins.

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [cls,M1] = spm_preproc_write8(res,tc,bf,df,mrf,cleanup,bb,vx)
% Write out VBM preprocessed data
% FORMAT [cls,M1] = spm_preproc_write8(res,tc,bf,df,mrf,cleanup,bb,vx)
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_preproc_write8.m 6881 2016-09-19 09:48:54Z john $

% Prior adjustment factor.
% This is a fudge factor to weaken the effects of the tissue priors.  The
% idea here is that the bias from the tissue priors probably needs to be
% reduced because of the spatial smoothing typically used in VBM studies.
% Having the optimal bias/variance tradeoff for each voxel is not the same
% as having the optimal tradeoff for weighted averages over several voxels.

if isfield(res,'mg')
    lkp = res.lkp;
    Kb  = max(lkp);
else
    Kb  = size(res.intensity(1).lik,2);
end
N   = numel(res.image);

if nargin<2, tc = true(Kb,4); end % native, import, warped, warped-mod
if nargin<3, bf = false(N,2); end % field, corrected
if nargin<4, df = false(1,2); end % inverse, forward
if nargin<5, mrf = 1;         end % MRF parameter
if nargin<6, cleanup = 1;     end % Run the ad hoc cleanup

% Read essentials from tpm (it will be cleared later)
tpm = res.tpm;
if ~isstruct(tpm) || ~isfield(tpm, 'bg1'),
    tpm = spm_load_priors8(tpm);
end
d1        = size(tpm.dat{1});
d1        = d1(1:3);
M1        = tpm.M;

% Define orientation and field of view of any "normalised" space
% data that may be generated (wc*.nii, mwc*.nii, rc*.nii & y_*.nii).
if nargin>=7 && any(isfinite(bb(:)))
    % If a bounding box is supplied, combine this with the closest
    % bounding box derived from the dimensions and orientations of
    % the tissue priors.
    if nargin<7, bb = NaN(2,3);   end % Default to TPM bounding box
    if nargin<8, vx = NaN;        end % Default to TPM voxel size
    [bb1,vx1] = spm_get_bbox(tpm.V(1), 'old');
    bb(~isfinite(bb)) = bb1(~isfinite(bb));
    if ~isfinite(vx), vx = abs(prod(vx1))^(1/3); end
    bb(1,:) = vx*round(bb(1,:)/vx);
    bb(2,:) = vx*round(bb(2,:)/vx);
    odim    = abs(round((bb(2,1:3)-bb(1,1:3))/vx))+1;

    mm  = [[bb(1,1) bb(1,2) bb(1,3)
            bb(2,1) bb(1,2) bb(1,3)
            bb(1,1) bb(2,2) bb(1,3)
            bb(2,1) bb(2,2) bb(1,3)
            bb(1,1) bb(1,2) bb(2,3)
            bb(2,1) bb(1,2) bb(2,3)
            bb(1,1) bb(2,2) bb(2,3)
            bb(2,1) bb(2,2) bb(2,3)]'; ones(1,8)];
    vx3 = [[1       1       1
            odim(1) 1       1
            1       odim(2) 1
            odim(1) odim(2) 1
            1       1       odim(3)
            odim(1) 1       odim(3)
            1       odim(2) odim(3)
            odim(1) odim(2) odim(3)]'; ones(1,8)];
    mat    = mm/vx3;
else
    % Use the actual dimensions and orientations of
    % the tissue priors.
    odim = tpm.V(1).dim;
    mat  = tpm.V(1).mat;
end

[pth,nam] = fileparts(res.image(1).fname);
ind  = res.image(1).n;
d    = res.image(1).dim(1:3);

[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3 = 1:d(3);

chan(N) = struct('B1',[],'B2',[],'B3',[],'T',[],'Nc',[],'Nf',[],'ind',[]);
for n=1:N
    d3         = [size(res.Tbias{n}) 1];
    chan(n).B3 = spm_dctmtx(d(3),d3(3),x3);
    chan(n).B2 = spm_dctmtx(d(2),d3(2),x2(1,:)');
    chan(n).B1 = spm_dctmtx(d(1),d3(1),x1(:,1));
    chan(n).T  = res.Tbias{n};

    % Need to fix writing of bias fields or bias corrected images, when the data used are 4D.
    [pth1,nam1] = fileparts(res.image(n).fname);
    chan(n).ind = res.image(n).n;

    if bf(n,2)
        chan(n).Nc      = nifti;
        chan(n).Nc.dat  = file_array(fullfile(pth1,['m', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nc.mat  = res.image(n).mat;
        chan(n).Nc.mat0 = res.image(n).mat;
        chan(n).Nc.descrip = 'Bias corrected';
        create(chan(n).Nc);
    end

    if bf(n,1),
        chan(n).Nf      = nifti;
        chan(n).Nf.dat  = file_array(fullfile(pth1,['BiasField_', nam1, '.nii']),...
                                     res.image(n).dim(1:3),...
                                     [spm_type('float32') spm_platform('bigend')],...
                                     0,1,0);
        chan(n).Nf.mat  = res.image(n).mat;
        chan(n).Nf.mat0 = res.image(n).mat;
        chan(n).Nf.descrip = 'Estimated Bias Field';
        create(chan(n).Nf);
    end
end

do_cls   = any(tc(:)) || nargout>=1;
tiss(Kb) = struct('Nt',[]);
for k1=1:Kb
    if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1,
        do_cls  = true;
    end
    if tc(k1,1),
        tiss(k1).Nt      = nifti;
        tiss(k1).Nt.dat  = file_array(fullfile(pth,['c', num2str(k1), nam, '.nii']),...
                                      res.image(n).dim(1:3),...
                                      [spm_type('uint8') spm_platform('bigend')],...
                                      0,1/255,0);
        tiss(k1).Nt.mat  = res.image(n).mat;
        tiss(k1).Nt.mat0 = res.image(n).mat;
        tiss(k1).Nt.descrip = ['Tissue class ' num2str(k1)];
        create(tiss(k1).Nt);
        do_cls = true;
    end
end

prm     = [3 3 3 0 0 0];
Coef    = cell(1,3);
Coef{1} = spm_bsplinc(res.Twarp(:,:,:,1),prm);
Coef{2} = spm_bsplinc(res.Twarp(:,:,:,2),prm);
Coef{3} = spm_bsplinc(res.Twarp(:,:,:,3),prm);

do_defs = any(df);
do_defs = do_cls | do_defs;
if do_defs
    if df(1)
        [pth,nam] =fileparts(res.image(1).fname);
        Ndef      = nifti;
        Ndef.dat  = file_array(fullfile(pth,['iy_', nam, '.nii']),...
                               [res.image(1).dim(1:3),1,3],...
                               [spm_type('float32') spm_platform('bigend')],...
                               0,1,0);
        Ndef.mat  = res.image(1).mat;
        Ndef.mat0 = res.image(1).mat;
        Ndef.descrip = 'Inverse Deformation';
        create(Ndef);
    end
    if df(2) || any(any(tc(:,[2,3,4]))) || nargout>=1,
        y = zeros([res.image(1).dim(1:3),3],'single');
    end
end

spm_progress_bar('init',length(x3),['Working on ' nam],'Planes completed');
M = M1\res.Affine*res.image(1).mat;

if do_cls
    Q = zeros([d(1:3),Kb],'single');
end

for z=1:length(x3)

    % Bias corrected image
    cr = cell(1,N);
    for n=1:N
        f          = spm_sample_vol(res.image(n),x1,x2,o*x3(z),0);
        bf         = exp(transf(chan(n).B1,chan(n).B2,chan(n).B3(z,:),chan(n).T));
        cr{n}      = bf.*f;
        if ~isempty(chan(n).Nc),
            % Write a plane of bias corrected data
            chan(n).Nc.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = cr{n};
        end
        if ~isempty(chan(n).Nf),
            % Write a plane of bias field
            chan(n).Nf.dat(:,:,z,chan(n).ind(1),chan(n).ind(2)) = bf;
        end
    end

    if do_defs
        % Compute the deformation (mapping voxels in image to voxels in TPM)
        [t1,t2,t3] = defs(Coef,z,res.MT,prm,x1,x2,x3,M);

        if exist('Ndef','var')
            % Write out the deformation to file, adjusting it so mapping is
            % to voxels (voxels in image to mm in TPM)
            tmp = M1(1,1)*t1 + M1(1,2)*t2 + M1(1,3)*t3 + M1(1,4);
            Ndef.dat(:,:,z,1,1) = tmp;
            tmp = M1(2,1)*t1 + M1(2,2)*t2 + M1(2,3)*t3 + M1(2,4);
            Ndef.dat(:,:,z,1,2) = tmp;
            tmp = M1(3,1)*t1 + M1(3,2)*t2 + M1(3,3)*t3 + M1(3,4);
            Ndef.dat(:,:,z,1,3) = tmp;
        end

        if exist('y','var')
            % If needed later, save in variable y
            y(:,:,z,1) = t1;
            y(:,:,z,2) = t2;
            y(:,:,z,3) = t3;
        end

        if do_cls
            % Generate variable Q if tissue classes are needed
           %msk = (f==0) | ~isfinite(f);

            if isfield(res,'mg')
                % Parametric representation of intensity distributions
                q   = zeros([d(1:2) Kb]);
                q1  = likelihoods(cr,[],res.mg,res.mn,res.vr);
                q1  = reshape(q1,[d(1:2),numel(res.mg)]);
                b   = spm_sample_priors8(tpm,t1,t2,t3);
                wp  = res.wp;
                s   = zeros(size(b{1}));
                for k1 = 1:Kb
                    b{k1} = wp(k1)*b{k1};
                    s     = s + b{k1};
                end
                for k1=1:Kb
                    q(:,:,k1) = sum(q1(:,:,lkp==k1),3).*(b{k1}./s);
                end
            else
                % Nonparametric representation of intensity distributions
                q   = spm_sample_priors8(tpm,t1,t2,t3);
                wp  = res.wp;
                s   = zeros(size(q{1}));
                for k1 = 1:Kb
                    q{k1} = wp(k1)*q{k1};
                    s     = s + q{k1};
                end
                for k1 = 1:Kb
                    q{k1} = q{k1}./s;
                end
                q   = cat(3,q{:});

                for n=1:N
                    tmp = round(cr{n}*res.intensity(n).interscal(2) + res.intensity(n).interscal(1));
                    tmp = min(max(tmp,1),size(res.intensity(n).lik,1));
                    for k1=1:Kb
                        likelihood = res.intensity(n).lik(:,k1);
                        q(:,:,k1)  = q(:,:,k1).*likelihood(tmp);
                    end
                end
            end
            Q(:,:,z,:) = reshape(q,[d(1:2),1,Kb]);
        end
    end
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');

cls   = cell(1,Kb);
if do_cls
    P = zeros([d(1:3),Kb],'uint8');
    if mrf==0
        % Normalise to sum to 1
        sQ = (sum(Q,4)+eps)/255;
        for k1=1:size(Q,4)
            P(:,:,:,k1) = uint8(round(Q(:,:,:,k1)./sQ));
        end
        clear sQ
    else
        % Use a MRF cleanup procedure
        nmrf_its = 10;
        spm_progress_bar('init',nmrf_its,['MRF: Working on ' nam],'Iterations completed');
        G   = ones([Kb,1],'single')*mrf;
        vx2 = 1./single(sum(res.image(1).mat(1:3,1:3).^2));
        %save PQG P Q G tiss Kb x3 ind
        for iter=1:nmrf_its
            spm_mrf(P,Q,G,vx2);
            spm_progress_bar('set',iter);
        end
    end
    clear Q

    if cleanup
        % Use an ad hoc brain cleanup procedure
        if size(P,4)>5
            P = clean_gwc(P,cleanup);
        else
            warning('Cleanup not done.');
        end
    end

    % Write tissues if necessary
    for k1=1:Kb
        if ~isempty(tiss(k1).Nt)
            for z=1:length(x3)
                tmp = double(P(:,:,z,k1))/255;
                tiss(k1).Nt.dat(:,:,z,ind(1),ind(2)) = tmp;
            end
        end
    end
    spm_progress_bar('clear');

    % Put tissue classes into a cell array...
    for k1=1:Kb
        if tc(k1,4) || any(tc(:,3)) || tc(k1,2) || nargout>=1
            cls{k1} = P(:,:,:,k1);
        end
    end
    clear P % ...and remove the original 4D array
end

clear tpm
M0  = res.image(1).mat;

if any(tc(:,2))
    % "Imported" tissue class images

    % Generate mm coordinates of where deformations map from
    x      = affind(rgrid(d),M0);

    % Generate mm coordinates of where deformation maps to
    y1     = affind(y,M1);

    % Procrustes analysis to compute the closest rigid-body
    % transformation to the deformation, weighted by the
    % interesting tissue classes.
    ind    = find(tc(:,2)); % Saved tissue classes
    [dummy,R]  = spm_get_closest_affine(x,y1,single(cls{ind(1)})/255);
    clear x y1

    mat0   = R\mat; % Voxel-to-world of original image space

    [pth,nam] = fileparts(res.image(1).fname);
    fwhm   = max(vx./sqrt(sum(res.image(1).mat(1:3,1:3).^2))-1,0.01);
    for k1=1:size(tc,1)
        if tc(k1,2)

            % Low pass filtering to reduce aliasing effects in downsampled images,
            % then reslice and write to disk
            tmp1    = decimate(single(cls{k1}),fwhm);
            Ni      = nifti;
            Ni.dat  = file_array(fullfile(pth,['rc', num2str(k1), nam, '.nii']),...
                                 odim,...
                                 [spm_type('float32') spm_platform('bigend')],...
                                 0,1,0);
            Ni.mat         = mat;
            Ni.mat_intent  = 'Aligned';
            Ni.mat0        = mat0;
            Ni.mat0_intent = 'Aligned';
            Ni.descrip     = ['Imported Tissue ' num2str(k1)];
            create(Ni);

            for i=1:odim(3)
                tmp           = spm_slice_vol(tmp1,M0\mat0*spm_matrix([0 0 i]),odim(1:2),[1,NaN])/255;
                Ni.dat(:,:,i) = tmp;
            end
            clear tmp1
        end
    end
end

if any(tc(:,3)) || any(tc(:,4)) || nargout>=1 || df(2)
    % Adjust stuff so that warped data (and deformations) have the
    % desired bounding box and voxel sizes, instead of being the same
    % as those of the tissue probability maps.
    M = mat\M1;
    for i=1:size(y,3)
        t1         = y(:,:,i,1);
        t2         = y(:,:,i,2);
        t3         = y(:,:,i,3);
        y(:,:,i,1) = M(1,1)*t1 + M(1,2)*t2 + M(1,3)*t3 + M(1,4);
        y(:,:,i,2) = M(2,1)*t1 + M(2,2)*t2 + M(2,3)*t3 + M(2,4);
        y(:,:,i,3) = M(3,1)*t1 + M(3,2)*t2 + M(3,3)*t3 + M(3,4);
    end
    M1 = mat;
    d1 = odim;
end

if any(tc(:,3)) || any(tc(:,4)) || nargout>=1

    if any(tc(:,3))
        C = zeros([d1,Kb],'single');
    end

    spm_progress_bar('init',Kb,'Warped Tissue Classes','Classes completed');
    for k1 = 1:Kb
        if ~isempty(cls{k1})
            c = single(cls{k1})/255;
            if any(tc(:,3))
                [c,w]       = spm_diffeo('push',c,y,d1(1:3));
                vx          = sqrt(sum(M1(1:3,1:3).^2));
                spm_field('boundary',1);
                C(:,:,:,k1) = spm_field(w,c,[vx  1e-6 1e-4 0  3 2]);
                clear w
            else
                c      = spm_diffeo('push',c,y,d1(1:3));
            end
            if nargout>=1
                cls{k1} = c;
            end
            if tc(k1,4)
                N      = nifti;
                N.dat  = file_array(fullfile(pth,['mwc', num2str(k1), nam, '.nii']),...
                                    d1,...
                                    [spm_type('float32') spm_platform('bigend')],...
                                    0,1,0);
                N.mat  = M1;
                N.mat0 = M1;
                N.descrip = ['Jac. sc. warped tissue class ' num2str(k1)];
                create(N);
                N.dat(:,:,:) = c*abs(det(M0(1:3,1:3))/det(M1(1:3,1:3)));
            end
            spm_progress_bar('set',k1);
        end
    end
    spm_progress_bar('Clear');

    if any(tc(:,3))
        spm_progress_bar('init',Kb,'Writing Warped Tis Cls','Classes completed');
        C = max(C,eps);
        s = sum(C,4);
        for k1=1:Kb
            if tc(k1,3)
                N      = nifti;
                N.dat  = file_array(fullfile(pth,['wc', num2str(k1), nam, '.nii']),...
                                    d1,'uint8',0,1/255,0);
                N.mat  = M1;
                N.mat0 = M1;
                N.descrip = ['Warped tissue class ' num2str(k1)];
                create(N);
                N.dat(:,:,:) = C(:,:,:,k1)./s;
            end
            spm_progress_bar('set',k1);
        end
        spm_progress_bar('Clear');
        clear C s
    end
end

if df(2)
    y         = spm_diffeo('invdef',y,d1,eye(4),M0);
    y         = spm_extrapolate_def(y,M1);
    N         = nifti;
    N.dat     = file_array(fullfile(pth,['y_', nam, '.nii']),...
                           [d1,1,3],'float32',0,1,0);
    N.mat     = M1;
    N.mat0    = M1;
    N.descrip = 'Deformation';
    create(N);
    N.dat(:,:,:,:,:) = reshape(y,[d1,1,3]);
end

return;

%==========================================================================
% function [x1,y1,z1] = defs(sol,z,MT,prm,x0,y0,z0,M)
%==========================================================================
function [x1,y1,z1] = defs(sol,z,MT,prm,x0,y0,z0,M)
iMT = inv(MT);
x1  = x0*iMT(1,1)+iMT(1,4);
y1  = y0*iMT(2,2)+iMT(2,4);
z1  = (z0(z)*iMT(3,3)+iMT(3,4))*ones(size(x1));
x1a = x0    + spm_bsplins(sol{1},x1,y1,z1,prm);
y1a = y0    + spm_bsplins(sol{2},x1,y1,z1,prm);
z1a = z0(z) + spm_bsplins(sol{3},x1,y1,z1,prm);
x1  = M(1,1)*x1a + M(1,2)*y1a + M(1,3)*z1a + M(1,4);
y1  = M(2,1)*x1a + M(2,2)*y1a + M(2,3)*z1a + M(2,4);
z1  = M(3,1)*x1a + M(3,2)*y1a + M(3,3)*z1a + M(3,4);
return;

%==========================================================================
% function t = transf(B1,B2,B3,T)
%==========================================================================
function t = transf(B1,B2,B3,T)
if ~isempty(T)
    d2 = [size(T) 1];
    t1 = reshape(reshape(T, d2(1)*d2(2),d2(3))*B3', d2(1), d2(2));
    t  = B1*t1*B2';
else
    t = zeros(size(B1,1),size(B2,1),size(B3,1));
end
return;

%==========================================================================
% function p = likelihoods(f,bf,mg,mn,vr)
%==========================================================================
function p = likelihoods(f,bf,mg,mn,vr)
K  = numel(mg);
N  = numel(f);
M  = numel(f{1});
cr = zeros(M,N);
for n=1:N
    if isempty(bf)
        cr(:,n) = double(f{n}(:));
    else
        cr(:,n) = double(f{n}(:).*bf{n}(:));
    end
end
p  = ones(numel(f{1}),K);
for k=1:K
    amp    = mg(k)/sqrt((2*pi)^N * det(vr(:,:,k)));
    d      = bsxfun(@minus,cr,mn(:,k)')/chol(vr(:,:,k));
    p(:,k) = amp*exp(-0.5*sum(d.*d,2)) + eps;
end
return;

%==========================================================================
% function dat = decimate(dat,fwhm)
%==========================================================================
function dat = decimate(dat,fwhm)
% Convolve the volume in memory (fwhm in voxels).
lim = ceil(2*fwhm);
x  = -lim(1):lim(1); x = spm_smoothkern(fwhm(1),x); x  = x/sum(x);
y  = -lim(2):lim(2); y = spm_smoothkern(fwhm(2),y); y  = y/sum(y);
z  = -lim(3):lim(3); z = spm_smoothkern(fwhm(3),z); z  = z/sum(z);
i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;
spm_conv_vol(dat,dat,x,y,z,-[i j k]);
return;

%==========================================================================
% function y1 = affind(y0,M)
%==========================================================================
function y1 = affind(y0,M)
y1 = zeros(size(y0),'single');
for d=1:3
    y1(:,:,:,d) = y0(:,:,:,1)*M(d,1) + y0(:,:,:,2)*M(d,2) + y0(:,:,:,3)*M(d,3) + M(d,4);
end
return;

%==========================================================================
% function x = rgrid(d)
%==========================================================================
function x = rgrid(d)
x = zeros([d(1:3) 3],'single');
[x1,x2] = ndgrid(single(1:d(1)),single(1:d(2)));
for i=1:d(3)
    x(:,:,i,1) = x1;
    x(:,:,i,2) = x2;
    x(:,:,i,3) = single(i);
end
return;

%==========================================================================
% function [P] = clean_gwc(P,level)
%==========================================================================
function [P] = clean_gwc(P,level)
if nargin<2, level = 1; end

b    = P(:,:,:,2);

% Build a 3x3x3 seperable smoothing kernel
%--------------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

th1 = 0.15;
if level==2, th1 = 0.2; end
% Erosions and conditional dilations
%--------------------------------------------------------------------------
niter  = 32;
niter2 = 32;
spm_progress_bar('Init',niter+niter2,'Extracting Brain','Iterations completed');
for j=1:niter
    if j>2, th=th1; else th=0.6; end  % Dilate after two its of erosion
    for i=1:size(b,3)
        gp       = double(P(:,:,i,1));
        wp       = double(P(:,:,i,2));
        bp       = double(b(:,:,i))/255;
        bp       = (bp>th).*(wp+gp);
        b(:,:,i) = uint8(round(bp));
    end
    spm_conv_vol(b,b,kx,ky,kz,-[1 1 1]);
    spm_progress_bar('Set',j);
end

% Also clean up the CSF.
if niter2 > 0
    c = b;
    for j=1:niter2
        for i=1:size(b,3)
            gp       = double(P(:,:,i,1));
            wp       = double(P(:,:,i,2));
            cp       = double(P(:,:,i,3));
            bp       = double(c(:,:,i))/255;
            bp       = (bp>th).*(wp+gp+cp);
            c(:,:,i) = uint8(round(bp));
        end
        spm_conv_vol(c,c,kx,ky,kz,-[1 1 1]);
        spm_progress_bar('Set',j+niter);
    end
end

th = 0.05;
for i=1:size(b,3)
    slices = cell(1,size(P,4));
    for k1=1:size(P,4)
        slices{k1} = double(P(:,:,i,k1))/255;
    end
    bp        = double(b(:,:,i))/255;
    bp        = ((bp>th).*(slices{1}+slices{2}))>th;
    slices{1} = slices{1}.*bp;
    slices{2} = slices{2}.*bp;

    if niter2>0
        cp        = double(c(:,:,i))/255;
        cp        = ((cp>th).*(slices{1}+slices{2}+slices{3}))>th;
        slices{3} = slices{3}.*cp;
    end
    slices{5} = slices{5}+1e-4; % Add a little to the soft tissue class
    tot       = zeros(size(bp))+eps;
    for k1=1:size(P,4)
        tot   = tot + slices{k1};
    end
    for k1=1:size(P,4)
        P(:,:,i,k1) = uint8(round(slices{k1}./tot*255));
    end
end
spm_progress_bar('Clear');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%MAKEFIGS_TRACKEXTRACT  Makes a sample figure for TRACKEXTRACT.

load tpjaos
use tpjaos

region=[-66 -41 52 67];
[lat,lon,ssh,mss,atd]=trackextract(lat,lon,ssh,mss,atd,region);

figure
axis(region)
topoplot continents
plot(lon,lat,'linewidth',2)
latratio(58.5)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [sh,pnl,pos] = marisa(x,N,M)
% this is a combo model, MA+RSI
thresh=55;
S=length(x);
e=ema(x,N);
% take the RSI of the *DETRENDED* series
r=rsi2(x-ema(x,15*M),M);
% bid/ask spreads
cost=0.01;     % BUND/BOBL
%cost=0.005;    % SCHATZ
%cost=1/64;     %Treasury futures
%cost = 0.0001;  %US/EUR
%cost = 0.05;   %EUR/JPY
%cost = 0.04;   %US/JPY
%cost = 0.0005;  %AUS/USD
%cost= 0.02;
%cost=1;        % eurostoxx

rpos=zeros(S,1);
% position of the ema
epos=sign(x-e);
% position of the rsi
I= (thresh-r)<0;
IU=[~1;I(1:end-1) & ~I(2:end)];
rpos(IU)=-1;
% crossing the lower threshold
I= ((100-thresh)-r)>0;
IL=[~1;I(1:end-1) & ~I(2:end)];
rpos(IL)=1;
% copy down previous position values
for i=2:S
    if rpos(i)==0;
        rpos(i)=rpos(i-1);
    end
end

% get the combined signal position
pos=(rpos+epos)/2;
pos( abs(pos) < 1)=0;

% pnl calculation
pnl= pos(1:end-1).*diff(x) - abs(diff(pos))*cost/2;
sh= sqrt(250)*mean(pnl)/std(pnl);


text:
function [x,tspan] = impliciteuler(f,tspan,x0)
%% IMPLICITEULER - Numerical ODE solver: Implicit Euler method
%
% Syntax:
%   [x,tspan] = euler(f,tspan,x0)
%
% In:
%   f      - Function handle, f(x,t)
%   tspan  - Time steps to simulate, [t0,...,tend]
%   x0     - Initial condition
%
% Out:
%   x      - Solved values
%   tspan  - Time steps
%
% Description:
%   Integrates the system of differential equations
%     x' = f(x,t),  for x(0) = x0
%   over the time interval defined in tspan.
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%%

  % Number of steps
  steps = numel(tspan);

  % Allocate space
  x = zeros(size(x0,1),steps);

  % Initial state
  x(:,1) = x0;

  % Iterate
  for k=2:steps

    % Time discretization
    dt = tspan(k)-tspan(k-1);

    % Step solve the algebraic problem
    x(:,k) = fsolve(@(z) x(:,k-1) + f(z,tspan(k))*dt - z, ...
        x(:,k-1), optimset('display','none'));

  end

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
classdef CEC2017_F24 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{12}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{12}.M_10;
            elseif obj.D < 50
                obj.D   = 30;
                obj.Mat = Data{12}.M_30;
            elseif obj.D < 100
                obj.D   = 50;
                obj.Mat = Data{12}.M_50;
            else
                obj.D   = 100;
                obj.Mat = Data{12}.M_100;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z = Z*obj.Mat';
            PopObj = max(abs(Z),[],2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z      = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z      = Z*obj.Mat';
            PopObj = max(abs(Z),[],2);
            PopCon(:,1) = sum(Z.^2,2) - 100*size(Z,2);
            PopCon(:,2) = abs(cos(PopObj)+sin(PopObj)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Fig. 6.82   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clf
clear all;
%close all;

numS=[1 11 10 0];
denS=[1 11 60 100];
w=logspace(-1,2,100);
sysS=tf(numS,denS);
[m,p]=bode(sysS,w);
loglog(w,squeeze(m));
axis([.1 100 .01 10])

hold on
% add line at mag = 1
w2=[.1 100];
mcl=[1 1];
loglog(w2,mcl,'r')

xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.82 Sensitivity Function for Example 6.24');
bodegrid;
hold off


text:
function [msh,lights,tenseMsh] = fs_FSMesh2MLR(fsMeshName,mmPerVox, host, id, varargin);
%
%   [msh,lights,tenseMsh] = fs_fs_FSMesh2MLR(fsMeshName,mmPerVox, host, id, varargin);
%
%
% Author: ARW (based on RFD's mrmBuildMesh)
% Purpose:
%  Take a mesh in Freesurfer sorf format and convert it into a mrLoadRet /
%  mrMesh-type mesh. Allows you to do the usual mrMesh type things like
%  relaxation...
%  Additional processing options can be set as well.  These
%  options include:
%     * 'RelaxIterations'- the next value specifies how many extra smoothing iterations.
%     * 'SavetenseMsh'- if you specify any RelaxIterations and add this option, then
% you will get two meshes- the relaxed mesh and the unrelaxed (tense) mesh.
%
% Returns:
% The mesh structure and a lights structure.
% You can also get the unrelaxed mesh out (tenseMsh).
%
% See Also
%  mrmMapVerticesToGray
%
mrGlobals;

summaryParams = 1;
meshName = '';
QueryFlag = 1;
relaxIter = 1;
saveTense = 0;
backColor = [0.2 0 0];  % CHanged from Gray - partly just to keep track of EMSE meshes...

if (ieNotDefined('fsMeshName'))
   [fsMeshName,fsPath]=uigetfile('*.*','Pick a freesurfer mesh file');
   fsMeshName=fullfile(fsPath,fsMeshName);

end
if ieNotDefined('mmPerVox'),
    mmPerVox=[1 1 1];
    disp('Assuming 1 1 1 mm per vox');
end
if ieNotDefined('host'), host = 'localhost'; end
if ieNotDefined('id'), id = 1; end

if(nargout>2),  saveTense = 1; end

% Parse the varargin values
for(ii=1:length(varargin))
    if    (strcmpi(varargin{ii}, 'RelaxIterations')), relaxIter = varargin{ii+1};
    elseif(strcmpi(varargin{ii}, 'QueryFlag')),       QueryFlag = varargin{ii+1};
    elseif(strcmpi(varargin{ii}, 'MeshName')),        meshName = varargin{ii+1};
    elseif (strcmpi(varargin{ii},'Background')),      backColor = varargin{ii+1};
    end
end

% Set initial parameters for the mesh.
msh = meshDefault(host,id,mmPerVox,relaxIter,meshName);

if QueryFlag,  msh = meshQuery(msh,summaryParams); end

% If the window is already open, no harm is done.
msh = mrmInitHostWindow(msh)

[msh, lights] = mrmInitMesh(msh,backColor);

disp('Importing mesh data from Freesurfer')
%

% ***************************************************************
% We need to construct a basic msh structure from the file data.

[path,name,ext,ver]=fileparts(fsMeshName);

if (strcmp(upper(ext),'.TRI'))
    [vertex,face]=freesurfer_read_tri(fsMeshName);
else

[vertex,face]=freesurfer_read_surf(fsMeshName);
end

% Now  I'm a little curious about this: It seemed to me that the faces in
% regular FS mesh files are in quad format (4 vertices per face). And so I
% don't understand how mrMesh can deal with this. But using mris_convert
% you can generate triangle mesh files from regular freesurfer meshes so I
% propose to do this and avoid all confusion.

size(vertex)
size(face)
vertex=vertex'+128;
face=face';
%face=face([3 2 1],:);

% Vertices come back rotated in an odd way. The entire cortex is rotated 90
% degrees about the L/R axis (i.e. the sag view is rotated 90degrees
% anticlockwise).
vertex=vertex([2 3 1],:);

v=squeeze(vertex([1,2],:));
v=v-128;

rotMat=[-1 0;0 -1];
v=v'*rotMat;
vertex([1,2],:)=v'+128;

nVerts=length(vertex);

p.vertices=(vertex);
p.triangles=face;
p.triangles = p.triangles - 1
p.class='mesh';

% Sometimes we pass in the whole array of meshes.  Mostly, just one,
% though.
host = meshGet(msh,'host');
windowID = meshGet(msh,'windowid');

if isempty(host), host = 'localhost'; end

if isempty(windowID), error('Mesh must specify a window'); end
[id, status, result] = mrMesh ('localhost', windowID, 'add_actor', p)
p.scale = meshGet(msh,'mmPerVox');
p = setSmooth(p,msh,1);
p = setDecimate(p,msh,1);
p.actor = meshGet(msh,'actor');
disp(p);

p.colors=ones(4,length(p.vertices))*255;

[a,b,c]=mrmesh('localhost',id,'set_mesh',p)

        p.scale = meshGet(msh,'mmPerVox');
        p = setSmooth(p,msh,1);
        p = setDecimate(p,msh,1);
        %p.actor = meshGet(msh,'actor');
        p.normals=mrmGet(msh,'normals');
        p.actor = meshGet(msh,'actor');

%disp('Building smoothed and decimated mesh for display...');
%mrmSet(msh,'buildMeshAndDecimateAndSmooth',voxels);
%[msh] = mrmSet(msh,'smooth')

% This is a little 'center object' routine.  We could put this into mrmSet,
% really.

vertices = mrmGet(msh,'vertices');
%keyboard
mrmSet(msh,'origin',-mean(vertices'));

% Save these unsmoothed data.
unSmoothedData = mrmGet(msh,'data');
msh.initVertices=vertices;
msh.grayLayers=3;
view=getSelectedVolume;
    disp('Finding vertex to gray map');
    % You have to do this before smoothing...

v2gMap = mrmMapVerticesToGray(vertices,view.nodes,[1 1 1],view.edges,4);

msh.vertexGrayMap=v2gMap;
msh.grayToVertexMap = mrmMapGrayToVertices(view.nodes,vertices, [1 1 1]);

% Now also try to find the other mapping: the mapping from the mesh to all
% the gray nodes: So that we can ask: for any arbitray gray node, which
% mesh point is it closest to?

% If we smooth, the mesh, it is done here.
if(relaxIter>0)

    mrmSet(msh,'smooth');

    % We get the curvature colors from the uninflated mesh
    % We should figure out the correct value to use as our threshold. The mean
    % isn't ideal, since it is moved around by the large areas with arbitrary
    % curvature (eg. corpus callosum). What we really want is the value that
    % corresponds to zero curvature.
    % maybe make the specific curvature map colors adjustable?
    % We now use the actual curvature values, so we know that zero is zeros
    % curvature.
    disp('Setting up the curvature colors');

    % Attach curvature data to the mesh.  We turn on the color later, I think.
    msh=mrmSet(msh,'curvature');
    p.colors=[repmat(msh.curvature,3,1); ones(1,length(msh.curvature))]*255;
    %[a,b,c]=mrmesh('localhost',id,'set_mesh',p)

    curvColorIntensity = 128*meshGet(msh,'curvatureModDepth'); % mesh.curvature_mod_depth;

    monochrome = uint8(round((double(msh.curvature>0)*2-1)*curvColorIntensity+127.5));

    msh = mrmSet(msh,'colors',monochrome);

    disp('Storing the smoothed mesh data computed by mrMesh...');
    data = mrmGet(msh,'data');
    msh = meshSet(msh,'data',data);
    msh = meshSet(msh,'connectionMatrix',1);

else
    % We don't smooth the data with mrMesh.  We just assign it.
    msh = meshSet(msh,'data',unSmoothedData);
    msh = meshSet(msh,'connectionMatrix',1);
end

% In either case, we return a version of the data without smoothing.  This
% mesh is used to register with the gray coordinates.  This is a little
% disorganized.  If we don't smooth, tenseMesh is identically msh.
tenseMsh = msh;
tenseMsh = meshSet(tenseMsh,'data',unSmoothedData);

return;

%---------------------------------------
function p = setSmooth(p,mesh,val)
if val
    p.do_smooth = 1;
    p.smooth_iterations = meshGet(mesh,'smoothiterations');
    p.smooth_relaxation = meshGet(mesh,'smoothrelaxation');
    p.smooth_sinc_method = meshGet(mesh,'smoothmethod');
    p.do_smooth_pre = meshGet(mesh,'smooth_pre');
else
    p.do_smooth = 0;
    p.do_smooth_pre = 0;
end

return;

%----------------------------
function p = setDecimate(p,mesh,val)

if val
    p.do_decimate = 1;
    p.decimate_reduction = meshGet(mesh,'decimatereduction');
    p.decimate_iterations = meshGet(mesh,'decimateiterations');
else
    p.do_decimate = 0;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function kdtree_compile()
% clc
localpath = fileparts(which('kdtree_compile'));
fprintf(1,'Compiling kdtree library [%s]...\n', localpath);

err = 0;
err = err | mex('-outdir',localpath, [localpath,'/kdtree_build.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_delete.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_k_nearest_neighbors.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_ball_query.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_nearest_neighbor.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_range_query.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_io_from_mat.cpp']);
err = err | mex('-outdir',localpath, [localpath,'/kdtree_io_to_mat.cpp']);

if err ~= 0,
   error('compile failed!');
else
   fprintf(1,'\bDone!\n');
end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function y = modefilt1( x, s )
% One-dimensional mode filtering.
%
% Applies a width s one-dimensional mode filter to vector x. That is each
% element of the output y(j) corresponds to the mode of x(j-r:j+r), where
% r~s/2. At boundary regions, y is calculated on smaller windows, for
% example y(1) is calculated over x(1:1+r).  Note that for this function to
% make sense x should take on only a small number of discrete values
% (running time is actually proportional to number of unique values of x).
% This function is modeled after medfilt1, which is part of Matlab's
% 'Signal Processing Toolbox'.
%
% USAGE
%  y = modefilt1( x, s )
%
% INPUTS
%  x   - length n vector
%  s   - filter size
%
% OUTPUTS
%  y   - filtered vector x
%
% EXAMPLE
%  x=[0 1 0 0 0 3 0 1 3 1 2 2 0 1]; s=3;
%  xmedian = medfilt1( x, s ); % may not be available
%  xmode = modefilt1( x, s );
%  [x; xmedian; xmode]
%
% See also MEDFILT1
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.35
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get unique values in x
[vals,~,inds]=unique(x(:)'); m=length(vals); n=length(x);
if(m>256), warning('x takes on large number of diff vals'); end %#ok<WNTAG>

% create quantized representation [H(i,j)==1 iff x(j)==vals(i)]
H=zeros(m,n); H(sub2ind2([m,n],[inds; 1:n]'))=1;

% create histogram [H(i,j) is count of x(j-r:j+r)==vals(i)]
H=localSum(H,[0 s],'same');

% compute mode for each j and map inds back to original vals
[~,inds]=max(H,[],1); y=vals(inds);

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [s,time] = modulation(x,Ts,Nos,Fc)
% Ts : Sampling period
% Nos: Oversampling factor
% Fc : Carrier frequency
Nx=length(x);  offset = 0;
if nargin<5
    scale = 1;
    T=Ts/Nos; % Scale and Oversampling period for Baseband
else
    scale = sqrt(2);
    T=1/Fc/2/Nos; % Scale and Oversampling period for Passband
end
t_Ts = [0:T:Ts-T];
time = [0:T:Nx*Ts-T]; % One sampling interval and whole interval
tmp = 2*pi*Fc*t_Ts+offset;
len_Ts=length(t_Ts);
cos_wct = cos(tmp)*scale;
sin_wct = sin(tmp)*scale;
%s = zeros(N*len_Ts,1);
for n = 1:Nx
   s((n-1)*len_Ts+1:n*len_Ts) = real(x(n))*cos_wct-imag(x(n))*sin_wct;
end

author: LyricYang


text:
function [mps2] = cmps22mps2(cmps2)
% Convert acceleration from centimeters per square centimeter to meters
% per second-squared.
% Chad A. Greene 2012
mps2 = cmps2*1e-2;
end

author: Sable

repo: mcbench-benchmarks


text:
function [last] = filepartsLast(fname)
% filepartsLast Return the last part of a filename (strip off directory and suffix)
% function filepartsLast(fname)
%
% Examples
% filepartsLast('C:/foo/bar') = 'bar'
% filepartsLast('C:/foo/bar.mat') = 'bar'
% filepartsLast('C:/foo/bar.mat.gz') = 'bar.mat'
% filepartsLast('bar.mat') = 'bar'

[pathstr,name,ext,versn] = fileparts(fname);
last = name;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef TestFastNlMeansDenoisingColored
    %TestFastNlMeansDenoisingColored

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','lena.jpg'), ...
                'Color',true, 'ReduceScale',2);
            out = cv.fastNlMeansDenoisingColored(img, 'H',10, 'HColor',10);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            % we use IMNOISE from Image Processing Toolbox
            if ~mexopencv.require('images')
                error('mexopencv:testskip', 'toolbox');
            end

            img = cv.imread(fullfile(mexopencv.root(),'test','lena.jpg'), ...
                'Color',true, 'ReduceScale',2);
            img = imnoise(img, 'gaussian');
            out = cv.fastNlMeansDenoisingColored(img, 'H',10, 'HColor',10);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.fastNlMeansDenoisingColored();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function out = rdivide( dividend, divisor )
%RDIVIDE for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% determine input types
inType = {class(dividend), class(divisor)};
if(strcmp(inType{1},'TRAFO') && strcmp(inType{2},'TRAFO'))
    [outA, idx] = flattenCellMatrix(dividend.data);
    [outB, idxB] = flattenCellMatrix(divisor.data);
    meta = dividend.meta;
    meta_b = divisor.meta;
    clear 'dividend' 'divisor'

    % idx must be consistent
    if(~isempty(idx) && ~isempty(idxB) && ~all(cellfun(@(x,y) isequal(x,y), idx, idxB)))
        error('TRAFO::rdivide: Nested cells must have the same size');
    end
    if(iscell(meta) && iscell(meta_b))
        if(~isempty(meta) && ~isempty(meta_b))
            if(iscell(meta{1}) && iscell(meta_b{1}))
                if(~all(cellfun(@(x,y) isequal(x,y), flattenCellMatrix(meta), flattenCellMatrix(meta_b))))
                    error('TRAFO::rdivide: Unequal reconstruction information');
                end
            else
                if(~all(cellfun(@(x,y) isequal(x,y), meta, meta_b)))
                    error('TRAFO::rdivide: Unequal reconstruction information');
                end
            end
        end
    else
        if ~isequal(meta, meta_b)
            error('TRAFO::rdivide: Unequal reconstruction information');
        end
    end
    out = cellfun(@(x,y) x./y, outA, outB, 'UniformOutput', false);

elseif(strcmp(inType{1},'TRAFO')) % assume divisor is double/uint/int
    [out, idx] = flattenCellMatrix(dividend.data);
    meta = dividend.meta;
    clear 'dividend'

    out = cellfun(@(x) x./divisor, out, 'UniformOutput', false);

elseif(strcmp(inType{2},'TRAFO'))
    [out, idx] = flattenCellMatrix(divisor.data);
    meta = divisor.meta;
    clear 'divisor'

    out = cellfun(@(x) dividend./x, out, 'UniformOutput', false);

else
    error('TRAFO::rdivide: Impossible constellation!');
end

out = reconFlatCellMatrix(out,idx);
out = TRAFO(out,meta); % return TRAFO object again without modifying input

end

author: thomaskuestner


text:
function value = r8_mach ( i )

%*****************************************************************************80
%
%% R8_MACH returns double precision real machine constants.
%
%  Discussion:
%
%    Assume that double precision real numbers are stored with a mantissa
%    of T digits in base B, with an exponent whose value must lie
%    between EMIN and EMAX.  Then for values of I between 1 and 5,
%    D1MACH will return the following values:
%
%      D1MACH(1) = B^(EMIN-1), the smallest positive magnitude.
%      D1MACH(2) = B^EMAX*(1-B^(-T)), the largest magnitude.
%      D1MACH(3) = B^(-T), the smallest relative spacing.
%      D1MACH(4) = B^(1-T), the largest relative spacing.
%      D1MACH(5) = log10(B)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Phyllis Fox, Andrew Hall, Norman Schryer,
%    Algorithm 528,
%    Framework for a Portable Library,
%    ACM Transactions on Mathematical Software,
%    Volume 4, Number 2, June 1978, page 176-188.
%
%  Parameters:
%
%    Input, integer I, chooses the parameter to be returned.
%    1 <= I <= 5.
%
%    Output, real VALUE, the value of the chosen parameter.
%
  if ( i < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 5.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0.0;
    error ( 'R8_MACH - Fatal error!' );
  elseif ( i == 1 )
    value = 1.112536929253601E-308;
  elseif ( i == 2 )
    value = 4.494232837155789E+307;
  elseif ( i == 3 )
    value = 1.110223024625157E-016;
  elseif ( i == 4 )
    value = 2.220446049250313E-016;
  elseif ( i == 5 )
    value = 0.301029995663981;
  elseif ( 5 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 5.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0.0;
    error ( 'R8_MACH - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Xw, Yw, Zw]=imageToWorld(Xi, Yi, camPar)

% if scalar, just up/downscale = orthographic
if camPar.ortho
    Xw=Xi*camPar.scale;
    Yw=Yi*camPar.scale;
    Zw=0;
else

    mGeo=camPar.mGeo;
    mExt=camPar.mExt;
    mInt=camPar.mInt;

    mTx=mExt.mTx;
    mTy=mExt.mTy;
    mTz=mExt.mTz;

    mT=[mExt.mTx;mExt.mTy;mExt.mTz];

    %% internal init
    sa = sin(mExt.mRx);
    ca = cos(mExt.mRx);
    sb = sin(mExt.mRy);
    cb = cos(mExt.mRy);
    sg = sin(mExt.mRz);
    cg = cos(mExt.mRz);

    mR11 = cb * cg;
    mR12 = cg * sa * sb - ca * sg;
    mR13 = sa * sg + ca * cg * sb;
    mR21 = cb * sg;
    mR22 = sa * sb * sg + ca * cg;
    mR23 = ca * sb * sg - cg * sa;
    mR31 = -sb;
    mR32 = cb * sa;
    mR33 = ca * cb;

    % 		/* convert from image to distorted sensor coordinates */
    Xd = mGeo.mDpx * (Xi - mInt.mCx) / mInt.mSx;
    Yd = mGeo.mDpy * (Yi - mInt.mCy);

    % 		/* convert from distorted sensor to undistorted sensor plane coordinates */
    [Xu Yu]=distortedToUndistortedSensorCoord (Xd, Yd, mInt.mKappa1);

    % 		/* calculate the corresponding xw and yw world coordinates	 */
    % 		/* (these equations were derived by simply inverting	 */
    % 		/* the perspective projection equations using Macsyma)	 */
    Zw=0;
    common_denominator = ((mR11 * mR32 - mR12 * mR31) * Yu + ...
        (mR22 * mR31 - mR21 * mR32) * Xu - ...
        mInt.mFocal * mR11 * mR22 + mInt.mFocal * mR12 * mR21);

    Xw = (((mR12 * mR33 - mR13 * mR32) * Yu + ...
        (mR23 * mR32 - mR22 * mR33) * Xu - ...
        mInt.mFocal * mR12 * mR23 + mInt.mFocal * mR13 * mR22) * Zw + ...
        (mR12 * mTz - mR32 * mTx) * Yu + ...
        (mR32 * mTy - mR22 * mTz) * Xu - ...
        mInt.mFocal * mR12 * mTy + mInt.mFocal * mR22 * mTx) / common_denominator;

    Yw = -(((mR11 * mR33 - mR13 * mR31) * Yu + ...
        (mR23 * mR31 - mR21 * mR33) * Xu - ...
        mInt.mFocal * mR11 * mR23 + mInt.mFocal * mR13 * mR21) * Zw + ...
        (mR11 * mTz - mR31 * mTx) * Yu + ...
        (mR31 * mTy - mR21 * mTz) * Xu - ...
        mInt.mFocal * mR11 * mTy + mInt.mFocal * mR21 * mTx) / common_denominator;
% else
%     error('imageToWorld: camera parameters format unknown');
end
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function s=getStatisticsStr(obj, format)
    format2printer=struct();
    format2printer.text=@obj2stat_str_text;
    format2printer.xml=@obj2stat_str_xml;

    f=format2printer.(format);
    s=f(obj);

function s=obj2stat_str_xml(obj)
    % do not return any stats
    s='';

function status=obj2stat_str_text(obj)
    label_verbosity=2;
    label2count=struct();

    for i=1:countTestOutcomes(obj)
        test_outcome=getTestOutcome(obj,i);

        if isSuccess(test_outcome)
            % no need to report statistics for successes
            continue;
        end

        test_label=getOutcomeStr(test_outcome,label_verbosity);

        if ~isfield(label2count,test_label);
            label2count.(test_label)=0;
        end

        label2count.(test_label)=label2count.(test_label)+1;
    end

    if wasSuccessful(obj)
        status='OK';
    else
        status='FAILED';
    end

    labels=fieldnames(label2count);
    n_labels=numel(labels);

    if n_labels>0
        parts=cell(1,n_labels);
        for k=1:n_labels
            label=labels{k};
            parts{k}=sprintf('%s=%d',label,label2count.(label));
        end

        status=sprintf('%s (%s)',status,moxunit_util_strjoin(parts,', '));
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [ a, seed ] = r8vec_uniform ( n, b, c, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM returns a scaled pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real B, C, the range of the pseudorandom values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM - Fatal error!' );
  end

  r = zeros ( n, 1 );

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    a(i) = b + ( c - b ) * seed * 4.656612875E-10;

  end

  return
end