text:
% reachableForwardRel
%
% A script for computing the forward reachable set, using an relative
% coordinate system. It seems to give roughly reasonable reasonable
% results, but there seems to be some numerical dissapation during the run.
% For example, points that we know are stable appear to be unstable, which
% I think is due to the feature size of the reachable set approaching the
% grid-spacing. This will occur for any grid spacing, since it seems like
% the reachable set becomes an very long and thin ellipsoid-like shape.
%

nReach = 40;  %Number of times at which to display level sets
if (tSpan(2)~=0), error('tSpan(2) must be equal to 0 for reachability calculations!'); end;
% tReach = linspace(tSpan(1),tSpan(2),nReach);
 tReach = linspace(tSpan(1),-1.5,nReach);

radius = 0.3;

nGrid = 400;  %Discritization of the state space in each dimension
grid.dim = 2;  %Two-dimensional problem
grid.min = 8*radius*[-1;-1];
grid.max = 8*radius*[1;1];
grid.dx = (1/(nGrid-1))*(grid.max-grid.min);
grid.bdry = @addGhostExtrapolate;
grid = processGrid(grid);

% Create initial conditions (a circle around the target)
center = [0;0];   %We're working in relative coordinates here!
data = shapeSphere(grid, center, radius); %Initialize level-set function
levelSet = 0;

% Set up the level-set problem:
schemeFunc = @termConvection;
schemeData.velocity = @reachableDynamicsRel;
schemeData.grid = grid;

schemeData.fit.x = xFit;
schemeData.fit.y = yFit;
schemeData.fit.u = uFit;
schemeData.fit.kx = kxFit;
schemeData.fit.ky = kyFit;

% Set up time approximation scheme.
integratorOptions = odeCFLset('factorCFL', 0.5, 'stats', 'on');

% Choose approximations at appropriate level of accuracy.
accuracy = 'medium';
switch(accuracy)
    case 'low'
        schemeData.derivFunc = @upwindFirstFirst;
        integratorFunc = @odeCFL1;
    case 'medium'
        schemeData.derivFunc = @upwindFirstENO2;
        integratorFunc = @odeCFL2;
    case 'high'
        schemeData.derivFunc = @upwindFirstENO3;
        integratorFunc = @odeCFL3;
    case 'veryHigh'
        schemeData.derivFunc = @upwindFirstWENO5;
        integratorFunc = @odeCFL3;
    otherwise
        error('Unknown accuracy level %s', accuracy);
end

%Data structure for storing the reachable set, as contour lines:
Reach(nReach).C = [];

%Store the initial contour line:
Reach(1).C = getContour(grid,data,levelSet);

%Run the backwards reachability calculation:
startTime = cputime;
for idx = 1:(nReach-1)

    %Display data to user:
    figure(12); clf; contour(grid.xs{1},grid.xs{2}, data);

    % Reshape data array into column vector for ode solver call.
    y0 = data(:);

    % How far to step?
    tSpanReach = [tReach(idx), tReach(idx+1)];

    % Take a timestep.
    [t, y] = feval(integratorFunc, schemeFunc, tSpanReach, y0,...
        integratorOptions, schemeData);
    if t~=tSpanReach(2)
        error('Something funny happened!')
    end

    % Get back the correctly shaped data array
    data = reshape(y, grid.shape);

    if min(min(data)) > 0, break;   end

    %Store the level curve:
    Reach(idx+1).C = getContour(grid,data,levelSet);

end

%Plot the results:
figure(7); clf; hold on;
colors = jet(nReach);
xNomReach = polyval(xFit,tReach);
yNomReach = polyval(yFit,tReach);
for i=1:nReach
    C = Reach(i).C;
    if ~isempty(Reach(i).C)
        for j=1:length(C)
            x = C(j).x + xNomReach(i);
            y = C(j).y +yNomReach(i);
            plot(x,y,'LineWidth',3,'color',colors(i,:));
        end
    end
end
plot(xNomReach,yNomReach,'k.','MarkerSize',20);
plot(xSol,ySol,'k')
xlabel('x');
ylabel('y');
zlabel('t');
title('Forward reachable set')

author: MatthewPeterKelly

repo: dscTutorials


text:
function test_old_ft_topoplotTFR

% MEM 2gb
% WALLTIME 00:20:00

% DEPENDENCY_FT_TOPOPLOTTFR
% This script tests the ft_topoplotTFR function and should display a figure
% with the ctf275 layout showing power decreases at the parietal lobes

% load time-frequency data
datadir = dccnpath('/home/common/matlab/fieldtrip/data/test/');
fprintf('loading data\n');
%load observe_comm_moves_freqmtmconvol.mat
load(fullfile(datadir,'observe_comm_moves_freqmtmconvol.mat'));

% topoplot a low frequency band
figure;
cfg              = [];
cfg.baselinetype = 'relchange';
cfg.baseline     = [-2 -.5];
cfg.xlim         = [.5 3];
cfg.ylim         = [8 12];
cfg.zlim         = [-1 1];
ft_topoplotTFR(cfg, obs_lo);

% topoplot a high frequency band
figure;
cfg              = [];
cfg.baselinetype = 'relchange';
cfg.baseline     = [-2 -.5];
cfg.xlim         = [.5 3];
cfg.ylim         = [30 40];
cfg.zlim         = [-1 1];
ft_topoplotTFR(cfg, obs_hi);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear; close all;
x = mandrill;
opt_filters.P = 2;

K = 10;
filters = morlet_filter_bank_2d_pyramid(opt_filters);

opt_wav.J = 5;
tic;
for k = 1:K
	[x_phi, x_psi, meta_phi, meta_psi] = wavelet_2d_pyramid(x, filters, opt_wav);
end
toc;
vx.signal = x_psi;
vx.meta = meta_psi;
ux = modulus_layer(vx);
img = image_scat_layer(ux,0,0);
imagesc(img);

%%
clear;
x = mandrill;
filt_opt.J = 5;
scat_opt.oversampling = 0;
options.L = 8;
K = 10;
[w,filters] = wavelet_factory_2d(size(x), filt_opt, scat_opt);

U{1}.signal{1} = x;
U{1}.meta.j = zeros(0,1);
tic;
for k = 1:K
	[ax,wx] = w{1}(U{1});
end
toc;
ux = modulus_layer(wx);
img = image_scat_layer(ux,0,0);
immac(img,2);
%%

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [c, ceq] = ma_inclinationConstraint(stateLog, eventID, lbInc, ubInc, bodyIDApply, celBodyData, maData)
%ma_semiEccentricityConstraint Summary of this function goes here
%   Detailed explanation goes here
    normFact = pi;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
%         hMAMainGUI = findall(0,'tag','ma_MainGUI');
%         maData = getappdata(hMAMainGUI,'ma_data');
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);

    bodyEventLog = eventLog(eventLog(:,8)==bodyIDApply,:);
    if(isempty(bodyEventLog))
        finalEntry = eventLog(end,:);
    else
        finalEntry = bodyEventLog(end,:);
    end

    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        gmu = bodyInfo.gm;
        rVect = finalEntry(2:4)';
        vVect = finalEntry(5:7)';

        [~, ~, inc, ~, ~, ~] = getKeplerFromState(rVect,vVect,gmu);

        if(lbInc == ubInc)
            c = [0 0];
            ceq(1) = inc - ubInc;
        else
            c(1) = inc - ubInc;
            c(2) = lbInc - inc;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = spm_coreg(varargin)
% Between modality coregistration using information theory
% FORMAT x = spm_coreg(VG,VF,flags)
% VG    - handle for reference image (see spm_vol).
% VF    - handle for source (moved) image.
% flags - a structure containing the following elements:
%          sep      - optimisation sampling steps (mm)
%                     default: [4 2]
%          params   - starting estimates (6 elements)
%                     default: [0 0 0  0 0 0]
%          cost_fun - cost function string:
%                       'mi'  - Mutual Information
%                       'nmi' - Normalised Mutual Information
%                       'ecc' - Entropy Correlation Coefficient
%                       'ncc' - Normalised Cross Correlation
%                     default: 'nmi'
%          tol      - tolerences for accuracy of each param
%                     default: [0.02 0.02 0.02 0.001 0.001 0.001]
%          fwhm     - smoothing to apply to 256x256 joint histogram
%                     default: [7 7]
%          graphics - display coregistration outputs
%                     default: ~spm('CmdLine')
%
% x     - the parameters describing the rigid body rotation, such that a
%         mapping from voxels in G to voxels in F is attained by:
%         VF.mat\spm_matrix(x(:)')*VG.mat
%
% At the end, the voxel-to-voxel affine transformation matrix is
% displayed, along with the histograms for the images in the original
% orientations, and the final orientations. The registered images are
% displayed at the bottom.
%__________________________________________________________________________
%
% The registration method used here is based on the work described in:
% A Collignon, F Maes, D Delaere, D Vandermeulen, P Suetens & G Marchal
% (1995) "Automated Multi-modality Image Registration Based On
% Information Theory". In the proceedings of Information Processing in
% Medical Imaging (1995).  Y. Bizais et al. (eds.).  Kluwer Academic
% Publishers.
%
% The original interpolation method described in this paper has been
% changed in order to give a smoother cost function.  The images are
% also smoothed slightly, as is the histogram.  This is all in order to
% make the cost function as smooth as possible, to give faster convergence
% and less chance of local minima.
%__________________________________________________________________________
% Copyright (C) 1994-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_coreg.m 6435 2015-05-14 09:59:54Z guillaume $

%--------------------------------------------------------------------------
% References
%==========================================================================
%
% Mutual Information
% -------------------------------------------------------------------------
% Collignon, Maes, Delaere, Vandermeulen, Suetens & Marchal (1995).
% "Automated multi-modality image registration based on information theory".
% In Bizais, Barillot & Di Paola, editors, Proc. Information Processing
% in Medical Imaging, pages 263--274, Dordrecht, The Netherlands, 1995.
% Kluwer Academic Publishers.
%
% Wells III, Viola, Atsumi, Nakajima & Kikinis (1996).
% "Multi-modal volume registration by maximisation of mutual information".
% Medical Image Analysis, 1(1):35-51, 1996.
%
% Entropy Correlation Coefficient
% -------------------------------------------------------------------------
% Maes, Collignon, Vandermeulen, Marchal & Suetens (1997).
% "Multimodality image registration by maximisation of mutual
% information". IEEE Transactions on Medical Imaging 16(2):187-198
%
% Normalised Mutual Information
% -------------------------------------------------------------------------
% Studholme, Hill & Hawkes (1998).
% "A normalized entropy measure of 3-D medical image alignment".
% in Proc. Medical Imaging 1998, vol. 3338, San Diego, CA, pp. 132-143.
%
% Optimisation
% -------------------------------------------------------------------------
% Press, Teukolsky, Vetterling & Flannery (1992).
% "Numerical Recipes in C (Second Edition)".
% Published by Cambridge.
%--------------------------------------------------------------------------

SVNid = '$Rev: 6435 $';

if nargin >= 4
    x = optfun(varargin{:});
    return;
end

%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);

def_flags          = spm_get_defaults('coreg.estimate');
def_flags.params   = [0 0 0  0 0 0];
def_flags.graphics = ~spm('CmdLine');
if nargin < 3
    flags = def_flags;
else
    flags = varargin{3};
    fnms  = fieldnames(def_flags);
    for i=1:length(fnms)
        if ~isfield(flags,fnms{i})
            flags.(fnms{i}) = def_flags.(fnms{i});
        end
    end
end

if nargin < 1
    VG = spm_vol(spm_select(1,'image','Select reference image'));
else
    VG = varargin{1};
    if ischar(VG), VG = spm_vol(VG); end
end
if nargin < 2
    VF = spm_vol(spm_select(Inf,'image','Select moved image(s)'));
else
    VF = varargin{2};
    if ischar(VF) || iscellstr(VF), VF = spm_vol(char(VF)); end;
end

if ~isfield(VG, 'uint8')
    VG.uint8 = loaduint8(VG);
    vxg      = sqrt(sum(VG.mat(1:3,1:3).^2));
    fwhmg    = sqrt(max([1 1 1]*flags.sep(end)^2 - vxg.^2, [0 0 0]))./vxg;
    VG       = smooth_uint8(VG,fwhmg); % Note side effects
end

sc = flags.tol(:)'; % Required accuracy
sc = sc(1:length(flags.params));
xi = diag(sc*20);
x = zeros(numel(VF),numel(flags.params));

for k=1:numel(VF)
    VFk = VF(k);
    if ~isfield(VFk, 'uint8')
        VFk.uint8 = loaduint8(VFk);
        vxf       = sqrt(sum(VFk.mat(1:3,1:3).^2));
        fwhmf     = sqrt(max([1 1 1]*flags.sep(end)^2 - vxf.^2, [0 0 0]))./vxf;
        VFk       = smooth_uint8(VFk,fwhmf); % Note side effects
    end

    xk  = flags.params(:);
    for samp=flags.sep(:)'
        xk     = spm_powell(xk(:), xi,sc,mfilename,VG,VFk,samp,flags.cost_fun,flags.fwhm);
        x(k,:) = xk(:)';
    end
    if flags.graphics
        display_results(VG(1),VFk(1),xk(:)',flags);
    end
end

fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#

%==========================================================================
% function o = optfun(x,VG,VF,s,cf,fwhm)
%==========================================================================
function o = optfun(x,VG,VF,s,cf,fwhm)
% The function that is minimised.
if nargin<6, fwhm = [7 7];   end
if nargin<5, cf   = 'mi';    end
if nargin<4, s    = [1 1 1]; end

% Voxel sizes
vxg = sqrt(sum(VG.mat(1:3,1:3).^2));sg = s./vxg;

% Create the joint histogram
H = spm_hist2(VG.uint8,VF.uint8, VF.mat\spm_matrix(x(:)')*VG.mat ,sg);

% Smooth the histogram
lim  = ceil(2*fwhm);
krn1 = spm_smoothkern(fwhm(1),-lim(1):lim(1)) ; krn1 = krn1/sum(krn1); H = conv2(H,krn1);
krn2 = spm_smoothkern(fwhm(2),-lim(2):lim(2))'; krn2 = krn2/sum(krn2); H = conv2(H,krn2);

% Compute cost function from histogram
H  = H+eps;
sh = sum(H(:));
H  = H/sh;
s1 = sum(H,1);
s2 = sum(H,2);

switch lower(cf)
    case 'mi'
        % Mutual Information:
        H   = H.*log2(H./(s2*s1));
        mi  = sum(H(:));
        o   = -mi;
    case 'ecc'
        % Entropy Correlation Coefficient of:
        % Maes, Collignon, Vandermeulen, Marchal & Suetens (1997).
        % "Multimodality image registration by maximisation of mutual
        % information". IEEE Transactions on Medical Imaging 16(2):187-198
        H   = H.*log2(H./(s2*s1));
        mi  = sum(H(:));
        ecc = -2*mi/(sum(s1.*log2(s1))+sum(s2.*log2(s2)));
        o   = -ecc;
    case 'nmi'
        % Normalised Mutual Information of:
        % Studholme,  Hill & Hawkes (1998).
        % "A normalized entropy measure of 3-D medical image alignment".
        % in Proc. Medical Imaging 1998, vol. 3338, San Diego, CA, pp. 132-143.
        nmi = (sum(s1.*log2(s1))+sum(s2.*log2(s2)))/sum(sum(H.*log2(H)));
        o   = -nmi;
    case 'ncc'
        % Normalised Cross Correlation
        i     = 1:size(H,1);
        j     = 1:size(H,2);
        m1    = sum(s2.*i');
        m2    = sum(s1.*j);
        sig1  = sqrt(sum(s2.*(i'-m1).^2));
        sig2  = sqrt(sum(s1.*(j -m2).^2));
        [i,j] = ndgrid(i-m1,j-m2);
        ncc   = sum(sum(H.*i.*j))/(sig1*sig2);
        o     = -ncc;
    otherwise
        error('Invalid cost function specified');
end

%==========================================================================
% function udat = loaduint8(V)
%==========================================================================
function udat = loaduint8(V)
% Load data from file indicated by V into an array of unsigned bytes.
if size(V.pinfo,2)==1 && V.pinfo(1) == 2
    mx = 255*V.pinfo(1) + V.pinfo(2);
    mn = V.pinfo(2);
else
    spm_progress_bar('Init',V.dim(3),...
        ['Computing max/min of ' spm_file(V.fname,'filename')],...
        'Planes complete');
    mx = -Inf; mn =  Inf;
    for p=1:V.dim(3)
        img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
        img = img(isfinite(img));
        mx  = max([max(img(:))+paccuracy(V,p) mx]);
        mn  = min([min(img(:)) mn]);
        spm_progress_bar('Set',p);
    end
end

% Another pass to find a maximum that allows a few hot-spots in the data.
spm_progress_bar('Init',V.dim(3),...
        ['2nd pass max/min of ' spm_file(V.fname,'filename')],...
        'Planes complete');
nh = 2048;
h  = zeros(nh,1);
for p=1:V.dim(3)
    img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
    img = img(isfinite(img));
    img = round((img+((mx-mn)/(nh-1)-mn))*((nh-1)/(mx-mn)));
    h   = h + accumarray(img,1,[nh 1]);
    spm_progress_bar('Set',p);
end
tmp = [find(cumsum(h)/sum(h)>0.9999); nh];
mx  = (mn*nh-mx+tmp(1)*(mx-mn))/(nh-1);

% Load data from file indicated by V into an array of unsigned bytes.
spm_progress_bar('Init',V.dim(3),...
    ['Loading ' spm_file(V.fname,'filename')],...
    'Planes loaded');
udat = zeros(V.dim,'uint8');
st = rand('state'); % st = rng;
rand('state',100); % rng(100,'v5uniform'); % rng('defaults');
for p=1:V.dim(3)
    img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
    acc = paccuracy(V,p);
    if acc==0
        udat(:,:,p) = uint8(max(min(round((img-mn)*(255/(mx-mn))),255),0));
    else
        % Add random numbers before rounding to reduce aliasing artifact
        r = rand(size(img))*acc;
        udat(:,:,p) = uint8(max(min(round((img+r-mn)*(255/(mx-mn))),255),0));
    end
    spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
rand('state',st); % rng(st);

%==========================================================================
% function acc = paccuracy(V,p)
%==========================================================================
function acc = paccuracy(V,p)
if ~spm_type(V.dt(1),'intt')
    acc = 0;
else
    if size(V.pinfo,2)==1
        acc = abs(V.pinfo(1,1));
    else
        acc = abs(V.pinfo(1,p));
    end
end

%==========================================================================
% function V = smooth_uint8(V,fwhm)
%==========================================================================
function V = smooth_uint8(V,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(V.uint8,V.uint8,x,y,z,-[i j k]);

%==========================================================================
% function display_results(VG,VF,x,flags)
%==========================================================================
function display_results(VG,VF,x,flags)
fig = spm_figure('FindWin','Graphics');
if isempty(fig), return; end;
set(0,'CurrentFigure',fig);
spm_figure('Clear','Graphics');

%txt = 'Information Theoretic Coregistration';
switch lower(flags.cost_fun)
    case 'mi',  txt = 'Mutual Information Coregistration';
    case 'ecc', txt = 'Entropy Correlation Coefficient Registration';
    case 'nmi', txt = 'Normalised Mutual Information Coregistration';
    case 'ncc', txt = 'Normalised Cross Correlation';
    otherwise, error('Invalid cost function specified');
end

% Display text
%--------------------------------------------------------------------------
ax = axes('Position',[0.1 0.8 0.8 0.15],'Visible','off','Parent',fig);
text(0.5,0.7, txt,'FontSize',16,...
    'FontWeight','Bold','HorizontalAlignment','center','Parent',ax);

Q = inv(VF.mat\spm_matrix(x(:)')*VG.mat);
text(0,0.5, sprintf('X1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(1,:)),'Parent',ax);
text(0,0.3, sprintf('Y1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(2,:)),'Parent',ax);
text(0,0.1, sprintf('Z1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(3,:)),'Parent',ax);

% Display joint histograms
%--------------------------------------------------------------------------
ax  = axes('Position',[0.1 0.5 0.35 0.3],'Visible','off','Parent',fig);
H   = spm_hist2(VG.uint8,VF.uint8,VF.mat\VG.mat,[1 1 1]);
tmp = log(H+1);
image(tmp*(64/max(tmp(:))),'Parent',ax');
set(ax,'DataAspectRatio',[1 1 1],...
    'PlotBoxAspectRatioMode','auto','XDir','normal','YDir','normal',...
    'XTick',[],'YTick',[]);
title('Original Joint Histogram','Parent',ax);
xlabel(spm_file(VG.fname,'short22'),'Parent',ax,'Interpreter','none');
ylabel(spm_file(VF.fname,'short22'),'Parent',ax,'Interpreter','none');

H   = spm_hist2(VG.uint8,VF.uint8,VF.mat\spm_matrix(x(:)')*VG.mat,[1 1 1]);
ax  = axes('Position',[0.6 0.5 0.35 0.3],'Visible','off','Parent',fig);
tmp = log(H+1);
image(tmp*(64/max(tmp(:))),'Parent',ax');
set(ax,'DataAspectRatio',[1 1 1],...
    'PlotBoxAspectRatioMode','auto','XDir','normal','YDir','normal',...
    'XTick',[],'YTick',[]);
title('Final Joint Histogram','Parent',ax);
xlabel(spm_file(VG.fname,'short22'),'Parent',ax,'Interpreter','none');
ylabel(spm_file(VF.fname,'short22'),'Parent',ax,'Interpreter','none');

% Display ortho-views
%--------------------------------------------------------------------------
spm_orthviews('Reset');
     spm_orthviews('Image',VG,[0.01 0.01 .48 .49]);
h2 = spm_orthviews('Image',VF,[.51 0.01 .48 .49]);
global st
st.vols{h2}.premul = inv(spm_matrix(x(:)'));
spm_orthviews('Space');

spm_print;


text:
function [TorF, vstr, rdate] = have_feature_lu_vec()
%HAVE_FEATURE_LU_VEC  Detect availability/version info for LU vector support
%
%   Feature detection function implementing 'lu_vec' tag for HAVE_FEATURE
%   to detect support for LU(..., 'vector') syntax.
%
%   See also HAVE_FEATURE, LU.

%   MIPS
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MIPS.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mips for more info.

vstr = '';
rdate = '';
if have_feature('matlab') && have_feature('matlab', 'vnum') < 7.003
    TorF = 0;     %% lu(..., 'vector') syntax not supported
else
    TorF = 1;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ alpha, beta, gamma ] = bdf_set ( n )

%*****************************************************************************80
%
%% BDF_SET sets weights for backward differentiation ODE weights.
%
%  Discussion:
%
%    GAMMA * Y(N+1) = Sum ( 1 <= I <= N ) ALPHA(I) * Y(N+1-I)
%                     + dX * BETA * Y'(X(N+1),Y(N+1))
%
%    This is equivalent to the backward differentiation corrector formulas.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%    N must be between 1 and 6.
%
%    Output, real ALPHA(N), BETA, GAMMA, the weights.
%
  alpha = zeros ( n, 1 );

  if ( n == 1 )
    beta =     1.0;
    gamma =    1.0;
    alpha(1) = 1.0;
  elseif ( n == 2 )
    beta =       2.0;
    gamma =      3.0;
    alpha(1) =   4.0;
    alpha(2) = - 1.0;
  elseif ( n == 3 )
    beta =       6.0;
    gamma =     11.0;
    alpha(1) =  18.0;
    alpha(2) = - 9.0;
    alpha(3) =   2.0;
  elseif ( n == 4 )
    beta =       12.0;
    gamma =      25.0;
    alpha(1) =   48.0;
    alpha(2) = - 36.0;
    alpha(3) =   16.0;
    alpha(4) =  - 3.0;
  elseif ( n == 5 )
    beta =        60.0;
    gamma =      137.0;
    alpha(1) =   300.0;
    alpha(2) = - 300.0;
    alpha(3) =   200.0;
    alpha(4) =  - 75.0;
    alpha(5) =    12.0;
  elseif ( n == 6 )
    beta =        60.0;
    gamma =      147.0;
    alpha(1) =   360.0;
    alpha(2) = - 450.0;
    alpha(3) =   400.0;
    alpha(4) = - 225.0;
    alpha(5) =    72.0;
    alpha(6) =  - 10.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BDF_SET - Fatal error!\n' );
    fprintf ( 1, '  Illegal order requested = %d\n', n );
    error ( 'BDF_SET - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = gt( x, y )

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '>' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrspwvt
%TFRSPWVT Unit test for the function TFRSPWV.

%       F. Auger, Dec. 1995 - O. Lemoine, March 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrspwv(sig1);
tfr2=tfrspwv(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrspwv test 1 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrspwv(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrspwv test 2 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrspwv(sig,1:N,N,[1]);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrspwv test 3 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrspwv(sig,1:N,N,[1]);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrspwv test 4 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrspwv(sig,N/2+2,N,tftb_window(11,'rect'),tftb_window(N+1,'rect'));
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>2.0*eps),
 error('tfrspwv test 5 failed');
end;

% A SPWVD with a Dirac time-smoothing window is a PWVD
sig=noisecg(N);
tfr1=tfrspwv(sig,1:N,N,1);
tfr2=tfrpwv(sig);
if any(any(abs(tfr1-tfr2)>sqrt(eps))),
 error('tfrspwv test 6 failed');
end;

N=127;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrspwv(sig1);
tfr2=tfrspwv(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrspwv test 7 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrspwv(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrspwv test 8 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrspwv(sig,1:N,N,[1]);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrspwv test 9 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrspwv(sig,1:N,N,[1]);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrspwv test 10 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrspwv(sig,round(N/2)+2,N,tftb_window(11,'rect'),tftb_window(N,'rect'));
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>sqrt(eps)),
 error('tfrspwv test 11 failed');
end;

% A SPWVD with a Dirac time-smoothing window is a PWVD
sig=noisecg(N);
tfr1=tfrspwv(sig,1:N,N,1);
tfr2=tfrpwv(sig);
if any(any(abs(tfr1-tfr2)>sqrt(eps))),
 error('tfrspwv test 12 failed');
end;

author: HeLiangHIT


text:
function [success] = convertInscopixIsxdToHdf5(inputFilePath,varargin)
	% Converts Inscopix proprietary ISXD to HDF5.
	% By default all dropped frames are included in the output as a frame full of zeros.
	% Biafra Ahanonu
	% started: 2019.01.15 [21:17:45]
	% inputs
		% inputMoviePath - char: path to ISXD file.
	% outputs
		%
	% changelog
		% 2022.01.26 [07:46:01] - Update +image package to +inscopix.

	[success] = ciapkg.inscopix.convertInscopixIsxdToHdf5(inputFilePath,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function B = sub(A, varargin)
% sub(A,i,j,k) = A(i;j;k)
% Written by Mo Chen (sth4nth@gmail.com).
assert(ndims(A)==numel(varargin));
sz = cellfun(@numel,varargin);
IDX = cell(1,ndims(A));
for i = 1:ndims(A)
    idx = varargin{i};
    shape = ones(1,ndims(A));
    shape(i) = sz(i);
    idx = reshape(idx,shape);
    shape = sz;
    shape(i) = 1;
    idx = repmat(idx,shape);
    IDX{i} = idx(:);
end
B = reshape(A(sub2ind(size(A),IDX{:})),sz);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [sym1,sym2] = extractSym(obj)
% extract left and right symmetry from various objects

if isa(obj,'symmetry')
  sym1 = obj;
  sym2 = obj;
elseif isa(obj,'orientation')
  sym1 = obj.SS;
  sym2 = obj.CS;
else
  sym2 = [];
  try
    sym1 = obj.CS;
  catch
    sym1 = [];
  end
end

end

author: mtex-toolbox

repo: mtex

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


text:
% See http://stackoverflow.com/questions/16146599/create-artificial-data-in-matlab
% and http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly

figure;
hold on;
dotsize = 12;
 colormap([1 0 .5;   % magenta
           0 0 .8;   % blue
           0 .6 0;   % dark green
           .3 1 0]); % bright green

subplot(231);
data = twospirals();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Two spirals');

subplot(232);
data = clusterincluster();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Cluster in cluster');

subplot(233);
data = corners();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Corners');

subplot(234);
data = halfkernel();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Half-kernel');

subplot(235);
data = crescentfullmoon();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Crescent & Full Moon');

subplot(236);
data = outlier();
scatter(data(:,1), data(:,2), dotsize, data(:,3)); axis equal;
title('Outlier');

repo: gspbox


text:
function [S, Serr, per, tau, exctn, lambda] = armode(A, C, th)
%ARMODE	Eigendecomposition of AR model.
%
%  [S,Serr,per,tau,exctn]=ARMODE(A,C,th) computes the
%  eigendecomposition of an AR(p) model that has been fitted using
%  ARFIT. The input arguments of ARMODE are output of ARFIT.
%
%  The columns of the output matrix S contain the estimated eigenmodes
%  of the AR model. The output matrix Serr contains margins of error
%  for the components of the estimated eigenmodes S, such that
%  (S +/- Serr) are approximate 95% confidence intervals for the
%  individual components of the eigenmodes.
%
%  The two-row matrices per and tau contain in their first rows the
%  estimated oscillation period per(1,k) and the estimated damping
%  time tau(1,k) of the eigenmode S(:,k). In their second rows, the
%  matrices per and tau contain margins of error for the periods and
%  damping times, such that
%     ( per(1,k) +/- per(2,k) )   and   ( tau(1,k) +/- tau(2,k) )
%  are approximate 95% confidence intervals for the period and damping
%  time of eigenmode S(:,k).
%
%  For a purely relaxatory eigenmode, the period is infinite (Inf).
%  For an oscillatory eigenmode, the periods are finite.
%
%  The excitation of an eigenmode measures its dynamical importance
%  and is returned as a fraction exctn that is normalized such that
%  the sum of the excitations of all eigenmodes equals one.
%
%  See also ARFIT, ARCONF.

%  Modified 13-Oct-00
%  Author: Tapio Schneider
%	   tapio@gps.caltech.edu

  ccoeff   = .95;                       % confidence coefficient
  m 	   = size(C,1);			% dimension of state space
  p 	   = size(A,2) / m; 		% order of model
  if p <= 0
    error('Order must be greater 0.');
  end

  % Assemble coefficient matrix of equivalent AR(1) model
  A1 	   = [A; eye((p-1)*m) zeros((p-1)*m,m)];

  % Eigenvalues and eigenvectors of coefficient matrix of equivalent
  % AR(1) model
  [BigS,d] = eig(A1);  			% columns of BigS are eigenvectors
  lambda   = diag(d);    		% vector containing eigenvalues
  lambda   = lambda(:); 	        % force lambda to be column vector

  % Warning if the estimated model is unstable
  if any(abs(lambda) > 1)
    warning(sprintf(['The estimated AR model is unstable.\n',...
		     '\t Some excitations may be negative.']))
  end

  % Fix phase of eigenvectors such that the real part and the
  % imaginary part of each vector are orthogonal
  BigS     = adjph(BigS);

  % Return only last m components of each eigenvector
  S 	   = BigS((p-1)*m+1:p*m, :);

  % Compute inverse of BigS for later use
  BigS_inv = inv(BigS);

  % Recover the matrix Uinv that appears in the asymptotic covariance
  % matrix of the least squares estimator (Uinv is output of AR)
  if (size(th,2) == m*p+1)
    % The intercept vector has been fitted by AR; in computing
    % confidence intervals for the eigenmodes, this vector is
    % irrelevant. The first row and first column in Uinv,
    % corresponding to elements of the intercept vector, are not
    % needed.
    Uinv   = th(3:size(th,1), 2:size(th,2));

  elseif (size(th,2) == m*p)
    %  No intercept vector has been fitted
    Uinv   = th(2:size(th,1), :);
  else
    error('Input arguments of ARMODE must be output of ARFIT.')
  end
  % Number of degrees of freedom
  dof 	 = th(1,1);
  % Quantile of t distribution for given confidence coefficient and dof
  t      = tquant(dof, .5+ccoeff/2);

  % Asymptotic covariance matrix of estimator of coefficient matrix A
  Sigma_A  = kron(Uinv, C);

  % Noise covariance matrix of system of relaxators and oscillators
  CovDcpld = BigS_inv(:, 1:m) * C * BigS_inv(:, 1:m)';

  % For each eigenmode j: compute the period per, the damping time
  % tau, and the excitation exctn; also get the margins of error for
  % per and tau
  for j=1:m*p				% eigenmode number
    a		= real(lambda(j)); 	% real part of eigenvalue j
    b 		= imag(lambda(j)); 	% imaginary part of eigenvalue j
    abs_lambda_sq= abs(lambda(j))^2;  	% squared absolute value of eigenvalue j
    tau(1,j) 	= -2/log(abs_lambda_sq);% damping time of eigenmode j

    % Excitation of eigenmode j
    exctn(j) 	= real(CovDcpld(j,j) / (1-abs_lambda_sq));

    % Assemble derivative of eigenvalue with respect to parameters in
    % the coefficient matrix A
    dot_lam 	= zeros(m^2*p, 1);
    for k=1:m
      dot_lam(k:m:k+(m*p-1)*m) = BigS_inv(j,k) .* BigS(1:m*p,j);
    end
    dot_a 	= real(dot_lam); 	% derivative of real part of lambda(j)
    dot_b 	= imag(dot_lam); 	% derivative of imag part of lambda(j)

    % Derivative of the damping time tau w.r.t. parameters in A
    phi 	= tau(1,j)^2 / abs_lambda_sq * (a*dot_a + b*dot_b);
    % Margin of error for damping time tau
    tau(2,j) 	= t * sqrt(phi'*Sigma_A*phi);

    % Period of eigenmode j and margin of error for period. (The
    % if-statements avoid warning messages that may otherwise result
    % from a division by zero)
    if (b == 0 & a >= 0)      % purely real, nonnegative eigenvalue
      per(1,j)	= Inf;
      per(2,j)  = 0;
    elseif (b == 0 & a < 0)   % purely real, negative eigenvalue
      per(1,j)	= 2;
      per(2,j)  = 0;
    else                      % complex eigenvalue
      per(1,j)	= 2*pi/abs(atan2(b,a));

      % Derivative of period with respect to parameters in A
      phi 	= per(1,j)^2 / (2*pi*abs_lambda_sq)*(b*dot_a-a*dot_b);
      % Margin of error for period
      per(2,j) 	= t * sqrt(phi'*Sigma_A*phi);
    end
  end

  % Give the excitation as `relative importance' that sums to one
  exctn 	= exctn/sum(exctn);

  % Compute confidence intervals for eigenmodes
  % -------------------------------------------
  % Shorthands for matrix products
  XX 	        = real(BigS)'*real(BigS);
  YY 	        = imag(BigS)'*imag(BigS);
  XY 	        = real(BigS)'*imag(BigS);

  % Need confidence intervals only for last m rows of BigS
  row1 	        = (p-1)*m+1; % first row for which confidence interval is needed

  mp = m*p;                           	% dimension of equivalent AR(1) model
  for k=1:mp        			% loop over columns of S
    for l=row1:mp  			% loop over rows of S

      % evaluate gradient of S_{lk}
      for ii=1:m
	for jj=1:mp
	  % compute derivative with respect to A(ii,jj)
	  zsum = 0;
	  zkkr = 0; 			% real part of Z_{kk}
	  zkki = 0; 			% imaginary part of Z_{kk}
	  for j=1:mp
	    if (j ~= k) 		% sum up elements that appear in Z_{kk}
	      zjk  = BigS_inv(j,ii)*BigS(jj,k)/(lambda(k)-lambda(j));
	      zjkr = real(zjk);
	      zjki = imag(zjk);
	      zkkr = zkkr+zjki*(XY(k,j)-XY(j,k))-zjkr*(XX(k,j)+YY(k,j));
	      zkki = zkki+zjki*(YY(k,j)-XX(k,j))-zjkr*(XY(k,j)+XY(j,k));
	      zsum = zsum+BigS(l,j)*zjk;
	    end
	  end
	  % now add Z_{kk}
	  zkki = zkki / (XX(k,k)-YY(k,k));
	  grad_S((jj-1)*m+ii) = zsum+BigS(l,k)*(zkkr+i*zkki);
	end
      end
      Serr(l,k) = t * ( sqrt( real(grad_S)*Sigma_A*real(grad_S)') ...
			 + i*sqrt( imag(grad_S)*Sigma_A*imag(grad_S)') );
    end
  end

  % Only return last m*p rows of Serr
  Serr = Serr(row1:m*p, :);

author: Sable

repo: mcbench-benchmarks


text:
function ex2 (n)
%EX2: create an n-by-n 2D mesh, four different ways

% Example:
%   ex2
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 1)
    n = 30 ;
end

subplot (1,2,1) ;

% method 1: create an n-by-n 2D mesh for the 2nd difference operator
tic
ii = zeros (5*n^2, 1) ;
jj = zeros (5*n^2, 1) ;
xx = zeros (5*n^2, 1) ;
k = 1 ;
for j = 0:n-1
    for i = 0:n-1
        s = j*n+i + 1 ;
        ii (k:k+4) = [(j-1)*n+i j*n+(i-1) j*n+i j*n+(i+1) (j+1)*n+i ] + 1 ;
        jj (k:k+4) = [s s s s s] ;
        xx (k:k+4) = [-1 -1 4 -1 -1] ;
        k = k + 5 ;
    end
end

% remove entries beyond the boundary
keep = find (ii >= 1 & ii <= n^2 & jj >= 1 & jj <= n^2) ;
ii = ii (keep) ;
jj = jj (keep) ;
xx = xx (keep) ;
A = sparse (ii,jj,xx) ;
t1 = toc ; disp (t1) ;
% subplot (2,2,1) ;
spy (A)
title (sprintf ('%d-by-%d 2D mesh\n', n, n)) ;

% method 2, using no for loops
tic
nn = 1:n^2 ;
i2 = [nn-n ; nn-1 ; nn ; nn+1 ; nn+n] ;
j2 = repmat (nn, 5, 1) ;
x2 = repmat ([-1 -1 4 -1 -1]', 1, n^2) ;
keep = find (i2 >= 1 & i2 <= n^2 & j2 >= 1 & j2 <= n^2) ;
i2 = i2 (keep) ;
j2 = j2 (keep) ;
x2 = x2 (keep) ;
C = sparse (i2,j2,x2) ;
t2 = toc ; disp (t2) ;

% subplot (2,2,2) ; plot (j2) ;
% title ('2D fast j2') ;
disp (A-C) ;

any (ii-i2)
any (jj-jj)

% method 3: create an n-by-n-by-n 3D mesh for the 2nd difference operator
tic
[A, keep, ii, jj, xx] = mesh3d1 (n) ;
ii = ii (keep) ;
jj = jj (keep) ;
xx = xx (keep) ;
t3 = toc ; disp (t3) ;
tic
E = sparse (ii,jj,xx) ;
t3b = toc ; disp (t3b) ;
subplot (1,2,2) ; spy (E) ;
title (sprintf ('%d-by-%d-by-%d 3D mesh\n', n, n, n)) ;

% method 4, using no for loops
tic
nn = 1:n^3 ;
i2 = [nn-n^2 ; nn-n ; nn-1 ; nn ; nn+1 ; nn+n ; nn+n^2] ;
j2 = repmat (nn, 7, 1) ;
x2 = repmat ([-1 -1 -1 6 -1 -1 -1]', 1, n^3) ;
keep = find (i2 >= 1 & i2 <= n^3 & j2 >= 1 & j2 <= n^3) ;
i2 = i2 (keep) ;
j2 = j2 (keep) ;
x2 = x2 (keep) ;
t4 = toc ; disp (t4) ;
tic
F = sparse (i2,j2,x2) ;
t4b = toc ; disp (t4b) ;
disp (E-F) ;


text:
% batch_computeEntropy.m

dirPath = '/Volumes/deasylab1/Data/NSCLC Study';

dirPath = '/Users/aptea/Documents/MSKCC/Projects/Entropy_test';

%structureName = 'totallung';
structureName = 'ptv';

%Find all CERR files
fileC = {};
if strcmpi(dirPath,'\') || strcmpi(dirPath,'/')
    filesTmp = getCERRfiles(dirPath(1:end-1));
else
    filesTmp = getCERRfiles(dirPath);
end
fileC = [fileC filesTmp];

filesNotConvertedC = {};
etrpy = {};

%Loop over CERR plans
for iFile = 1:length(fileC)

    drawnow

    global planC stateS

    try
        planC = loadPlanC(fileC{iFile},tempdir);
        planC = updatePlanFields(planC);
        indexS = planC{end};
    catch
        continue
    end

    % Combine Iliacs
    structure_index = find(strcmpi(structureName,{planC{indexS.structures}.structureName}));

    if length(structure_index) > 1
        structure_index = structure_index(1);
    end

    if isempty(structure_index)

        filesNotConvertedC = [filesNotConvertedC, fileC{iFile}];
        continue;

    end

    doseNum = 1;

    % Check Units
    if any(strcmpi(planC{indexS.dose}(doseNum).doseUnits,{'cgy','cgrays','cgray'}))
        planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray / 100;
    end

    % Get scan associated with this structure
    scanNum = getStructureAssociatedScan(structure_index, planC);

    % Get structure mask
    [rV, cV, sV] = getUniformStr(structure_index, planC);

    % Get uniformized coords
    [xV, yV, zV] = getUniformScanXYZVals(planC{indexS.scan});

    % Get dose at structure locations
    dosesV = getDoseAt(doseNum, xV(cV), yV(rV), zV(sV), planC);

    % Get DVH
    %[planC, doseBinsV, volsHistV] = getDVHMatrix(planC, structure_index, doseNum);

    % Calculate Entropy
    [jnk,fName] = fileparts(fileC{iFile});
    etrpy{iFile,1} = fName;
    % etrpy{iFile,2} = entropy(dosesV./max(dosesV));
    % Manual Entrop calculation
    dose_hist = hist(dosesV,256);
    dose_hist(dose_hist == 0) = [];
    dose_hist = dose_hist./numel(dosesV);
    etrpy{iFile,2} = -sum(dose_hist.*log2(dose_hist));

    clear global planC stateS

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% The COBRA Toolbox: testFluxSplits.m
%
% Purpose:
%     - tests the predictFluxSplits function
%
% Authors:
%     - Original file: Hulda Haraldsdottir
%     - CI integration: Laurent Heirendt January 2017
%
% Note:
%     - This test only runs with solvers that can solve LP and QP problems

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testFluxSplits'));
cd(fileDir);

% define a toy model with single internal loop
model.mets = {'A'; 'B'; 'C'};
model.rxns = {'R1'; 'R2'; 'R3'; 'U'; 'S'};
model.S = [-1  0 -1 -1  0;
            1 -1  0  0  0;
            0  1  1  0 -1];
model.b = [0; 0; 0];
model.lb = [-1000; -1000; -1000; -10; 0];
model.ub = [1000; 1000; 1000; 0; 10];
model.c = [0; 0; 0; 0; 1];

% define inputs
obj = 'S';
met2test = {'C'};
samples = {'model'};
ResultsAllCellLines.model.modelPruned = model;

% define parameteres for tests
tol = 1e-6;
v_ref = [10/3; 10/3; 20/3; -10; 10]; % reference flux distribution

% list of solver packages
solverPkgs = {'tomlab_cplex', 'gurobi'};

for k = 1:length(solverPkgs)

    fprintf(' -- Running testFluxSplits using the solver interface: %s ... ', solverPkgs{k});

    s1 = changeCobraSolver(solverPkgs{k}, 'LP', 0);
    s2 = changeCobraSolver(solverPkgs{k}, 'QP', 0);

    if s1 == 1 && s2 == 1
        % Test of production
        p.rxn = {'R3'}; % max contributing reaction
        p.flux = [20/3, 10, 200/3]; % contributing fluxes
        [BMall, ResultsAllCellLines, ~, maximum_contributing_rxn, maximum_contributing_flux, ~] = predictFluxSplits(model, obj, met2test, samples, ResultsAllCellLines, 1);

        assert(norm(BMall - v_ref) < tol);
        assert(strcmp(maximum_contributing_rxn,p.rxn));
        assert(norm(maximum_contributing_flux - p.flux) < tol);

        % Test of consumption
        c.rxn = {'S'}; % max contributing reaction
        c.flux = [10, 10, 100]; % contributing fluxes
        [BMall, ResultsAllCellLines, ~, maximum_contributing_rxn, maximum_contributing_flux, ~] = predictFluxSplits(model, obj, met2test, samples, ResultsAllCellLines, 0);

        assert(norm(BMall - v_ref) < tol);
        assert(strcmp(maximum_contributing_rxn,c.rxn));
        assert(norm(maximum_contributing_flux - c.flux) < tol);

        % set a status message
        fprintf('Done.\n');
    else
        warning('The test testFluxSplits cannot run using the solver interface: %s. The solver interface is not installed or not configured properly.\n', solverPkgs{k});
    end
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [z,mz,onrim,Hd_record] = cg_steihaug(grad,hess,Delta,epsilon,maxCG,y,Hd_back,quiet)
% Helper function for trustregion_newton_cg

savemem = true;

grad = grad(:);
onrim = false;
if exist('Hd_back','var') && ~isempty(Hd_back)
    Hd_record = Hd_back;
    backtrack = true;
else
    backtrack = false;
    if savemem
        Hd_record = zeros(length(grad),maxCG,'single');
    else
        Hd_record = zeros(length(grad),maxCG);
    end
end

z = zeros(size(grad));  % z is step, we start at origin
mz = 0;                 % 2nd order prediction for objective change at step z
r = grad;               % 2nd order residual at z:  r = H*z-grad
d = -r;                 % we're going down

residual = sqrt(r'*r)/epsilon;
if residual <= 1
    Hd_record = [];
    fprintf('CG 0: as far as I can see, this should never happen\n');
    fprintf('CG 0: converged with zero step, residual = %g\n',residual);
    return
end

j=0;
while true
    if backtrack && j+1 <= size(Hd_back,2)
        Hd = double(Hd_back(:,j+1));
    else
        Hd = hess(d);
        if j+1 <= size(Hd_record)
            Hd_record(:,j+1) = Hd;
        end
    end
    dHd = d'*Hd;
    if dHd <=0 %region is non-convex in direction d
        a = d'*d;
        b = 2*z'*d;
        c = z'*z-Delta^2;
        discr = sqrt(b^2-4*a*c);
        tau1 = (-b - discr)/(2*a);
        tau2 = (-b + discr)/(2*a);
        model = @(tau) tau*grad'*d + 0.5*tau^2*dHd;
        if model(tau1) < model(tau2)
            tau = tau1;
        else
            tau = tau2;
        end
        z = z + tau*d;
        mz = mz + tau*d'*grad + 0.5*tau^2*dHd;
        onrim = true;
        if ~quiet, fprintf('CG %i: curv=%g, jump to trust region boundary\n',j,dHd); end
        break;
    end

    alpha = r'*r/dHd;
    old_z = z;
    old_mz = mz;
    z = z + alpha*d;
    mz = mz + alpha*d'*grad + 0.5*alpha^2*dHd;

    radius = sqrt(z'*z)/Delta;
    if radius > 1
        a = d'*d;
        b = 2*z'*d;
        c = old_z'*old_z-Delta^2;
        discr = sqrt(b^2-4*a*c);
        tau = (-b + discr)/(2*a);
        z = old_z + tau*d;
        mz = old_mz + tau*d'*grad + 0.5*tau^2*dHd;
        onrim = true;
        if ~quiet, fprintf('CG %i: curv=%g, terminate on trust region boundary, model=%g\n',j,dHd,-mz); end
        break;
    end

    old_r = r;
    r = r + alpha*Hd;
    residual = sqrt(r'*r)/epsilon;
    if residual <= 1
        if ~quiet,
            fprintf('CG %i: curv=%G, converged inside trust region; radius = %g, residual=%g, model=%g\n',j,dHd,radius,residual,-mz);
        end
        break;
    end

    overshot = ~isempty(y) && y+mz <0;
    if overshot
        if ~quiet,
            fprintf('CG %i: curv=%G, overshot inside trust region; radius = %g, residual=%g, model=%g\n',j,dHd,radius,residual,-mz);
        end
        break;
    end

    %stopped = backtrack && j+1 >= maxCG;
    stopped = j+1 >= maxCG;
    if stopped
        if ~quiet,
            fprintf('CG %i: curv=%G, stopped inside trust region; radius = %g, residual=%g, model=%g\n',j,dHd,radius,residual,-mz);
        end
        break;
    end

    beta = (r'*r)/(old_r'*old_r);
    d = -r + beta*d;
    if ~quiet, fprintf('CG %i: curv=%g, radius = %g, residual=%g, model=%g\n',j,dHd,radius,residual,-mz); end
    j = j+1;
end

if j+1<size(Hd_record,2)
    Hd_record = Hd_record(:,1:j+1);
end
return

author: nesl


text:
function res = ball(size, varargin)
%BALL Generate a ball in a matrix in 2 or 3 dimensions
%
% parameters
%  size: radius of the ball, center exclude
%  dim:  dimensions of the ball (2 or 3)
%
% author : Gaetan Lehmann
%

%	HISTORY
%	13/03/2007: update doc (DL)

if ~isempty(varargin)
    dim = varargin{1};
else
    dim = 2;
end

if dim == 2
    if length(size) == 1
        xSize = size;
        ySize = size;
    else
        xSize = size(1);
        ySize = size(2);
    end

    realXSize = xSize*2+1;
    realYSize = ySize*2+1;

    res = zeros([realXSize, realYSize]);
    center = [xSize ySize] + 1;
    radius2 = ([xSize ySize] + 0.5) .^ 2;

    for x = 1:realXSize
        for y = 1:realYSize
            d = [x y] - center;
            d2 = d .^ 2;
            if sum(d2./radius2) <= 1
                res(x, y) = 1;
            end
        end
    end
else
    if length(size) == 1
        xSize = size;
        ySize = size;
        zSize = size;
    else
        xSize = size(1);
        ySize = size(2);
        zSize = size(3);
    end

    realXSize = xSize*2+1;
    realYSize = ySize*2+1;
    realZSize = zSize*2+1;

    res = zeros([realXSize, realYSize, realZSize]);
    center = [xSize ySize zSize] + 1;
    radius2 = ([xSize ySize zSize] + 0.5) .^ 2;

    for x = 1:realXSize
        for y = 1:realYSize
            for z = 1:realZSize
                d = [x y z] - center;
                d2 = d .^ 2;
                if sum(d2./radius2) <= 1
                    res(x, y, z) = 1;
                end
            end
        end
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [oa ua pa K confu]=confusion(true_label,estim_label)
%
% function confusion(true_label,estim_label)
%
% This function compute the confusion matrix and extract the OA, AA
% and the Kappa coefficient.
%
% INPUT
% easy! just read

l=length(true_label);
nb_c=max(true_label);

confu=zeros(nb_c,nb_c);

for i=1:l
  confu(estim_label(i),true_label(i))= confu(estim_label(i),true_label(i))+1;
end

oa=trace(confu)/sum(confu(:)); %overall accuracy
ua=diag(confu)./sum(confu,2);  %class accuracy
pa=diag(confu)./sum(confu,1)';

Po=oa;
Pe=(sum(confu)*sum(confu,2))/(sum(confu(:))^2);

K=(Po-Pe)/(1-Pe);%kappa coefficient

%http://kappa.chez-alice.fr/kappa_intro.htm

author: BehnoodRasti


text:
function value = local_sum(input, window_size)
% LOCAL_SUM Fast computation of 2D local sum using integral images.
%
% Algorithm is described in Viola and Jones, "Robust real-time face
% detection", International Journal of Computer Vision, 2004.  Computation
% time is independent of window size. Precision may be an issue for very
% large images, but a few thousand by a few thousand images seems to work
% with no issues.
%
% Suprisingly, for small window sizes, MATLAB's conv2 actually seems to be
% a faster way to do this, but this implementation gets much faster as
% window sizes get larger.  MATLAB Central's fastrunmean demonstrates
% another way of computing fast means/sums.  Fastrunmean is not quite as
% fast as this integral image implementation, but doesn't have the same
% potential precision issues for very large images.
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

padsize = floor(window_size/2);
padded = padarray(input, padsize+1, 0, 'pre'); % Zeropad to handle local sums at edges
padded = padarray(padded, padsize, 0, 'post'); % Zeropad to handle local sums at edges
padded = cumsum(cumsum(double(padded),1),2); % Create integral image
value = padded(1:size(input,1),1:size(input,2))-... % Compute sum over local window
    padded(1:size(input,1),(window_size(2)+1):end)-...
    padded((window_size(1)+1):end,1:size(input,2))+...
    padded((window_size(1)+1):end,(window_size(2)+1):end);

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%% IMAGE_SCRIPT is a script to call IMAGE_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'IMAGE_SCRIPT\n' );
  fprintf ( 1, '  Filter an image using 3 workers\n' );

  [ x, y, z ] = image_fun ( );


text:
%LOADPCD Load a point cloud from a PCD format file
%
% P = LOADPCD(FNAME) is a set of points loaded from the PCD format
% file FNAME.
%
% For an unorganized point cloud the columns of P represent the 3D points,
% and the rows are: x, y, z, r, g, b, a depending on the FIELDS in the file.
%
% For an organized point cloud P is a 3-dimensional matrix (HxWxN) where the
% N planes are: x, y, z, r, g, b, a depending on the FIELDS in the file.  This
% format is useful since the planes z, r, g, b, a can be considered as images.
%
% Notes::
% - Only the x y z field format are currently supported
% - The file can be in ascii or binary format, binary_compressed is not
%   supported
%
% See also pclviewer, lspcd, loadpcd.
%
% Copyright (C) 2013, by Peter I. Corke

% TODO
% - handle binary_compressed

function points = loadpcd(fname)

    verbose = true;

    fp = fopen(fname, 'r');

    while true
        line = fgetl(fp);

        if line(1) == '#'
            continue;
        end

        [field,remain] = strtok(line, ' \t');
        remain = strtrim(remain);

        switch field
            case 'VERSION'
                continue;
            case 'FIELDS'
                FIELDS = remain;
            case 'TYPE'
                TYPE = remain;
            case 'SIZE'
                sizes = str2num(remain);
            case 'WIDTH'
                width = str2num(remain);
            case 'HEIGHT'
                height = str2num(remain);
            case 'POINTS'
                npoints = str2num(remain);
            case 'COUNT'
                count = str2num(remain);
            case 'DATA'
                mode = remain;
                break;
            otherwise
                warning('unknown field %s\n', field);
        end
    end

    % parse out details of the fields
    %  numFields    the number of fields
    %  sizes        vector of field sizes (in bytes)
    %  types        vector of type identifiers (I U F)
    %  fields       vector of field names (x y z rgb rgba etc)
    numFields = numel(sizes);
    types = cell2mat(regexp(TYPE,'\s+','split'));
    fields = regexp(FIELDS,'\s+','split');

    if verbose
        % the doco says height > 1 means organized, but some old files have
        % height = 1 and width > 1
        if height > 1 && width > 1
            organized = true;
            org = 'organized';
        else
            organized = false;
            org = 'unorganized';
        end
        fprintf('%s: %s, %s, <%s> %dx%d\n', ...
            fname, mode, org, FIELDS, width, height);
        fprintf('  %s; %s\n', TYPE, num2str(sizes));
    end

    if any(count > 1)
        error('can only handle 1 element per dimension');
    end

    switch mode
        case 'ascii'
            format = '';
            for j=1:numFields
                switch types(j)
                    case 'I', typ = 'd';
                    case 'U', typ = 'u';
                    case 'F', typ = 'f';
                end
                format = [format '%' typ num2str(sizes(j)*8)];
            end
            c = textscan(fp, format, npoints);
            points = [];
            for j=1:length(c)
                points = [points; c{j}'];
            end
            if size(points,2) ~= npoints
                error('incorrect number of points in file: was %d, should be %d', ...
                    size(points,2), npoints);
            end

        case 'binary'
            format = '';

            if true || all(types == types(1)) && all(sizes == sizes(1))
                % simple case where all fields have the same length and type

                % map IUF -> int, uint, float
                switch types(1)
                    case 'I'
                        fmt = 'int';
                    case 'U'
                        fmt = 'uint';
                    case 'F'
                        fmt = 'float';
                end

                format = [format '*' fmt num2str(sizes(1)*8)];
                points = fread(fp, [numFields npoints], format);

            else

                % more complex case where fields have different length and type
                % code contributed by Will

                startPos_fp = ftell(fp);

                % Just initialize the xyz portion for now
                points = zeros(3, npoints);

                for i=1:numFields
                    % map each field sequentially, using fread() skip functionality,
                    % essentially interleaved reading of the file

                    % map IUF -> int, uint, float
                    switch types(i)
                        case 'I'
                            fmt = 'int';
                        case 'U'
                            fmt = 'uint';
                        case 'F'
                            fmt = 'float';
                    end

                    format = ['*' fmt num2str(sizes(i)*8)];
                    fseek(fp, startPos_fp + sum(sizes(1:i-1)), 'bof');
                    data = fread(fp, [1 npoints], format, sum(sizes)-sizes(i));

                    switch fields{i}
                        case 'x'
                            points(1,:) = data;
                        case 'y'
                            points(2,:) = data;
                        case 'z'
                            points(3,:) = data;
                        case {'rgb', 'rgba'}
                            points(4,:) = data;
                    end
                end

            end

        case 'binary_compressed'
            % binary part of the file contains:
            %  compressed size of data (uint32)
            %  uncompressed size of data (uint32)
            %  compressed data
            %  junk
            compressed_size = fread(fp, 1, 'uint32');
            uncompressed_size = fread(fp, 1, 'uint32');
            compressed_data = fread(fp, compressed_size, 'uint8')';

            uncompressed_data = lzfd(compressed_data);
            if length(uncompressed_data) ~= uncompressed_size
                error('decompression error');
            end

            % the data is stored unpacked, that is one field for all points,
            % then the next field for all points, etc.

            start = 1;
            for i=1:numFields
                len = sizes(i)*npoints;
                switch types(1)
                    case 'I'
                        fmt = 'int32';
                    case 'U'
                        fmt = 'uint32';
                    case 'F'
                        fmt = 'single';
                end
                field = typecast(uncompressed_data(start:start+len-1), fmt);
                start = start + len;

                switch fields{i}
                    case 'x'
                        points(1,:) = field;
                    case 'y'
                        points(2,:) = field;
                    case 'z'
                        points(3,:) = field;
                    case {'rgb', 'rgba'}
                        points(4,:) = field;
                end
            end

        otherwise
            error('unknown DATA mode: %s', mode);
    end

    if size(points,1) > 3
        % convert RGB from float to rgb
        rgb = typecast(points(4,:), 'uint32');
        switch FIELDS
            case 'x y z rgb'
                R = double(bitand(255, bitshift(rgb, 16))) /255;
                G = double(bitand(255, bitshift(rgb, 8))) /255;
                B = double(bitand(255, rgb)) /255;
                points = [points(1:3,:); R; G; B];

            case 'x y z rgba'
                R = double(bitand(255, bitshift(rgb, 24))) /255;
                G = double(bitand(255, bitshift(rgb, 16))) /255;
                B = double(bitand(255, bitshift(rgb, 8))) /255;
                A = double(bitand(255, rgb)) /255;
                points = [points(1:3,:); R; G; B; A];
        end
    end

    if organized
        % data is an organized point cloud, rearrange it into planes

        points = permute( reshape( shiftdim(points, 1), width, height, []), [2 1 3]);
    end

    fclose(fp);

author: petercorke

repo: machinevision-toolbox-matlab


text:
%SAD	Stereo matching using Sum of Squared Differences metric
%
%	[DISP SCORES] = SSD(LEFT, RIGHT, LEFT_RIGHT, WIN_SIZE, DISP_RANGE)
%
%       DISP = Disparity map
%       SCORES = SSD scores at each disparity
%       LEFT, RIGHT = images
%       LEFT_RIGHT = 'l' if wrt left image, 'r' if wrt right image
%       WIN_SIZE = window size or [x_size y_size]
%       DISP_RANGE = [START_DISPARITY END_DISPARITY]
% Author: Jasmine E. Banks (jbanks@ieee.org)

% Copyright in this software is owned by CSIRO.  CSIRO grants permission to
% any individual or institution to use, copy, modify, and distribute this
% software, provided that:
%
% (a)     this copyright and permission notice appears in its entirety in or
% on (as the case may be) all copies of the software and supporting
% documentation;
%
% (b)     the authors of papers that describe software systems using this
% software package acknowledge such use by citing the paper as follows:
%
%     "Quantitative Evaluation of Matching Methods and Validity Measures for
%     Stereo Vision" by J. Banks and P. Corke, Int. J. Robotics Research,
%     Vol 20(7), 2001; and
%
% (c)     users of this software acknowledge and agree that:
%
%   (i) CSIRO makes no representations about the suitability of this software
%   for any purpose;
%
%   (ii) that the software is provided "as is" without express or implied
%   warranty; and
%
%   (iii) users of this software use the software entirely at their own risk.


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 Hand, Omega=(0,20)x(0,25), level=7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             MI
%   - regularizer          mbElastic
%   - optimizer            lBFGS
% ===============================================================================

% setup data and initialize image viewer
setup2DhandData;
level = 7; omega = ML{level}.omega; m = ML{level}.m;

% initialize the interpolation scheme and coefficients
imgModel('reset','imgModel','splineInter');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

% initialize distance measure
distance('set','distance','MImex','nT',8,'nR',8);

% initialize regularization, note: yc-yRef is regularized, elastic is staggered
regularizer('reset','regularizer','mbElastic','alpha',1e-2,'mu',1,'lambda',0);
y0   = getStaggeredGrid(omega,m); yRef = y0; yStop = y0;

% setup and initialize plots
FAIRplots('set','mode','NPIR-Gauss-Newton','omega',omega,'m',m,'fig',1,'plots',1);
FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

% build objective function, note: T coefficients of template, Rc sampled reference
fctn = @(yc) NPIRBFGSobjFctn(T,Rc,omega,m,yRef,yc); fctn([]); % report status

% -- solve the optimization problem -------------------------------------------
[yc,his] = lBFGS(fctn,y0,'maxIter',500,'Plots',@FAIRplots,'yStop',yStop);
% report results
iter      = size(his.his,1)-2;
reduction = 100*(fctn(y0)-fctn(yc))/abs(fctn(y0));
fprintf('reduction = %s%% after %d iterations\n',num2str(reduction),iter);

%==============================================================================


text:
function [ ux, uy, uz ] = LKPR3D( image1, image2, r, numlevels, iterations, sigma )
%This function estimates deformation between two subsequent 3-D images using
%Lucas-Kanade optical flow equation with pyramidal approach.
%
%   Description :
%
%   -image1, image2 : two subsequent images or frames.
%   -r : redius of neighbourhood, default value is 2.
%   -numlevels : the number of levels in pyramid, default value is 2.
%   -iterations : number of iterations in refinement, default value is 1.
%   -sigma : standard deviation of Gaussian function, default value is 0.5.
%
%   Reference:
%   Lucas, B. D., Kanade, T., 1981. An iterative image registration
%   technique with an application to stereo vision. In: Proceedings of the
%   7th international joint conference on Artificial intelligence - Volume 2.
%   Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, pp. 674-679.
%
%   Author: Mohammad Mustafa
%   By courtesy of The University of Nottingham and Mirada Medical Limited,
%   Oxford, UK
%
%   Published under a Creative Commons Attribution-Non-Commercial-Share Alike
%   3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
%   June 2012

% Default parameters
if nargin==5
    sigma=0.5;
elseif nargin==4
    iterations=1; sigma=0.5;
elseif nargin==3
    numlevels=2; iterations=1; sigma=0.5;
elseif nargin==2
    r=2; numlevels=2; iterations=1; sigma=0.5;
end

currentImage1=image1; currentImage2=image2;

% s contains sizes of each pyramid levels
s=zeros(1,3,numlevels);
s(:,:,1)=(size(image1));

% Each of pyramid levels will have resized image
pyrIm1=zeros(size(image1,1),size(image1,2),size(image1,3),numlevels);
pyrIm2=pyrIm1;
pyrIm1(:,:,:,1)=currentImage1;
pyrIm2(:,:,:,1)=currentImage2;

% Building pyramid by downsampling
for i=2:numlevels
    currentImage1 = sampling3D(currentImage1,1);
    currentImage2 = sampling3D(currentImage2,1);
    % Adjusting the size
    pyrIm1(1:size(currentImage1,1),1:size(currentImage1,2),...
            1:size(currentImage1,3),i)=currentImage1;
    pyrIm2(1:size(currentImage2,1),1:size(currentImage2,2),...
            1:size(currentImage2,3),i)=currentImage2;
    s(:,:,i)=size(currentImage1);
end

% Base operation
currentImage1=pyrIm1(1:s(1,1,numlevels),1:s(1,2,numlevels),...
                        1:s(1,3,numlevels),numlevels);
currentImage2=pyrIm2(1:s(1,1,numlevels),1:s(1,2,numlevels),...
                        1:s(1,3,numlevels),numlevels);

[ux,uy,uz]=LKW3D(currentImage1,currentImage2,r,sigma);

% Refining flow vectors
if iterations>0
    for i=1:iterations
    [ux,uy,uz]=refinedLK3D(ux,uy,uz,currentImage1,currentImage2,r,sigma);
    end
end

% Operations at higher levels of pyramids

for i=(numlevels-1):-1:1
    % Size and magnitudes of flow vectors are upsampled
    temp=2 * sampling3D(ux,2);
    uxInitial=temp(1:s(1,1,i),1:s(1,2,i),1:s(1,3,i));
    temp=2 * sampling3D(uy,2);
    uyInitial=temp(1:s(1,1,i),1:s(1,2,i),1:s(1,3,i));
    temp=2 * sampling3D(uz,2);
    uzInitial=temp(1:s(1,1,i),1:s(1,2,i),1:s(1,3,i));

    currentImage1=pyrIm1(1:s(1,1,i),1:s(1,2,i),1:s(1,3,i),i);
    currentImage2=pyrIm2(1:s(1,1,i),1:s(1,2,i),1:s(1,3,i),i);

    [ux,uy,uz]=refinedLK3D(uxInitial,uyInitial,uzInitial,...
                            currentImage1,currentImage2,r,sigma);
    if iterations>0
        for j=1:iterations
            [ux,uy,uz]=refinedLK3D(ux,uy,uz,currentImage1,currentImage2,2,sigma);
        end
    end
end

end


text:
function obj = setParamValue(obj, varargin)
    % set the actual value of the system parameters
    %
    % Parameters:
    % varargin: the parameter values @type struct

    % convert to struct
    params = struct(varargin{:});

    update_fields = intersect(fieldnames(obj.params_),fieldnames(params));
    nfields = numel(update_fields);
    if nfields > 0
        for i=1:nfields
            field = update_fields{i};
            obj.params_.(field) = params.(field);
        end
    end

end

author: ayonga

repo: frost-dev

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


text:
fprintf(1,'Saving the stereo calibration results in Calib_Results_stereo.mat...\n');

string_save = 'save Calib_Results_stereo om R T recompute_intrinsic_right recompute_intrinsic_left calib_name_left format_image_left type_numbering_left image_numbers_left N_slots_left calib_name_right format_image_right type_numbering_right image_numbers_right N_slots_right fc_left cc_left kc_left alpha_c_left KK_left fc_right cc_right kc_right alpha_c_right KK_right active_images dX dY nx ny n_ima active_images_right active_images_left inconsistent_images center_optim_left est_alpha_left est_dist_left est_fc_left est_aspect_ratio_left center_optim_right est_alpha_right est_dist_right est_fc_right est_aspect_ratio_right history param param_error sigma_x om_error T_error fc_left_error cc_left_error kc_left_error alpha_c_left_error fc_right_error cc_right_error kc_right_error alpha_c_right_error';

for kk = 1:n_ima,
    if active_images(kk),
        string_save = [string_save ' X_left_' num2str(kk)  ' omc_left_' num2str(kk) ' Tc_left_' num2str(kk) ' omc_left_error_' num2str(kk) ' Tc_left_error_' num2str(kk)  ' n_sq_x_' num2str(kk) ' n_sq_y_' num2str(kk)];
    end;
end;
eval(string_save);

repo: ekfmonoslam


text:
function [Cnorm, Curv, Mu, Sigma] = dtiComputeFiberRadiusofCurvatureDistribution(fg,displayDist)

% Calculate the mean curvature of individual fiber first, and then make a
% summary across fibers.
%
% [Cnorm, Curv, Mu, Sigma]= dtiComputeFiberRadiusofCurvatureDistribution(fg,displayDist)
%
% INPUTS:
% fibers   = A fiber group.
% displayDist = 1: Display the distribution plot. 0: Not display.
%
% Outputs:
% Cnorm    = Normalized radius of curvature of each fiber, in units of z-score.
% Curvature      = The mean radius of curvature of each fiber.
%
% Written by Hiromasa Takemura (c) Stanford University 2014

if notDefined('fg'), error('Fiber group required'); end

% Compute the curvature of fibers in individual nodes.
[fibercurvature] = dtiComputeFiberCurvature(fg);

% Calculate the mean of curvature in individual fibers.

for i = 1:length(fibercurvature)
   Curv(i) = 1./mean(fibercurvature{i});

end

% Z-score the curvature
% Here we take the log first so that the distribution of Curv is 'more'
% gaussian and because it it not clipped at 0.
[Cnorm, Mu, Sigma] = zscore(log10(Curv));

% Show a histagram
if displayDist==1
    mrvNewGraphWin('Distibutions of fiber lengths');
    hold on;
    [y, x] = hist(Curv,round(length(fg.fibers)*0.1));
    bar(x,y,'FaceColor','k','EdgeColor','k');
    axis([min(x) max(x) 0 max(y)]);
    xlabel('Curvature');
    plot(10.^[Mu Mu],[0 max(y)],'r','linewidth',2);
    plot(10.^[Mu-Sigma   Mu-Sigma],[0 max(y)],'--r');
    plot(10.^[Mu+Sigma   Mu+Sigma],[0 max(y)],'--r');
    plot(10.^[Mu-2*Sigma Mu-2*Sigma],[0 max(y)],'--r');
    plot(10.^[Mu+2*Sigma Mu+2*Sigma],[0 max(y)],'--r');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function textscatter(varargin)
% textscatter(cx,cy,markersize,markerstyle,textcolor, textoffsetx, textoffsety)
% Uses scatter plot and each marker is labeled by the number of the row in
% the [cx, cy] matrix. Text is offste by [textoffsetx, textoffsety] from
% the markers.
% Examples:
% textscatter(100*rand(1,30), 100*rand(1,30))
% textscatter(100*rand(1,30), 100*rand(1,30),[],'d')
% textscatter(100*rand(1,30), 100*rand(1,30),100,'gx','r', 2,3)

args = varargin;
if nargin<5 % textcolor
%     args{5} = [0 0 0]; % black
    args{5} = [1 0 0]; % red
end
if nargin<6 % textoffset
    args{6}= .5;
    args{7}= .5;
end

cx = args{1};
cy = args{2};
textcolor = args{5};
textoffsetx = args{6};
textoffsety = args{7};

scatter(args{1:min(nargin,4)})
text(cx+textoffsetx, cy+textoffsety, num2str([1:length(cx)]' ),'color' , textcolor)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [R] = quat2mat(q)
  % QUAT2MAT convert quaternions to 3d rotation matrices
  %
  % [R] = quat2mat(q)
  %
  % Input:
  %   q is an m by 4 list of normalized quaternions, [1 i j k]
  % Output:
  %   R is a 3 by 3 by m list of rotation matrices
  %
  % See:
  %   http://en.wikipedia.org/wiki/
  %   Conversion_between_quaternions_and_Euler_angles#Rotation_matrices
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %

  assert(size(q,2) == 4);

  R = zeros(3,3,size(q,1));
  R(1,1,:) = q(:,1).^2 + q(:,2).^2 - q(:,3).^2 - q(:,4).^2;
  R(1,2,:) = 2*(q(:,2).*q(:,3) - q(:,1).*q(:,4));
  R(1,3,:) = 2*(q(:,1).*q(:,3) + q(:,2).*q(:,4));
  R(2,1,:) = 2*(q(:,2).*q(:,3) + q(:,1).*q(:,4));
  R(2,2,:) = q(:,1).^2 - q(:,2).^2 + q(:,3).^2 - q(:,4).^2;
  R(2,3,:) = 2*(q(:,3).*q(:,4) - q(:,1).*q(:,2));
  R(3,1,:) = 2*(q(:,2).*q(:,4) - q(:,1).*q(:,3));
  R(3,2,:) = 2*(q(:,1).*q(:,2) + q(:,3).*q(:,4));
  R(3,3,:) = q(:,1).^2 - q(:,2).^2 - q(:,3).^2 + q(:,4).^2;

  %yy2 = 2.0 .* q(:,2) .* q(:,2);
  %xy2 = 2.0 .* q(:,1) .* q(:,2);
  %xz2 = 2.0 .* q(:,1) .* q(:,3);
  %yz2 = 2.0 .* q(:,2) .* q(:,3);
  %zz2 = 2.0 .* q(:,3) .* q(:,3);
  %wz2 = 2.0 .* q(:,4) .* q(:,3);
  %wy2 = 2.0 .* q(:,4) .* q(:,2);
  %wx2 = 2.0 .* q(:,4) .* q(:,1);
  %xx2 = 2.0 .* q(:,1) .* q(:,1);
  %R(1,1,:) = - yy2 - zz2 + 1;
  %R(1,2,:) = xy2 + wz2;
  %R(1,3,:) = xz2 - wy2;
  %R(2,1,:) = xy2 - wz2;
  %R(2,2,:) = - xx2 - zz2 + 1;
  %R(2,3,:) = yz2 + wx2;
  %R(2,1,:) = xz2 + wy2;
  %R(2,2,:) = yz2 - wx2;
  %R(2,3,:) = - xx2 - yy2 + 1;
end

author: alecjacobson

repo: gptoolbox


text:
function VLMinstall(varargin)
%
% Downloads a set of folders from the LabelMe dataset.
%
% 1) Install all the videos (it can take a long time to complete)
%  LMinstall(HOMEFRAMES, HOMEANNOTATIONS);
%
% 2) install only videos in the specified folder
%  LMinstall(folder, HOMEFRAMES, HOMEANNOTATIONS);
%
% HOMEFRAMES and HOMEANNOTATIONS point to your local destination folders.
%
% This function is useful if you do not want to download the entire
% dataset. You can first browse the dataset, and when you have decided for
% a set of folders that you want to use, you can use this function to
% download only those folders.
%
% Contribute to the dataset by labeling few images:
% http://labelme.csail.mit.edu/
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Video LabelMe is a WEB-based image annotation tool and a Matlab toolbox that allows
% researchers to label images and share the annotations with the rest of the community.
%    Copyright (C) 2009  MIT, Computer Science and Artificial
%    Intelligence Laboratory. Jenny Yuen, Antonio Torralba, Bryan Russell
%
%    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/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If this function stops working, you might need to download the newest
% version of the toolbox in order to update the web addresses:

webpageimg = 'http://labelme.csail.mit.edu/LabelMeVideo/VLMFrames/';
webpageanno = 'http://labelme.csail.mit.edu/LabelMeVideo/VLMAnnotations/';

Narguments = length(varargin);

switch Narguments
    case {0,1}
        error('Not enough input arguments.')
    case 2
        disp('Install full database (more than 10Gb required)')
        HOMEIMAGES = varargin{1};
        HOMEANNOTATIONS = varargin{2};

        disp('generating folder list...')
        list_fr = folderlist(webpageimg);
        list_an = folderlist(webpageanno);
     case 3
         list = varargin{1};
         HOMEIMAGES = varargin{2};
         HOMEANNOTATIONS = varargin{3};

         disp('generating folder list...')
         list_fr = folderlist(webpageimg, list);
         list_an = folderlist(webpageanno, list);
         %     case 4
%         folders = varargin{1};
%         filelist = varargin{2};
%         list = unique(folders);
%         HOMEIMAGES = varargin{3};
%         HOMEANNOTATIONS = varargin{4};
end

Nfolders = length(list_fr);

% create folders:
disp('Create folders...');
for i = 1:length(list_fr)
    mkdir(HOMEIMAGES, list_fr{i});
end
for i = 1:length(list_an)
    mkdir(HOMEANNOTATIONS, list_an{i});
end

disp('Download annotations...')
for f = 1:length(list_an)
    disp(sprintf('Downloading folder %s (%d/%d)...',  list_an{f}, f, Nfolders))
    wpa = [webpageanno '/' list_an{f}];

    annotations = urldir(wpa, 'txt');
    if ~isempty(annotations)
        annotations = {annotations(:).name};
    end

    Nanno = length(annotations);
    for i = 1:Nanno
        disp(sprintf('    Downloading annotation %s (%d/%d)...',  annotations{i}, i, Nanno))
        [F,STATUS] = urlwrite([wpa '/' annotations{i}], fullfile(HOMEANNOTATIONS,list_an{f},annotations{i}));
        if STATUS == 0
            disp(sprintf('annotation file %s does not exist', annotations{i}))
        end
    end
end

disp('Download frames...')
for f = 1:Nfolders
    disp(sprintf('Downloading folder %s (%d/%d)...',  list_fr{f}, f, Nfolders))

    wpi = [webpageimg  '/' list_fr{f}];

    images = urldir(wpi, 'img');
    if ~isempty(images)
        images = {images(:).name};
    end

    Nimages = length(images);
    for i = 1:Nimages
        disp(sprintf('    Downloading image %s (%d/%d)...',  images{i}, i, Nimages))
        [F,STATUS] = urlwrite([wpi '/' images{i}], fullfile(HOMEIMAGES,list_fr{f},images{i}));
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [post,nlZ,dnlZ] = infExact_fastrobust(kmax, hyp, mean, cov, lik, x, y, s)

% Fast, robust inference for a GP with Gaussian likelihood. Compute a
% parametrization of the posterior, the negative log marginal likelihood
% and its derivatives w.r.t. the hyperparameters. See also "help infMethods".
%
% To be used in combination with 'fast' covariance functions (e.g.
% covSEard_fast and such).
%
% This implementation also handles badly conditioned covariance matrices,
% which sometimes become non-positive definite and cause a failure of the
% Cholesky decomposition. In these cases, the covariance matrix is corrected
% to the nearest symmetric positive semidefinite (SPD) matrix in Frobenius
% norm (see [1]; based on an implementation by John d'Errico).
%
% [1] Higham NJ, "Computing a nearest symmetric positive semidefinite
% matrix", Linear Algebra Appl, 1988.
% http://www.sciencedirect.com/science/article/pii/0024379588902236
%
% Added support for 'fast' inference and robust Cholesky decomposition by
% Luigi Acerbi, 2016-01-03.
% Original code by Carl Edward Rasmussen and Hannes Nickisch, 2015-07-13.
%
% See also INFMETHODS.M.

if iscell(lik), likstr = lik{1}; else likstr = lik; end
if ~ischar(likstr), likstr = func2str(likstr); end
if strcmp(likstr,'likGauss')    % NOTE: no explicit call to likGauss
    henoise_flag = false;
elseif strcmp(likstr,'likGaussHe')
    henoise_flag = true;        % Heteroskedastic noise
else
    error('Exact inference only possible with Gaussian likelihood');
end

if isempty(kmax); kmax = 5; end

[n, D] = size(x);
if nargout > 2                                      % do we want derivatives?
  [K,dK] = feval(cov{:}, hyp.cov, x);                 % evaluate covariance matrix and derivatives
else
  K = feval(cov{:}, hyp.cov, x);                      % evaluate covariance matrix
end
m = feval(mean{:}, hyp.mean, x);                          % evaluate mean vector

sn2_base = exp(2*hyp.lik);                               % noise variance of likGauss

if henoise_flag && ~isempty(s)
    % if isempty(s); error('Input-dependent vector S is empty.'); end
    sn2 = sn2_base + s.^2;                               % Vector of observation variance
else
    sn2 = sn2_base;
end

Lchol = min(sn2) >= 1e-6;       % tiny sn2 can lead to numerical trouble

% smallnoise = sn2 < 1e-6;
if Lchol
    if isscalar(sn2)
        sn2div = sn2;
        sn2_mat = eye(n);
    else
        sn2div = min(sn2);
        sn2_mat = diag(sn2/sn2div);
    end

    M = K/sn2div+sn2_mat;       % B matrix
else
    if isscalar(sn2)
        sn2_mat = sn2*eye(n);
    else
        sn2_mat = diag(sn2);
    end
    M = K+sn2_mat;           % Covariance with noise
end

[L,p] = chol(M);                % Try computing Cholesky factor

if p~=0     % Failed Cholesky decomposition, compute nearest SPD matrix
    if kmax <= 0; error('Cannot compute Cholesky decomposition.'); end
    % Mold = M;

    M = (M + M')/2;         % Ensure M is symmetric
    [U,Sigma,V] = svd(M);
    H = V*Sigma*V';         % Symmetric polar factor H is SPD
    M = (M+H)/2;
    M = (M + M')/2;         % Ensure symmetry again
    [L,p] = chol(M);        % Retry Cholesky decomposition
    k = 0;
    while p ~= 0            % Failed again, add a small diagonal component
        k = k + 1;
        % We do not want to change the matrix too much
        if k > kmax-1; error('Cannot compute Cholesky decomposition.'); end
        lambda = eig(M);
        mineig = min(lambda);
        kappa = 0.05;        % Sometimes even a small nudge is sufficient
        M = M + kappa*abs(mineig)*k.^2*eye(size(M));
        [L,p] = chol(M);
    end

    % [sort(eig(M))'; sort(eig(Mold))']

end

if Lchol
    sl = sn2div;
    pL = L;                         % L = chol(eye(n)+sW*sW'.*K)
else
    sl = 1;
    pL = -solve_chol(L,eye(n));     % L = -inv(K+inv(sW^2))
end

alpha = solve_chol(L,y-m)/sl;

post.alpha = alpha;                            % return the posterior parameters
post.sW = ones(n,1)/sqrt(min(sn2));            % sqrt of noise precision vector
post.L = pL;
post.Lchol = Lchol;

if nargout>1                               % do we want the marginal likelihood?
  nlZ = (y-m)'*alpha/2 + sum(log(diag(L))) + n*log(2*pi*sl)/2;   % -log marg lik
  if nargout>2                                         % do we want derivatives?
    dnlZ = hyp;                                 % allocate space for derivatives
    Q = solve_chol(L,eye(n))/sl - alpha*alpha';     % precompute for convenience
    for i = 1:numel(hyp.cov)
      dnlZ.cov(i) = sum(sum(Q.*dK(:,:,i)))/2;
    end
    dnlZ.lik = sn2_base*trace(Q);
    for i = 1:numel(hyp.mean)
      dnlZ.mean(i) = -feval(mean{:}, hyp.mean, x, i)'*alpha;
    end
  end
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [x,state] = struct_sum(z,task,dim)
%STRUCT_SUM Sum of elements.
%   [x,state] = struct_sum(z,[],dim) computes x as sum(z,dim), or as
%   sum(z(:)) if dim is the empty matrix [] or not supplied. The structure
%   state stores information which is reused in computing the right and
%   left Jacobian-vector products.
%
%   struct_sum(z,task,dim) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_prod.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
if nargin < 3, dim = []; end
state = [];

if isempty(task)
    if isempty(dim), x = sum(z(:));
    else x = sum(z,dim);
    end
elseif ~isempty(task.r)
    if isempty(dim), x = sum(task.r(:));
    else x = sum(task.r,dim);
    end
elseif ~isempty(task.l)
    if isempty(dim), siz = size(z);
    else siz = ones(1,ndims(z)); siz(dim) = size(z,dim);
    end
    x = repmat(task.l,siz);
end

end

author: andrewssobral

repo: lrslibrary


text:
classdef Plot < handle
%Plot OlyDat Browser Plot
%   A Plot object represents a plot/analysis in the OlyDat Browser. A plot
%   has various responsibilities (all except for the actual plotting are
%   "as appropriate"):
%
%   * Custom data preprocessing
%   * The actual plotting
%   * Interactivity: sending and responding to experiment-selection signals
%   * Math/statistical analysis
%
%   See also OlyDat.Plot/preprocessData, OlyDat.Plot/doPlot,
%   OlyDat.Browser, OlyDat.BrowserStat, OlyDat.ExperimentCoordinator.

    properties (Abstract)
        Name;
        AllowsMultiStats; % (logical) If true, this plot handles multiple statistics
        UsesAuxVar; % Either true, false, or 'custom'.
                    % If true: this plot uses the auxiliary pulldown, as configured at Browser startup.
                    % If false: this plot does not use the aux pulldown. (It will be inactivated.)
                    % If 'custom': This plot configures the auxiliary pulldown itself (using groupVarChanged()).
		    %	           pltCfg.auxVar will be empty when doPlot() is called. This class can use
		    % 		   auxVarChanged() to cache the current auxvar setting.
		    %		   TODO Requiring this class to cache the the current auxvar is unnecessary.

        UsesGroupVar2; % Either true or false
        GroupVarLabel; % Text label for Group Var 1
        GroupVar2Label;

        AutoReplotAfterStatChange; % (logical) If true, stat selection causes replot
        AutoReplotAfterGroupChange;
    end

    methods

        function data = preprocessData(obj,data,plotCfg) %#ok<INUSD,MANU>
            %data = preprocessData(obj,data,plotCfg)
            %   Plots should override this method if they need custom data
            %   manipulation/flattening for their plot. Typical uses
            %   include flattening statistics in the .Scores field to the
            %   top level, averaging scores over tubes, normalizing stats,
            %   etc.
            %
            %   plotCfg is a struct with fields .stat, .grp, .grpVals,
            %   .aux.

            % default implementation: no-op
        end

        function groupChanged(obj,grpStat,allData,hCtrls) %#ok<INUSD,MANU>
            % Called when user selects a new grouping variable. (When the
            % plot is first selected, the current grouping variable is
            % considered to be selected.)
            %
            % grpStat: OlyDat.BrowserStat for new grouping var.
            % allData: struct array, all data loaded in browser.
            % hCtrls: scalar struct:
            %         * hCtrls.pumAuxVar is the auxiliary var pulldown.

            % This method is used when a plot needs to eg customize what
            % auxiliary vars are available, depending on the grouping var.
            % Default implementation: no-op
        end

        function auxVarChanged(obj,auxVar,allData,hCtrls) %#ok<MANU,INUSD>
            % Called when a user selects a new auxiliary variable. (When
            % the plot is first selected, the current aux variable is
            % considered to be selected.)
            %
            % auxVar: OlyDat.BrowserStat for new aux var. if
            % UsesAuxVar=='custom', auxVar will be [].
            % allData: struct array, all data loaded in browser.
            % hCtrls: scalar struct:
            %           * hCtrls.pumAuxVar is the auxiliary var pulldown.

            % Default impl: no-op
        end

    end

    methods (Abstract)
        %detailStr = doPlot(obj,ax,data,stat,grp,plotCfg,expCoordinator)
        %   * ax: axes object
        %   * data: preprocessed data struct for current dataset
        %   * stat: Nrows x Nstats numeric matrix. As a convenience the
        %   Browser puts the statistic(s) of interest into a numeric array
        %   for plotting. Nrows is the number of elements in the
        %   preprocessed data struct.
        %   * grp: Nrows x 1 column vector, either numeric or cellstr. As a
        %   convenience the Browser puts the value of the grouping variable
        %   (if any) in a column vector for plotting.
        %   * plotCfg: The plot configuration specified by the user within
        %   the Browser. This is a struct with fields .stat, .grp,
        %   .grpVals, and .aux.
        %   * expCoordinator: ExperimentCoordinator object used by Browser.
        %   Plots that want to participate in plot interactivity should
        %   register appropriately with the coordinator.
        %   * detailStr: Plots should return a string/cellstr that will be
        %   displayed in the Browser. Typical uses include a description of
        %   the plot, or math details of any statistics performed, etc.
        %
        % OlyDat calls this method when drawing a plot in an axis. This can
        % occur either in the main Browser window, or in standalone figures
        % (eg in the case of 'Compare plots'). If doPlot() performs any
        % actions that require cleanup, that cleanup should be put in ax's
        % DeleteFcn. The Browser will ensure this method is called when a
        % plot is going away. Unregistering from the expCoordinator is an
        % example of such a cleanup action.
        %
        % TODO: unregistering from the expCoordinator is a common action,
        % we should make this easier for users.
        str = doPlot(obj,ax,data,bstat,grp,pltCfg,expCoordinator)
    end

    %% Utilities

    methods (Static)

      function messageAxisCenter(ax,txt)
        % Put a message in the center of axis ax.

        axes(ax);
        xl = xlim(ax);
        yl = ylim(ax);
        hTxt = text(mean(xl),mean(yl),txt);
        set(hTxt,'horizontalalignment','center','interpreter','none');
      end

    end

end

author: kristinbranson

repo: JAABA


text:
function [errorCode,glvar] = pco_camera_resetmem(glvar,segment)
%reset memory size to 0 of all or distinct segment
%
%   [errorCode,glvar] = pco_camera_resetmem(glvar,segment)
%
% * Input parameters :
%    glvar                   structure to hold status info
%    segment                 segment to reset
% * Output parameters :
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%    glvar                   structure to hold status info
%
%structure glvar is used to set different modes for
%load/unload library
%open/close camera SDK
%
%glvar.do_libunload: 1 unload lib at end
%glvar.do_close:     1 close camera SDK at end
%glvar.camera_open:  open status of camera SDK
%glvar.out_ptr:      libpointer to camera SDK handle
%
%if glvar does not exist,
%the library is loaded at begin and unloaded at end
%the SDK is opened at begin and closed at end
%
%if segment does not exist all camera memory is set to segment 1
%
%function workflow
%parameters are checked
%set camera memory structure
%return

% Test if library is loaded
if (~libisloaded('PCO_CAM_SDK'))
    % make sure the dll and h file specified below resides in your current
    % folder
	loadlibrary('SC2_Cam','SC2_CamMatlab.h' ...
                ,'addheader','SC2_CamExport.h' ...
                ,'alias','PCO_CAM_SDK');
	disp('PCO_CAM_SDK library is loaded!');
end

if((exist('glvar','var'))&& ...
   (isfield(glvar,'do_libunload'))&& ...
   (isfield(glvar,'do_close'))&& ...
   (isfield(glvar,'camera_open'))&& ...
   (isfield(glvar,'out_ptr')))
 unload=glvar.do_libunload;
 cam_open=glvar.camera_open;
 do_close=glvar.do_close;
else
 unload=1;
 cam_open=0;
 do_close=1;
end

ph_ptr = libpointer('voidPtrPtr');

%libcall PCO_OpenCamera
if(cam_open==0)
 [errorCode,out_ptr] = calllib('PCO_CAM_SDK', 'PCO_OpenCamera', ph_ptr, 0);
 if(errorCode == 0)
  disp('PCO_OpenCamera done');
  cam_open=1;
  if((exist('glvar','var'))&& ...
     (isfield(glvar,'camera_open'))&& ...
     (isfield(glvar,'out_ptr')))
   glvar.camera_open=1;
   glvar.out_ptr=out_ptr;
  end
 else
   pco_errdisp('PCO_OpenCamera',errorCode);
  if(unload)
   unloadlibrary('PCO_CAM_SDK');
   disp('PCO_CAM_SDK unloadlibrary done');
  end
  return ;
 end
else
 if(isfield(glvar,'out_ptr'))
  out_ptr=glvar.out_ptr;
 end
end

pco_camera_load_defines();

%get Camera Description to test if camera has internal memory
cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
pco_errdisp('PCO_GetCameraDescription',errorCode);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_RECORDER)~=0)
 errorCode=PCO_ERROR_APPLICATION_FUNCTIONNOTSUPPORTED;
 disp('function is not supported because the selected camera does not have internal memory');
 return;
end

act_recstate = uint16(10);
[errorCode,~,act_recstate] = calllib('PCO_CAM_SDK', 'PCO_GetRecordingState',out_ptr,act_recstate);
pco_errdisp('PCO_GetRecordingState',errorCode);

if(act_recstate~=0)
 errorCode=PCO_ERROR_APPLICATION_FUNCTIONNOTSUPPORTED;
 disp('function is not supported while camera is in recording state ON');
 return;
end

dwSegment=uint32(zeros(1,4));

if(exist('segment','var'))
 if((segment<1)||(segment>4))
  errorCode=PCO_ERROR_APPLICATION|PCO_ERROR_WRONGVALUE;
  disp('value of segment must be in the range from 1 to 4');
  return;
 end

 wActSeg=uint16(0);
 [errorCode,~,wActSeg]=calllib('PCO_CAM_SDK', 'PCO_GetActiveRamSegment',out_ptr,wActSeg);
 if(segment==wActSeg)
  errorCode=PCO_ERROR_APPLICATION|PCO_ERROR_WRONGVALUE;
  disp('cannot clear active segment ',int2str(segment));
  return;
 end

 disp(['clear segment ',int2str(segment)]);
 if(errorCode==PCO_NOERROR)
  [errorCode,~,dwSegment]=calllib('PCO_CAM_SDK', 'PCO_GetCameraRamSegmentSize',out_ptr,dwSegment);
  pco_errdisp('PCO_GetCameraRamSegmentSize',errorCode);
 end
 dwSegment(segment)=0;
else
 disp('reset all segments to default state');
end

if((dwSegment(1)==0)&&(dwSegment(2)==0)&&(dwSegment(3)==0)&&(dwSegment(4)==0))
 dwRamSize=uint32(0);
 wPageSize=uint16(0);
 if(errorCode==PCO_NOERROR)
  [errorCode,~,dwRamSize]=calllib('PCO_CAM_SDK', 'PCO_GetCameraRamSize',out_ptr,dwRamSize,wPageSize);
  pco_errdisp('PCO_GetCameraRamSize',errorCode);
 end
 dwSegment(1)=dwRamSize;
 if(errorCode==PCO_NOERROR)
  [errorCode,~]=calllib('PCO_CAM_SDK', 'PCO_SetActiveRamSegment',out_ptr,1);
  pco_errdisp('PCO_SetActiveRamSegment',errorCode);
 end
end

if(errorCode==PCO_NOERROR)
 [errorCode,~]=calllib('PCO_CAM_SDK', 'PCO_SetCameraRamSegmentSize',out_ptr,dwSegment);
 pco_errdisp('PCO_SetCameraRamSegmentSize',errorCode);
end

if(errorCode==PCO_NOERROR)
 [errorCode,~]=calllib('PCO_CAM_SDK', 'PCO_ArmCamera',out_ptr);
 pco_errdisp('PCO_ArmCamera',errorCode);
end

if(exist('glvar','var'))
 glvar = close_camera(out_ptr,unload,do_close,cam_open,glvar);
else
 close_camera(out_ptr,unload,do_close,cam_open);
end

end

function [glvar] = close_camera(out_ptr,unload,do_close,cam_open,glvar)
 if((do_close==1)&&(cam_open==1))
  errorCode = calllib('PCO_CAM_SDK', 'PCO_CloseCamera',out_ptr);
  if(errorCode)
   pco_errdisp('PCO_CloseCamera',errorCode);
  else
   disp('PCO_CloseCamera done');
   cam_open=0;
   if((exist('glvar','var'))&& ...
      (isfield(glvar,'camera_open'))&& ...
      (isfield(glvar,'out_ptr')))
    glvar.out_ptr=[];
    glvar.camera_open=0;
   end
  end
 end
 if((unload==1)&&(cam_open==0))
  unloadlibrary('PCO_CAM_SDK');
  disp('PCO_CAM_SDK unloadlibrary done');
  commandwindow;
 end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function gf=gradloglikGaPExp(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt_r, cx, cy, Wxk] = loglikGaPreadparam(x,varargin);

Hkt=exp(Hkt_r); %nonnegativity constrains
[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);

P=Wxkbg*Hktbg; %current approximation

%linear grasdient shifted by cx
xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx');
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk;
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% gfHkt=(alpha-1)*1./Hkt - 1/beta + Wxk'*(Vxt./P)-1;
gfHkt= Hkt.*((Wxk'*(Vxt./P)-1)-(alpha-1)./Hkt - ones(size(Hkt))*1/beta); %without background (->not Wxkgb) and d(log(L)/dcx)
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');
% d(log(L))/dcy:
gfcy=diag(Wxtcy'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');

gf = [reshape(gfHkt,1,peval.nt*peval.ncomp), gfcx', gfcy'];
% gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];
% gf = [gfcx', gfcy'];
gf=-gf; %conjugate gradient is minimizing!
end

function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
        xxp=double(xx(sizevec, 'corner')); %linear function - pixels
    case 'yy'
        xxp=double(yy(sizevec, 'corner')); %linear function - pixels
    otherwise
        error('Wrong dir')
end
xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
end