text:
function M = spherefactory(n, m, gpuflag)
% Returns a manifold struct to optimize over unit-norm vectors or matrices.
%
% function M = spherefactory(n)
% function M = spherefactory(n, m)
% function M = spherefactory(n, m, gpuflag)
%
% Manifold of n-by-m real matrices of unit Frobenius norm.
% By default, m = 1, which corresponds to the unit sphere in R^n. The
% metric is such that the sphere is a Riemannian submanifold of the space
% of nxm matrices with the usual trace inner product, i.e., the usual
% metric.
%
% Set gpuflag = true to have points, tangent vectors and ambient vectors
% stored on the GPU. If so, computations can be done on the GPU directly.
%
% See also: obliquefactory spherecomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   Oct. 8, 2016 (NB)
%       Code for exponential was simplified to only treat the zero vector
%       as a particular case.
%
%   Oct. 22, 2016 (NB)
%       Distance function dist now significantly more accurate for points
%       within 1e-7 and less from each other.
%
%   July 20, 2017 (NB)
%       Following conversations with Bruno Iannazzo and P.-A. Absil,
%       the distance function is now even more accurate.
%
%   Sep. 7, 2017 (NB)
%       New isometric vector transport available in M.isotransp,
%       contributed by Changshuo Liu.
%
%   April 17, 2018 (NB)
%       ehess2rhess: Used to compute projection of ehess, then subtract a
%       multiple of u (which is assumed tangent.) Now, similarly to what
%       happens in stiefelfactory, we first subtract the multiple of u from
%       ehess, then we project. Mathematically, these operations are the
%       same. Numerically, the former version used to be better because tCG
%       in trustregions had some drift near fine convergence. Now that the
%       drift in tCG has been fixed, it is reasonable to apply the
%       projection last, to ensure best tangency of the output.
%
%   July 18, 2018 (NB)
%       Added the inverse retraction (M.invretr) for the sphere.
%
%   Aug. 3, 2018 (NB)
%       Added GPU support: just set gpuflag = true.
%
%   Jan. 8, 2021 (NB)
%       Added tangent2ambient/tangent2ambient_is_identity pair.

    if ~exist('m', 'var') || isempty(m)
        m = 1;
    end
    if ~exist('gpuflag', 'var') || isempty(gpuflag)
        gpuflag = false;
    end

    % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones)
    % are created directly on the GPU; otherwise, they are created in the
    % usual way (in double precision).
    if gpuflag
        array_type = 'gpuArray';
    else
        array_type = 'double';
    end

    if m == 1
        M.name = @() sprintf('Sphere S^%d', n-1);
    else
        M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
    end

    M.dim = @() n*m-1;

    M.inner = @(x, d1, d2) d1(:)'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @dist;
    function d = dist(x, y)

        % The following code is mathematically equivalent to the
        % computation d = acos(x(:)'*y(:)) but is much more accurate when
        % x and y are close.

        chordal_distance = norm(x - y, 'fro');
        d = real(2*asin(.5*chordal_distance));

        % Note: for x and y almost antipodal, the accuracy is good but not
        % as good as possible. One way to improve it is by using the
        % following branching:
        % % if chordal_distance > 1.9
        % %     d = pi - dist(x, -y);
        % % end
        % It is rarely necessary to compute the distance between
        % almost-antipodal points with full accuracy in Manopt, hence we
        % favor a simpler code.

    end

    M.typicaldist = @() pi;

    M.proj = @(x, d) d - x*(x(:)'*d(:));

    M.tangent = M.proj;

    M.tangent2ambient_is_identity = true;
    M.tangent2ambient = @(X, U) U;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(x, egrad, ehess, u)
        rhess = M.proj(x, ehess - (x(:)'*egrad(:))*u);
    end

    M.exp = @exponential;

    M.retr = @retraction;
    M.invretr = @inverse_retraction;

    M.log = @logarithm;
    function v = logarithm(x1, x2)
        v = M.proj(x1, x2 - x1);
        di = M.dist(x1, x2);
        % If the two points are "far apart", correct the norm.
        if di > 1e-6
            nv = norm(v, 'fro');
            v = v * (di / nv);
        end
    end

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() random(n, m, array_type);

    M.randvec = @(x) randomvec(n, m, x, array_type);

    M.zerovec = @(x) zeros(n, m, array_type);

    M.lincomb = @matrixlincomb;

    M.transp = @(x1, x2, d) M.proj(x2, d);

    % Isometric vector transport of d from the tangent space at x1 to x2.
    % This is actually a parallel vector transport, see Ch. 5 in
    % http://epubs.siam.org/doi/pdf/10.1137/16M1069298
    % "A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization
    %  on Manifolds", by Hosseini and Uschmajew, SIOPT 2017
    M.isotransp = @(x1, x2, d) isometricTransp(x1, x2, d);
    function Td = isometricTransp(x1, x2, d)
        v = logarithm(x1, x2);
        dist_x1x2 = norm(v, 'fro');
        if dist_x1x2 > 0
            u = v / dist_x1x2;
            utd = u(:)'*d(:);
            Td = d + (cos(dist_x1x2)-1)*utd*u ...
                    -  sin(dist_x1x2)  *utd*x1;
        else
            % x1 == x2, so the transport is identity
            Td = d;
        end
    end

    M.pairmean = @pairmean;
    function y = pairmean(x1, x2)
        y = x1+x2;
        y = y / norm(y, 'fro');
    end

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
    M.vecmatareisometries = @() true;

    % Automatically convert a number of tools to support GPU.
    if gpuflag
        M = factorygpuhelper(M);
    end

end

% Exponential on the sphere
function y = exponential(x, d, t)

    if nargin == 2
        % t = 1
        td = d;
    else
        td = t*d;
    end

    nrm_td = norm(td, 'fro');
    y = x*cos(nrm_td) + td*sinxoverx(nrm_td);

end

% Retraction on the sphere
function y = retraction(x, d, t)

    if nargin == 2
        % t = 1;
        td = d;
    else
        td = t*d;
    end

    y = x + td;
    y = y / norm(y, 'fro');

end

% Given x and y two points on the manifold, if there exists a tangent
% vector d at x such that Retr_x(d) = y, this function returns d.
function d = inverse_retraction(x, y)

    % Since
    %   x + d = y*||x + d||
    % and x'd = 0, multiply the above by x' on the left:
    %   1 + 0 = x'y * ||x + d||
    % Then solve for d:

    d = y/(x(:)'*y(:)) - x;

end

% Uniform random sampling on the sphere.
function x = random(n, m, array_type)

    x = randn(n, m, array_type);
    x = x / norm(x, 'fro');

end

% Random normalized tangent vector at x.
function d = randomvec(n, m, x, array_type)

    d = randn(n, m, array_type);
    d = d - x*(x(:)'*d(:));
    d = d / norm(d, 'fro');

end

author: NicolasBoumal

repo: manopt


text:
function c = connectpoly(x,y)
%CONNECTPOLY Connects the vertices of a polygon with straight lines.
%   C = CONNECTPOLY(X,Y) connects the points in X and Y with straight
%   lines, where X and Y are vectors containing the row and column
%   coordinates of the polygon vertices. These points are assumed to be
%   ordered in the clockwise or counterclockwise direction. The output,
%   C, is np-by-2 matrix whose rows are the (row,col) coordinates of the
%   boundary of the polygon, which are in the same direction as the
%   input vertices. The last point in the sequence is equal to the
%   first, thus producing a polygon that is both fully connected and
%   closed.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

v = [x(:),y(:)];

% Close the polygon.
if ~isequal(v(end,:),v(1,:))
   v(end + 1,:) = v(1,:);
end

% Connect the vertices.
segments = cell(1,length(v) - 1);
for k = 2:length(v)
   [x,y] = intline(v(k - 1,1),v(k,1),v(k - 1,2),v(k,2));
   segments{k - 1} = [x,y];
end

c = cat(1,segments{:});

author: dipum

repo: dipum-toolbox

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


text:
function VO = spm_write_sn(V,prm,flags,extras)
% Write out warped images
% FORMAT VO = spm_write_sn(V,prm,flags,msk)
% V         - Images to transform (filenames or volume structure).
% prm       - Transformation information (filename or structure).
% flags     - flags structure, with fields...
%           interp   - interpolation method (0-7)
%           wrap     - wrap edges (e.g., [1 1 0] for 2D MRI sequences)
%           vox      - voxel sizes (3 element vector - in mm)
%                      Non-finite values mean use template vox.
%           bb       - bounding box (2x3 matrix - in mm)
%                      Non-finite values mean use template bb.
%           preserve - either 0 or 1.  A value of 1 will "modulate"
%                      the spatially normalised images so that total
%                      units are preserved, rather than just
%                      concentrations.
%           prefix   - Prefix for normalised images. Defaults to 'w'.
% msk       - An optional cell array for masking the spatially
%             normalised images (see below).
%
% Warped images are written prefixed by "w".
%
% Non-finite vox or bounding box suggests that values should be derived
% from the template image.
%
% Don't use interpolation methods greater than one for data containing
% NaNs.
%__________________________________________________________________________
%
% FORMAT msk = spm_write_sn(V,prm,flags,'mask')
% V          - Images to transform (filenames or volume structure).
% prm        - Transformation information (filename or structure).
% flags      - flags structure, with fields...
%            wrap - wrap edges (e.g., [1 1 0] for 2D MRI sequences)
%            vox  - voxel sizes (3 element vector - in mm)
%                   Non-finite values mean use template vox.
%            bb   - bounding box (2x3 matrix - in mm)
%                   Non-finite values mean use template bb.
% msk        - a cell array for masking a series of spatially normalised
%              images.
%
%
%_________________________________________________________________________
%
% FORMAT VO = spm_write_sn(V,prm,'modulate')
% V         - Spatially normalised images to modulate (filenames or
%             volume structure).
% prm       - Transformation information (filename or structure).
%
%  After nonlinear spatial normalization, the relative volumes of some
%  brain structures will have decreased, whereas others will increase.
%  The resampling of the images preserves the concentration of pixel
%  units in the images, so the total counts from structures that have
%  reduced volumes after spatial normalization will be reduced by an
%  amount proportional to the volume reduction.
%
%  This routine rescales images after spatial normalization, so that
%  the total counts from any structure are preserved.  It was written
%  as an optional step in performing voxel based morphometry.
%
%__________________________________________________________________________
% Copyright (C) 1996-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_write_sn.m 4878 2012-08-31 12:07:30Z guillaume $

SVNid = '$Rev: 4878 $';

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

%-Get images
%--------------------------------------------------------------------------
if isempty(V), return; end

if ischar(prm), prm = load(prm);  end
if ischar(V),   V   = spm_vol(V); end

%-Modulate
%--------------------------------------------------------------------------
if nargin==3 && ischar(flags) && strcmpi(flags,'modulate')
    if nargout==0
        modulate(V,prm);
    else
        VO = modulate(V,prm);
    end
    return;
end

%-Get parameters
%--------------------------------------------------------------------------
def_flags        = spm_get_defaults('old.normalise.write');

if nargin < 3
    flags = def_flags;
else
    fnms = fieldnames(def_flags);
    for i=1:length(fnms)
        if ~isfield(flags,fnms{i})
            flags.(fnms{i}) = def_flags.(fnms{i});
        end
    end
end

%-Mask
%--------------------------------------------------------------------------
[x,y,z,mat] = get_xyzmat(prm,flags.bb,flags.vox);

if nargin==4
    if ischar(extras) && strcmpi(extras,'mask')
        VO = get_snmask(V,prm,x,y,z,flags.wrap);
        return;
    end
    if iscell(extras)
        msk = extras;
    end
end

if nargout>0 && numel(V)>8
    error('Too many images to save in memory');
end

%-Warp images
%--------------------------------------------------------------------------
if ~exist('msk','var')
    msk = get_snmask(V,prm,x,y,z,flags.wrap);
end

if nargout==0
    if isempty(prm.Tr)
        affine_transform(V,prm,x,y,z,mat,flags,msk);
    else
        nonlin_transform(V,prm,x,y,z,mat,flags,msk);
    end
else
    if isempty(prm.Tr)
        VO = affine_transform(V,prm,x,y,z,mat,flags,msk);
    else
        VO = nonlin_transform(V,prm,x,y,z,mat,flags,msk);
    end
end

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

%==========================================================================
%-function VO = affine_transform(V,prm,x,y,z,mat,flags,msk)
%==========================================================================
function VO = affine_transform(V,prm,x,y,z,mat,flags,msk)

[X,Y] = ndgrid(x,y);
d     = [flags.interp*[1 1 1]' flags.wrap(:)];

spm_progress_bar('Init',numel(V),'Resampling','volumes/slices completed');
for i=1:numel(V)
    VO     = make_hdr_struct(V(i),x,y,z,mat, flags.prefix);
    if flags.preserve
        VO.fname = spm_file(VO.fname, 'prefix','m');
    end
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    if flags.preserve, VO.pinfo(1:2,:) = VO.pinfo(1:2,:)/detAff; end

    %Dat= zeros(VO.dim(1:3));
    Dat = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    C = spm_bsplinc(V(i),d);

    for j=1:length(z)   % Cycle over planes
        [X2,Y2,Z2]  = mmult(X,Y,z(j),V(i).mat\prm.VF(1).mat*prm.Affine);
        dat         = spm_bsplins(C,X2,Y2,Z2,d);
        if flags.preserve, dat = dat*detAff; end
        dat(msk{j}) = NaN;

        Dat(:,:,j) = single(dat);

        if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end
    end
    if nargout~=0
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    else
        spm_write_vol(VO, Dat);
    end
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

%==========================================================================
%-function VO = nonlin_transform(V,prm,x,y,z,mat,flags,msk)
%==========================================================================
function VO = nonlin_transform(V,prm,x,y,z,mat,flags,msk)

[X,Y] = ndgrid(x,y);
Tr = prm.Tr;
BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
if flags.preserve
    DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
    DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
    DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');
end
d  = [flags.interp*[1 1 1]' flags.wrap(:)];

spm_progress_bar('Init',numel(V),'Resampling','volumes completed');
for i=1:numel(V)
    VO     = make_hdr_struct(V(i),x,y,z,mat, flags.prefix);
    if flags.preserve
        VO.fname = spm_file(VO.fname, 'prefix','m');
    end
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);

    % Accumulate data
    %Dat= zeros(VO.dim(1:3));
    Dat = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    C = spm_bsplinc(V(i),d);

    for j=1:length(z)   % Cycle over planes
        % Nonlinear deformations
        %------------------------------------------------------------------
        tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
        ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
        tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);
        X1 = X    + BX*tx*BY';
        Y1 = Y    + BX*ty*BY';
        Z1 = z(j) + BX*tz*BY';

        [X2,Y2,Z2]  = mmult(X1,Y1,Z1,V(i).mat\prm.VF(1).mat*prm.Affine);
        dat         = spm_bsplins(C,X2,Y2,Z2,d);
        dat(msk{j}) = NaN;

        if ~flags.preserve
            Dat(:,:,j) = single(dat);
        else
            j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
            j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
            j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;

            % The determinant of the Jacobian reflects relative volume changes
            %-----------------------------------------------------------------
            dat       = dat .* (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
            Dat(:,:,j) = single(dat);
        end
        if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end
    end
    if nargout==0
        if flags.preserve, VO = rmfield(VO,'pinfo'); end
        VO = spm_write_vol(VO,Dat);
    else
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    end
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

%==========================================================================
%-function VO = modulate(V,prm)
%==========================================================================
function VO = modulate(V,prm)

spm_progress_bar('Init',numel(V),'Modulating','volumes completed');
for i=1:numel(V)
    VO          = V(i);
    VO          = rmfield(VO,'pinfo');
    VO.fname    = spm_file(VO.fname, 'prefix','m');
    detAff      = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    %Dat        = zeros(VO.dim(1:3));
    Dat         = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;
    [x,y,z,mat] = get_xyzmat(prm,NaN,NaN,VO);

    if sum((mat(:)-VO.mat(:)).^2)>1e-7, error('Orientations not compatible'); end

    Tr = prm.Tr;

    if isempty(Tr)
        for j=1:length(z),   % Cycle over planes
            dat        = spm_slice_vol(V(i),spm_matrix([0 0 j]),V(i).dim(1:2),0);
            Dat(:,:,j) = single(dat);
            if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end
        end
    else
        BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
        BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
        BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
        DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
        DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
        DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');

        for j=1:length(z)   % Cycle over planes

            tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
            ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
            tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);

            j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
            j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
            j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;

            % The determinant of the Jacobian reflects relative volume changes
            %-----------------------------------------------------------------
            dat        = spm_slice_vol(V(i),spm_matrix([0 0 j]),V(i).dim(1:2),0);
            dat        = dat .* (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
            Dat(:,:,j) = single(dat);

            if numel(V)<5, spm_progress_bar('Set',i-1+j/length(z)); end
        end
    end

    if nargout==0
        VO = spm_write_vol(VO,Dat);
    else
        VO.pinfo  = [1 0]';
        VO.dt     = [spm_type('float32') spm_platform('bigend')];
        VO.dat    = Dat;
    end
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

%==========================================================================
%-function VO   = make_hdr_struct(V,x,y,z,mat,prefix)
%==========================================================================
function VO   = make_hdr_struct(V,x,y,z,mat,prefix)
VO            = V;
VO.fname      = spm_file(V.fname, 'prefix',prefix);
VO.mat        = mat;
VO.dim(1:3)   = [length(x) length(y) length(z)];
VO.pinfo      = V.pinfo;
VO.descrip    = 'spm - 3D normalized';

%==========================================================================
%-function T2 = get_2Dtrans(T3,B,j)
%==========================================================================
function T2 = get_2Dtrans(T3,B,j)
d   = [size(T3) 1 1 1];
tmp = reshape(T3,d(1)*d(2),d(3));
T2  = reshape(tmp*B(j,:)',d(1),d(2));

%==========================================================================
%-function Mask = getmask(X,Y,Z,dim,wrp)
%==========================================================================
function Mask = getmask(X,Y,Z,dim,wrp)
% Find range of slice
tiny = 5e-2;
Mask = true(size(X));
if ~wrp(1), Mask = Mask & (X >= (1-tiny) & X <= (dim(1)+tiny)); end
if ~wrp(2), Mask = Mask & (Y >= (1-tiny) & Y <= (dim(2)+tiny)); end
if ~wrp(3), Mask = Mask & (Z >= (1-tiny) & Z <= (dim(3)+tiny)); end

%==========================================================================
%-function [X2,Y2,Z2] = mmult(X1,Y1,Z1,Mult)
%==========================================================================
function [X2,Y2,Z2] = mmult(X1,Y1,Z1,Mult)
if length(Z1) == 1
    X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + (Mult(1,3)*Z1 + Mult(1,4));
    Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + (Mult(2,3)*Z1 + Mult(2,4));
    Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + (Mult(3,3)*Z1 + Mult(3,4));
else
    X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
    Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
    Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
end

%==========================================================================
%-function msk = get_snmask(V,prm,x,y,z,wrap)
%==========================================================================
function msk = get_snmask(V,prm,x,y,z,wrap)
% Generate a mask for where there is data for all images
%--------------------------------------------------------------------------
msk = cell(length(z),1);
t1 = cat(3,V.mat);
t2 = cat(1,V.dim);
t  = [reshape(t1,[16 length(V)])' t2(:,1:3)];
Tr = prm.Tr;
[X,Y] = ndgrid(x,y);

BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);

if numel(V)>1 && any(any(diff(t,1,1)))
    spm_progress_bar('Init',length(z),'Computing available voxels','planes completed');
    for j=1:length(z)   % Cycle over planes
        Count = zeros(length(x),length(y));
        if isempty(Tr)
            % Generate a mask for where there is data for all images
            %--------------------------------------------------------------
            for i=1:numel(V),
                [X2,Y2,Z2] = mmult(X,Y,z(j),V(i).mat\prm.VF(1).mat*prm.Affine);
                Count      = Count + getmask(X2,Y2,Z2,V(i).dim(1:3),wrap);
            end
        else
            % Nonlinear deformations
            %--------------------------------------------------------------
            X1 = X    + BX*get_2Dtrans(Tr(:,:,:,1),BZ,j)*BY';
            Y1 = Y    + BX*get_2Dtrans(Tr(:,:,:,2),BZ,j)*BY';
            Z1 = z(j) + BX*get_2Dtrans(Tr(:,:,:,3),BZ,j)*BY';

            % Generate a mask for where there is data for all images
            %--------------------------------------------------------------
            for i=1:numel(V)
                [X2,Y2,Z2] = mmult(X1,Y1,Z1,V(i).mat\prm.VF(1).mat*prm.Affine);
                Count      = Count + getmask(X2,Y2,Z2,V(i).dim(1:3),wrap);
            end
        end
        msk{j} = uint32(find(Count ~= numel(V)));
        spm_progress_bar('Set',j);
    end
    spm_progress_bar('Clear');
else
    for j=1:length(z), msk{j} = uint32([]); end
end

%==========================================================================
%-function [x,y,z,mat] = get_xyzmat(prm,bb,vox,VG)
%==========================================================================
function [x,y,z,mat] = get_xyzmat(prm,bb,vox,VG)
% The old voxel size and origin notation is used here.
% This requires that the position and orientation
% of the template is transverse.  It would not be
% straitforward to account for templates that are
% in different orientations because the basis functions
% would no longer be seperable.  The seperable basis
% functions mean that computing the deformation field
% from the parameters is much faster.

% bb  = sort(bb);
% vox = abs(vox);

if nargin<4
    VG = prm.VG(1);

    if all(~isfinite(bb(:))) && all(~isfinite(vox(:)))
        x   = 1:VG.dim(1);
        y   = 1:VG.dim(2);
        z   = 1:VG.dim(3);
        mat = VG.mat;
        return;
    end
end

[bb0 vox0] = spm_get_bbox(VG, 'old');
if ~all(isfinite(vox(:))), vox = vox0; end
if ~all(isfinite(bb(:))),  bb  = bb0;  end

msk       = find(vox<0);
bb        = sort(bb);
bb(:,msk) = flipud(bb(:,msk));

% Adjust bounding box slightly - so it rounds to closest voxel.
% Comment out if not needed.
%bb(:,1) = round(bb(:,1)/vox(1))*vox(1);
%bb(:,2) = round(bb(:,2)/vox(2))*vox(2);
%bb(:,3) = round(bb(:,3)/vox(3))*vox(3);

M   = prm.VG(1).mat;
vxg = sqrt(sum(M(1:3,1:3).^2));
if det(M(1:3,1:3))<0, vxg(1) = -vxg(1); end
ogn = M\[0 0 0 1]';
ogn = ogn(1:3)';

% Convert range into range of voxels within template image
x   = (bb(1,1):vox(1):bb(2,1))/vxg(1) + ogn(1);
y   = (bb(1,2):vox(2):bb(2,2))/vxg(2) + ogn(2);
z   = (bb(1,3):vox(3):bb(2,3))/vxg(3) + ogn(3);

og  = -vxg.*ogn;

% Again, chose whether to round to closest voxel.
%of  = -vox.*(round(-bb(1,:)./vox)+1);
of = bb(1,:)-vox;

M1  = [vxg(1) 0 0 og(1) ; 0 vxg(2) 0 og(2) ; 0 0 vxg(3) og(3) ; 0 0 0 1];
M2  = [vox(1) 0 0 of(1) ; 0 vox(2) 0 of(2) ; 0 0 vox(3) of(3) ; 0 0 0 1];
mat = prm.VG(1).mat*inv(M1)*M2;

LEFTHANDED = true;
if (LEFTHANDED && det(mat(1:3,1:3))>0) || (~LEFTHANDED && det(mat(1:3,1:3))<0),
    Flp = [-1 0 0 (length(x)+1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
    mat = mat*Flp;
    x   = flipud(x(:))';
end

%==========================================================================
%-function VO = write_dets(P,bb,vox)
%==========================================================================
function VO = write_dets(P,bb,vox)
if nargin==1
    job = P;
    P   = job.P;
    bb  = job.bb;
    vox = job.vox;
end

spm_progress_bar('Init',numel(P),'Writing','volumes completed');

for i=1:numel(V)
    prm = load(deblank(P{i}));
    [x,y,z,mat] = get_xyzmat(prm,bb,vox);
    Tr = prm.Tr;
    BX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1);
    BY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1);
    BZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1);
    DX = spm_dctmtx(prm.VG(1).dim(1),size(Tr,1),x-1,'diff');
    DY = spm_dctmtx(prm.VG(1).dim(2),size(Tr,2),y-1,'diff');
    DZ = spm_dctmtx(prm.VG(1).dim(3),size(Tr,3),z-1,'diff');

    VO = struct(...
        'fname',   spm_file(P{i}, 'prefix','jy_'),...
        'dim',     [numel(x),numel(y),numel(z)],...
        'dt',      [spm_type('float32') spm_platform('bigend')],...
        'pinfo',   [1 0 0]',...
        'mat',     mat,...
        'n',       1,...
        'descrip', 'Jacobian determinants');
    VO     = spm_create_vol(VO);
    detAff = det(prm.VF(1).mat*prm.Affine/prm.VG(1).mat);
    Dat    = single(0);
    Dat(VO.dim(1),VO.dim(2),VO.dim(3)) = 0;

    for j=1:length(z)   % Cycle over planes
        % Nonlinear deformations
        tx = get_2Dtrans(Tr(:,:,:,1),BZ,j);
        ty = get_2Dtrans(Tr(:,:,:,2),BZ,j);
        tz = get_2Dtrans(Tr(:,:,:,3),BZ,j);

        %------------------------------------------------------------------
        j11 = DX*tx*BY' + 1; j12 = BX*tx*DY';     j13 = BX*get_2Dtrans(Tr(:,:,:,1),DZ,j)*BY';
        j21 = DX*ty*BY';     j22 = BX*ty*DY' + 1; j23 = BX*get_2Dtrans(Tr(:,:,:,2),DZ,j)*BY';
        j31 = DX*tz*BY';     j32 = BX*tz*DY';     j33 = BX*get_2Dtrans(Tr(:,:,:,3),DZ,j)*BY' + 1;
        % The determinant of the Jacobian reflects relative volume changes.
        %------------------------------------------------------------------
        dat       = (j11.*(j22.*j33-j23.*j32) - j21.*(j12.*j33-j13.*j32) + j31.*(j12.*j23-j13.*j22)) * detAff;
        Dat(:,:,j) = single(dat);
        if numel(P)<5, spm_progress_bar('Set',i-1+j/length(z)); end
    end
    VO = spm_write_vol(VO,Dat);
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Classes representing Tensors
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function [inside, outside] = find_inside_vol(pos, vol);

% FIND_INSIDE_VOL locates dipole locations inside/outside the source
% compartment of a volume conductor model.
%
% [inside, outside] = find_inside_vol(pos, vol)
%
% This function is obsolete and its use in other functions should be replaced
% by inside_vol

% Copyright (C) 2003-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: find_inside_vol.m 2885 2011-02-16 09:41:58Z roboos $

inside  = ft_inside_vol(pos, vol);
% replace boolean vector with indexing vectors
outside = find(~inside);
inside  = find(inside);

author: PatternRecognition

repo: OpenBMI


text:
function [A,AE] = adjacency_incident_angle_matrix(V,E)
  % ADJACENCY_INCIDENT_ANGLE_MATRIX  Compute an adjacency matrix between *edges*
  % where non-zeros are the angles between adjacent edges (edges that share a
  % vertex).
  %
  % [A] = adjacency_incident_angle_matrix(V,E)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   E  #E by 2 list of edge indices into rows of V
  % Outputs:
  %   A  #E by #E adjacency matrix so that A(i,j) = h ≠ 0 implies that edges
  %     E(i,:) and E(j,:) meet at a coincident vertex with angle h.
  %   AE  #E by #E adjacency_matrix so that A(i,j) ≠ 0 implies that edges
  %     E(i,:) and E(j,:) meet at a coincident vertex (just in case h==0)
  %

  assert(size(E,1)==size(unique(sort(E,2),'rows'),1));

  V2E = sparse(E,repmat(1:size(E,1),2,1)',1,size(V,1),size(E,1));
  spy(V2E)
  AE = V2E'*V2E;
  spy(AE)
  [I,J] = find(triu(AE,1));
  C = sparse(repmat(1:numel(I),4,1)',[E(I,:) E(J,:)],1,numel(I),size(V,1));
  MI = (1:numel(I))';
  % get center vertex
  [~,M2] = max(C==2,[],2);
  C(sub2ind(size(C),MI,M2)) = 0;
  % get other vertex (who cares which order)
  [~,M1] = max(C,[],2);
  C(sub2ind(size(C),MI,M1)) = 0;
  % get other vertex (who cares which order)
  [~,M3] = max(C,[],2);
  V12 = normalizerow(V(M1,:)-V(M2,:));
  V23 = normalizerow(V(M3,:)-V(M2,:));
  H = acos(sum(V12.*V23,2));
  % rebuild A
  A = sparse([I J],[J I],[H H],size(E,1),size(E,1));
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% Set path to iFEM
%
% Add all subdirectories under pathstr to search path.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

addpath(genpath(pwd),'-begin');
savepath;

repo: Aero-matlab


text:
function B = biharmonic(f)
%BIHARMONIC   Biharmonic operator of a SEPARABLEAPPROX.
%   B = BIHARMONIC(F) returns a SEPARABLEAPPROX representing the biharmonic
%   operator applied to F.
%
% See also SEPARABLEAPPROX/BIHARM.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% biharmonic(f) = f_xxxx + f_yyyy + 2*f_xxyy:
B = diff(f, 4, 2) + diff(f, 4, 1) + 2*diff(diff(f, 2, 1), 2, 2);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [signal, state] = flt_sourceLocalize(varargin)
% Localize sources for a given head model.
% Estimated sources and inverse operator will be stored in
% signal.CSD and signal.invOperator
%
% Author: Tim Mullen, Jan 2013, SCCN/INC/UCSD
%         Alejandro Ojeda, Jan 2013, SCCN/INC/UCSD
%         Christian Kothe, Jan 2013, SCCN/INC/UCSD

% FIXME: Need to add a for loop over the calls to Loreta (or add it in the
% flt_loreta function) which loops over sub-partitions of the data block,
% estimating CSD for each partition (i.e. assuming local stationarity of
% the inverse operator within each partition). The invOperator results are
% stored in a 3D array and the mean invOperator is then returned in state
% obj as well as in signal obj.
% FIXME: Implement option to only update inverse operatior infrequently,
% otherwise using the inv op stored in state object (i.e. mean over last
% block).

if ~exp_beginfun('filter'), return; end

declare_properties('name','SourceLocalization', 'experimental',true, 'follows',{'flt_ica','flt_fir','flt_iir'}, 'cannot_precede',{'flt_clean_settings'},'independent_channels',false, 'independent_trials',false);

headmodel_default = 'resources:/headmodels/standard-Colin27-385ch.mat';

% try to get ROI names supported by head model (for use in arg definition)
ROINames = {};
if ~onl_isonline
    hmObj = arg_extract(varargin,{'hmObj','HeadModelObject'},[],headmodel_default);
    if isstruct(hmObj) || (ischar(hmObj) && ~isempty(hmObj))
        hmObj = hlp_validateHeadModelObject(hmObj);
        if ~isempty(hmObj)
            % get the unique ROI names
            [tmp, idx] = unique(hmObj.atlas.label); %#ok<ASGLU>
            ROINames = hmObj.atlas.label(sort(idx))';
        end
    end
end

arg_define(varargin, ...
    arg_norep({'signal','Signal'},[],[],'Signal structure. Must contain .data field with channel data'), ...
    arg({'hmObj','HeadModelObject'},headmodel_default,[],'Head model object generated by MOBILAB. Usually a filename, or a cell array of the form {query,file,query,file,...} where a query is a wildcard pattern that, if matching the sourcedata filename, determines that the subsequent file in the list shall be used  See MOBILAB''s headModel class.','type','expression'), ...
    arg({'channels','Channels'},[],[],'Cell array of channel labels to retain. If empty, all channels in signal.chanlocs.labels will be used. Channels must be a subset of hmObj.getChannelLabels()','type','cellstr','shape','row'), ...
    arg_subswitch({'invMethod','InverseMethod'},{'LORETA'}, ...
    [{...
      'LORETA', @flt_loreta, ...
      'Evidence Approximation', @flt_evidence_approx_tikhonov_reg, ...
      'Variational LORETA', @flt_vbloreta, ...
      'Beamforming', @flt_beamforming}, ...
      quickif(exist('flt_champagne','file'),{'Champagne',@flt_champagne},{})], ...
     'Source Localization Algorithm'), ...
    arg({'sourceAtlasLabels','SourceAtlasLabels','SourceAtlasROI'},ROINames,ROINames,'Source regions of interest (atlas labels). This is a cell array of strings corresponding to a subset of the labels stored in hmObj.atlas.label. Current source density will be estimated only for these ROIs.','type','logical','rangecheck',true), ...
    arg({'colRoiCsd','CollapseRoiCsd'},'sum',{'none','mean','sum','max','maxmag','median'},'Method for computing ROI CSD. Return the (mean, integral, max) current source density within each ROI. signal.srcpot matrix will be reduced to [num_rois x num_samples].'), ...
    arg({'keepFullCsd','KeepFullCsd'},true,[],'Keep a copy of the csd for each vertex. This is stored in srcpot_all'), ...
    arg({'roiAtlasLabels','ROIAtlasLabels','AtlasROI','ROI'},ROINames,ROINames,'Regions of interest atlas labels. This is a cell array of strings corresponding to a subset of the labels stored in hmObj.atlas.label. If CollapseRoiCsd is set, current source density will be integrated over these ROIs. If empty the all ROIs in SourceAtlasLabels are selected.','type','logical','rangecheck',true), ...
    arg({'combineROIs','CombineROI'},[],[],'Combine ROIs in ROIAtlasLabels. This is a cell array of (name, value) pairs. The name entry is a string denoting the name of a new ROI. The value part is a cell array of strings indicating which of the ROIs in ROIAtlasLabels should be combined to form a new ROI. Any ROIs that were merged into a new ROI are removed from the original set of ROIs. To preserve an ROI that was merged into a new ROI, specify a new ROI containing only this ROI. All new ROIs are then appended to the list of existing (unmerged) ROIs in ROIAtlasLabels. If CollapseRoiCsd is set, current source density will be integrated over all ROIs.','type','expression','shape','row'), ...
    arg({'combineMode','CombineMode'},'merge',{'replace','append','merge'},'Mode for combining MetaROIs with original ROIs'), ...
    arg_nogui({'appendROI','AppendROI'},false,[],'Append superROIs or replace'), ...
    arg({'reference','Reference'},[],[],'Reference scheme. This is a cell array of channel labels defining the reference. If more than one channel, average reference is assumed','type','cellstr','shape','row'), ...
    arg_nogui({'roiVertices','ROIVertices'},[],[],'Regions of interest (vertices). This is supplementary to the ''RegionsOfInterest'' option. ''roiVertices'' is a cell array where each cell contains the indices of all vertices within a region of interest. The vertex indices correspond to respective columns of the lead field matrix. These ROIs will be added to the set of ROIs defined in the ''RegionsOfInterest'' option.'), ...
    arg_nogui({'roiVerticesLabels','ROIVerticesLabels'},[],[],'Labels for the supplementary ROIs. This is a cell array of same dimension as ''roiVertices'' with labels for each ROI','type','cellstr'), ...
    arg({'normalCSDWeight','NormalCSDWeight'},0.7,[],'Weight of normal direction in CSD. When calculating the CSD at a given vertex from a lead-field matrix that is not restricted to the surface normal this is the weighting of the current normal to the surface (and tangential current will be weighted by 1 minus this value).','guru',true),...
    arg({'leadfieldOrientation','LeadfieldOrientation'},'axisparallel',{'axisparallel','tangentspace'},'Orientation of local leadfield components. This is for a given source vertex (if orientations are flexible). If axisparallel, the projections will point in X/Y/Z directions, if tangentspace, the projections will point in surface normal, tangent and binormal directions.','guru',true),...
    arg({'makeDipfitStruct','MakeDipfitStruct'},true,[],'Make dipfit structure. If selected, a dipfit structure will be created in signal.dipfit containing the locations (.posxyz) and moments (.mom) of each dipole or center of mass of ROI.'), ...
    arg({'do_transform','TransformData','transform'},false,[],'Transform the data rather than annotate. By default, source reconstructions are added as annotations (.srcpot) to the data set.'),...
    arg({'usecar','CAR'},true,[],'Rereference to common average. This should always be enabled, unless you are really sure about it'), ...
    arg({'verb','Verbosity'},0,[],'Verbosity level.','typecheck',false), ...
    arg_nogui({'state','State'}));

if isempty(roiAtlasLabels) %#ok
    roiAtlasLabels = ROINames;
end
if isempty(sourceAtlasLabels) %#ok
    sourceAtlasLabels = ROINames;
end

% Intialization block
% -------------------------------------------------------------------------
if isempty(state)
    % || (isfield(state,'roiAtlasLabels') && ~isequal(state.roiAtlasLabels,roiAtlasLabels))

    if iscellstr(hmObj)
        % per-subject headmodels given as a cell array of query-model pairsl resolve which one to use
        loader = utl_find_filter(signal,'io_loadset');
        filename = loader.parts{1};
        for k=1:2:length(hmObj)
            pattern = regexptranslate('wildcard',hmObj{k});
            if regexp(filename,pattern)
                hmObj = hmObj{k+1};
                break;
            end
        end
        if iscellstr(hmObj)
            error('No matching headmodel found for the given data filename: %s (patterns: %s)',filename,hlp_tostring(hmObj(1:2:end))); end
    end

    % Either we are initializing our adaptive estimator or user-specified
    % ROIs have changed.
    state = hlp_microcache('sourcespace',@build_state,hmObj,channels,sourceAtlasLabels,roiAtlasLabels,roiVertices,roiVerticesLabels,makeDipfitStruct,verb,signal.chanlocs,combineROIs,combineMode,usecar,normalCSDWeight,leadfieldOrientation);
end

% Estimate Current Source Density
% -------------------------------------------------------------------------

% Current source density is returned in signal.srcpot
% [num_vertices x num_samples]
% Inverse operator is returned in signal.srcweights
% [num_vertices x num_channels]

if usecar
    % apply common average reference
    signal.data = bsxfun(@minus,signal.data,mean(signal.data));
end
switch lower(invMethod.arg_selection)
    case 'loreta'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_loreta,                    ...
                                    'signal',signal,                ...
                                    invMethod,                      ...
                                    'state',state.solverState,      ...
                                    'K',state.leadFieldMatrix,      ...
                                    'L',state.laplacianOperator,    ...
                                    'verb',max(0,verb-1),'arg_direct',true);
        signal.srcweights = state.solverState.srcweights;
    case 'evidence approximation'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_evidence_approx_tikhonov_reg,                  ...
                                    'signal',signal,                ...
                                    invMethod,                      ...
                                    'state',state.solverState,      ...
                                    'K',state.leadFieldMatrix,      ...
                                    'L',state.laplacianOperator,    ...
                                    'verb',verb,'arg_direct',true);
        signal.srcweights = state.solverState.srcweights;
    case 'variational loreta'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_vbloreta,                  ...
                                    'signal',signal,                ...
                                    invMethod,                      ...
                                    'state',state.solverState,      ...
                                    'K',state.leadFieldMatrix,      ...
                                    'L',state.laplacianOperator,    ...
                                    'verb',max(0,verb-1),'arg_direct',true);
        signal.srcweights = state.solverState.srcweights;
    case 'dynamic loreta'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_dynamicloreta,             ...
                                    'signal',signal,                ...
                                    invMethod,                      ...
                                    'state',state.solverState,      ...
                                    'K',state.leadFieldMatrix,      ...
                                    'L',state.laplacianOperator,    ...
                                    'verb',max(0,verb-1),'arg_direct',true);
        signal.srcweights = state.solverState.srcweights;
    case 'beamforming'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_beamforming,                  ...
                                    'signal',signal,                ...
                                    invMethod, ...
                                    'lead_field', state.leadFieldMatrix, ...
                                    'state',state.solverState,      ...
                                    'arg_direct',true);
    case 'champagne'
        [signal, state.solverState] = hlp_scope({'disable_expressions',true}, ...
                                    @flt_champagne,                  ...
                                    'signal',signal,                ...
                                    invMethod,                      ...
                                    'state',state.solverState,      ...
                                    'head_file',struct('leadfield',reshape(state.leadFieldMatrix,state.nChannels,state.nVertices,state.orientationDegrees), ...
                                                       'chanlabels',{{signal.chanlocs.labels}}), ...
                                    'arg_direct',true);

    otherwise
        error('flt_sourceLocalize:badInverseMethod','Invalid inverse method %s',invMethod.arg_selection);
end

% collapse srcpot from 3d to 1d based on state.srcweightCollapseDoFWeights
if state.orientationDegrees == 3
    siz = size(signal.srcpot);
    signal.srcpot = reshape(sum(reshape(bsxfun(@times,signal.srcpot,state.srcweightCollapseDoFWeights),[3 siz(1)/3 siz(2:end)]),1),[siz(1)/3 siz(2:end)]);
end

% Collapse CSD within each ROI
% -------------------------------------------------------------------------
% FIXME: We may wish to replace avg with the surface integral over the ROI
% We can use numerical double integration i.e. quad2d() or more crudely:
% sum(signal.CSD(x,:)) * dx * dy;
% where dx, dy are the voxel x,y dimensions.
% If the csd is not constrained to the surface of a mesh, then we need to
% use triple integration i.e. triplequad(). Or
% sum(signal.CSD(x,:)) * dx * dy * dz

if ~strcmp(colRoiCsd,'none') && ~isempty(signal.srcpot)
    if verb
        fprintf('Computing %s CSD for each ROI \n', colRoiCsd);
    end

    if keepFullCsd
        % store backup with all vertices
        % also, register this as a time-series field
        signal = utl_register_field(signal,'timeseries','srcpot_all',signal.srcpot);
        signal.srcweights_all = signal.srcweights;
    end

    % collapse current density
    signal.srcpot = hlp_colsrc(signal.srcpot,state.roiVerticesReduced,colRoiCsd);
    % collapse weights (note: if weights are time-varying, then whatever snapshot this is isn't
    % going to be too useful)
    if isfield(signal,'srcweights')
        signal.srcweights = hlp_colsrc(signal.srcweights,state.roiVerticesReduced,colRoiCsd); end
else
    signal.srcpot_all       = [];
    signal.srcweights_all   = [];
end

% Store additional outputs in signal
% -------------------------------------------------------------------------
if makeDipfitStruct
    signal.dipfit   = state.dipfit;
end
signal.roiLabels          = state.roiLabels;
signal.roiVertices        = state.roiVertices;
signal.roiVerticesReduced = state.roiVerticesReduced;
signal.rmIndices          = state.rmIndices;

if do_transform
    signal.data = signal.srcpot;
    signal.nbchan = size(signal.data,1);
    signal.chanlocs = struct('labels',cellfun(@num2str,num2cell(1:signal.nbchan,1),'UniformOutput',false));
end
signal.hmObj = hmObj;
signal.leadFieldMatrix = state.leadFieldMatrix;

exp_endfun;

function state = build_state(hmObj,channels,sourceAtlasLabels,roiAtlasLabels,roiVertices,roiVerticesLabels,makeDipfitStruct,verb,chanlocs,combineROIs,combineMode,usecar,normalCSDWeight,leadfieldOrientation)
% validate the head model and construct the source space

orilen = length(roiAtlasLabels);
[res ia] = intersect(roiAtlasLabels,sourceAtlasLabels);
roiAtlasLabels = roiAtlasLabels(sort(ia));  % FIXME: double-check this order
if length(roiAtlasLabels)~=orilen
    warn_once('BCILAB:AtlasLabelMismatch','Some elements of ROIAtlasLabels were not contained in SourceAtlasLabels. These were removed.');
end

% Validate Inputs
hmObj     = hlp_validateHeadModelObject(hmObj);
if isempty(hmObj)
    error('HeadModelObject was improperly defined. Exiting');
end
if ~isempty(roiVerticesLabels) && ~isempty(roiVertices) ...
        && length(roiVertices) ~= length(roiVerticesLabels)
    fprintf('WARNING: The number of ROI labels in ''roiVerticesLabels'' does not match the number of ROIs in ''roiVertices''. Reverting to generic labeling scheme {''ROI1'' ''ROI2'' ... ''ROIN''}.');
    roiVerticesLabels = [];
end

if verb
    fprintf('Initializing source reconstructor \n');
end

% Prune lead field matrix and laplacian operator
% to contain only regions of interest
% If there are no specific ROIs selected, this will "open" the
% surface mesh by removing a few vertices
% (the surface mesh cannot be a closed surface)
if verb
    fprintf('Constructing source space \n');
end

% remove undesired brain structures
% note: we assume that the leadFieldMatrix produced by this step is 2d, but
% contain stacked pages for X, Y and Z dipole components (so that a reshape
% to [#ch,#vox,3] would yield the correct result.
brainStructsToRemove = setdiff(unique(hmObj.atlas.label),sourceAtlasLabels);
[   reducedSpace,                      ...
    state.leadFieldMatrix,             ...
    state.laplacianOperator            ...
    state.rmIndices                    ...
    ] = getSourceSpace4PEB(hmObj,brainStructsToRemove,roiVertices);

state.nVertices = size(reducedSpace.vertices,1);
state.nChannels = length(chanlocs);

% determine the number of orientation degrees of freedom
if size(state.leadFieldMatrix,2) == state.nVertices
    state.orientationDegrees = 1;
elseif size(state.leadFieldMatrix,2) == 3*state.nVertices
    state.orientationDegrees = 3;
else
    error('Your lead-field matrix has a # of columns that does not correspond to the # of source vertices.');
end

% project leadfield matrix onto tangent space
if state.orientationDegrees==3
    % get tangent space at each vertex
    if isprop(hmObj,'surfNormal')
        state.surfNormals = blk_diag(hmObj.surfNormal,1)';
        % prune removed indices
        state.surfNormals(state.rmIndices,:)=[];
    else
        % FIXME: these data contains NaN's (presumably near deleted triangles)
        state.surfNormals = geometricTools.getSurfaceNormals(reducedSpace.vertices,reducedSpace.faces,true);
    end
    state.surfTangents = cross(state.surfNormals,ones(length(state.surfNormals),1)*[1 0 0],2);
    state.surfBinormals = cross(state.surfNormals,state.surfTangents,2);

    % also calc final projection operator
    if strcmp(leadfieldOrientation,'axisparallel')
        state.srcweightCollapseDoFWeights = vec((state.surfNormals*normalCSDWeight + (state.surfTangents+state.surfBinormals)*(1-normalCSDWeight))'); % does assume axis-parallel lead-field
    elseif strcmp(leadfieldOrientation,'tangentspace')
        state.srcweightCollapseDoFWeights = vec((ones(length(state.surfNormals),1)*[normalCSDWeight (1-normalCSDWeight)/2 (1-normalCSDWeight)/2])');  % assumes rotated lead-field
        % turn into triplet order
        tmp = permute(reshape(state.leadFieldMatrix,state.nChannels,state.nVertices,state.orientationDegrees),[1 3 2]);
        % transpose for recombination
        tmp = tmp(:,:)'; siz = size(tmp);
        % recombine into tangent space in pages
        new_matrix = permute(cat(3,reshape(sum(reshape(bsxfun(@times,tmp,vec(state.surfNormals')),[3 siz(1)/3 siz(2:end)]),1),[siz(1)/3 siz(2:end)]), ...
            reshape(sum(reshape(bsxfun(@times,tmp,vec(state.surfTangents')),[3 siz(1)/3 siz(2:end)]),1),[siz(1)/3 siz(2:end)]), ...
            reshape(sum(reshape(bsxfun(@times,tmp,vec(state.surfBinormals')),[3 siz(1)/3 siz(2:end)]),1),[siz(1)/3 siz(2:end)])),[2 1 3]);
        % flatten again
        state.leadFieldMatrix = new_matrix(:,:);
    else
        error('Unsupported LeadfieldOrientation specified: %s',leadfieldOrientation);
    end
end

% We also store the indices of the vertices of each ROI (in the full
% source space) in a cell array. This allows us to obtain dipole
% centroids for each ROI
for k=1:length(roiAtlasLabels)
    state.roiVertices{k} = find(indices4Structure(hmObj,roiAtlasLabels{k}))';
    %         [~,state.roiVertices{k}] = removeStructureFromSourceSpace(hmObj,state.brainStructsToRemove{k});
end
% append any additional ROI defined by 'roiVertices' input argument
state.roiVertices = [state.roiVertices roiVertices];

% get ROI indices into reduced source space. These are used for
% integration over current source density (CSD) within each ROI
LFMcols = 1:length(hmObj.atlas.colorTable);
LFMcols(state.rmIndices) = [];
for k=1:length(state.roiVertices)
    state.roiVerticesReduced{k} = find(ismember(LFMcols,state.roiVertices{k}));
end

% set up the labeling for any vertex-index-defined ROIs
if isempty(roiVerticesLabels) && ~isempty(roiVertices)
    % use generic labeling scheme {'ROI1' 'ROI2' ... 'ROIN'}
    roiVerticesLabels = cellstr(num2str((1:length(roiVertices))'))';
    roiVerticesLabels = cellfun(@(x)['ROI' x],roiVerticesLabels,'UniformOutput',false);
end
state.roiLabels = [roiAtlasLabels roiVerticesLabels];

% if needed, combine ROIs according to rules in combineROIs
if ~isempty(combineROIs)
    if mod(length(combineROIs),2)
        error('BCILAB:BadNameValuePair','CombineROI must be a cell array of name,value pairs');
    end
    newRoiLabels = combineROIs(1:2:end-1);
    roiSubLabels = combineROIs(2:2:end);
    numroi = length(newRoiLabels);
    newRoiVerticesReduced = cell(1,numroi);
    newRoiVertices        = cell(1,numroi);
    roiIdxInAtlas  = cell(1,numroi);
    for k=1:numroi
        matchedLabels = ismember(roiSubLabels{k},roiAtlasLabels);
        if ~all(matchedLabels)
            error('BCILAB:BadROINames', ...
                  'The following ROIs in ''CombineROIs'' are not in the ROIAtlasLabels list: %s', ...
                          hlp_tostring(roiSubLabels{k}(~matchedLabels)));
        end
        roiIdxInAtlas{k}  = find(ismember(roiAtlasLabels,roiSubLabels{k}));
        A = cell2mat(state.roiVerticesReduced(roiIdxInAtlas{k}));
        [tmp, idx] = unique(A);
        newRoiVerticesReduced{k} = A(sort(idx));
        % do the same for roiVertices
        A = cell2mat(state.roiVertices(roiIdxInAtlas{k}));
        [tmp, idx] = unique(A);
        newRoiVertices{k} = A(sort(idx));
    end

    switch combineMode
        case 'merge'
            % remove all merged ROIs ...
            mergedROI = unique(cell2mat(roiIdxInAtlas));
            state.roiVerticesReduced(mergedROI) = [];
            state.roiVertices(mergedROI)        = [];
            state.roiLabels(mergedROI)          = [];

            % ... and append the new ROIs
            state.roiVerticesReduced = [state.roiVerticesReduced newRoiVerticesReduced];
            state.roiVertices        = [state.roiVertices newRoiVertices];
            state.roiLabels          = [state.roiLabels newRoiLabels];
        case 'replace'
            % keep only the new rois
            state.roiVerticesReduced = newRoiVerticesReduced;
            state.roiVertices        = newRoiVertices;
            state.roiLabels          = newRoiLabels;
        case 'append'
            % append the new ROIs
            state.roiVerticesReduced = [state.roiVerticesReduced newRoiVerticesReduced];
            state.roiVertices        = [state.roiVertices newRoiVertices];
            state.roiLabels          = [state.roiLabels newRoiLabels];
    end
end

% prune the lead field matrix to contain only desired channels
hmChanlabels  = lower(hmObj.channelLabel);
if isempty(channels)
    channels = lower({chanlocs.labels});
end

% use only selected channels that are in the head model
% chaninds is also an integer permutation vector so that rows of
% LFM match the ordering of 'channels' variable
[~, ia, ib] = intersect(lower(hmChanlabels),lower(channels),'stable');
if max(ib)>length(ia)
    error('The following channels could not be matched to the headmodel: %s', ...
          hlp_tostring(channels(ib>length(ia))));
end
chaninds = ia; %ia(ib);
% if nnz(chaninds)~=length(channels)
%     error('Some channels could not be matched to the headmodel');
% end
state.leadFieldMatrix = state.leadFieldMatrix(chaninds,:,:);

if usecar
    % apply common-average rereference
    state.leadFieldMatrix = bsxfun(@minus,state.leadFieldMatrix,mean(state.leadFieldMatrix));
end

% load the original source space (non-reduced)
tmp = load(hmObj.surfaces);
fn  = fieldnames(tmp);
state.sourceSpace = tmp.(fn{1})(end);  % dim end = cortical surface

% make dipfit structure containing centroids (dipfit.model.posxyz) and
% surface mesh (dipfit.model.surfmesh) for each ROI as well as complete
% surface mesh (dipfit.surfmesh). We also store the indices of each ROI
% into the complete surface mesh (dipfit.model.meshVertices)
if makeDipfitStruct
    state.dipfit = hlp_makeDipfitStruct(state.sourceSpace,state.roiVertices,reducedSpace);
end

% initialize the state of the solver
state.solverState = [];

repo: REST


text:
classdef test_Correlation < matlab.unittest.TestCase
   %UNTITLED Summary of this class goes here
   %   Detailed explanation goes here

   properties
   end

   methods(Test)

      % Cases for adjusttrig
      function TestAdjustTrigDefaults(testCase)
         c = correlation.demo();
         c = c.adjustrig();
      end
      function TestAdjustTrigTimeshift(testCase)
         c = correlation.demo();
         c = c.adjustrig('index', 10);
      end
      function TestAdjustTrigMin(testCase)
         c = correlation.demo();
         c = c.adjustrig('min');
      end
      function TestAdjustTrigMedian(testCase)
         c = correlation.demo();
         c = c.adjustrig('median');
      end
      function TestAdjustTrigMaxLag(testCase)
         c = correlation.demo();
         c = c.adjustrig('min', 1);
      end
      function TestAdjustTrigIndex(testCase)
         c = correlation.demo();
         c = c.adjustrig('index');
      end
      function TestAdjustTrigIndexRelativeToSpecificTrace(testCase)
         c = correlation.demo();
         c = c.adjustrig('index', 10);
      end
      function TestAdjustTrigLeastSquares(testCase)
         c = correlation.demo();
         c = c.adjustrig('lsq');

         % TEST RESULT
      end

      % tests for agc
      function TestAutoGainControl(testCase)
      end
      function TestAlign(testCase)
      end
      function testButter(testCase)
      end
      function testCat(testCase)
      end
      function testCheck(testCase)
      end
      function testCluster(testCase)
      end
      function testColormap(testCase)
      end
      function testConv(testCase)
      end
      function testCrop(testCase)
      end
      function testDeconv(testCase)
      end
      function testDemean(testCase)
      end
      function testDetrend(testCase)
      end
      function testDiff(testCase)
      end
      function testFind(testCae)
      end

   end
end

author: geoscience-community-codes

repo: GISMO


text:
% VL_IHASHFIND  Find labels in an integer hash table
%   SEL = VL_IHASHFIND(ID, NEXT, K, X) returns a vector SEL of the
%   entires in the hash table ID,NEXT,K corresponding to the labels
%   stored as columns of X.
%
%   The format is the same as for the functino VL_IHASHSUM().
%
%   See also: VL_IHASHSUM().

% Author: Andrea Vedaldi

% Copyright (C) 2008-12 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function y=logsum(x,d,k)
%LOGSUM logsum(x,d,k)=log(sum(k.*exp(x),d))
%
% Usage: (1) y=logsum(x) % log(sum(exp(x)))
%        (2) f=0.1*log(10); y=logsm(x*f)/f;  % add powers in dB
%
% Inputs:  X(M,N,...) data matrix to sum
%          D          optional dimension to sum along [1st non-singular dimension]
%          K(M,N,...) optional scaling matrix. It must either be idential
%                     in size to X, or else be a vector whose length is
%                     equal to the size of dimension D of X
%
% Outputs: Y(1,N,...) = log(sum(exp(X).*K,D))
%
% This routine evaluates the given expression for Y but takes care to avoid
% overflow or underflow.

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

if nargin==1 || ~numel(d)
    d=[find(size(x)-1) 1];
    d=d(1);
end
n=size(x,d);
if n<=1,            % use efficient computation if only one term in the sum
    if nargin<3
        y=x;
    else
        y=x+log(k);
    end
    return;
end
s=size(x);
p=[d:ndims(x) 1:d-1];
z=reshape(permute(x,p),n,prod(s)/n);
q=max(z,[],1);              % we subtract y from each row to avoid dynamic range problems
a=(q==Inf)|(q==-Inf);       % check for infinities
if nargin<3
    y=q+log(sum(exp(z-q(ones(n,1),:)),1));
elseif numel(k)==n
    y=q+log(sum(exp(z-q(ones(n,1),:)).*repmat(k(:),1,prod(s)/n),1));
else
    y=q+log(sum(exp(z-q(ones(n,1),:)).*reshape(permute(k,p),n,prod(s)/n),1));
end
y(a)=q(a);                  % correct any column whose max is +-Inf
s(d)=1;                     % update the dimension of the summed component
y=ipermute(reshape(y,s(p)),p);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function LU = extract_bounds_from_norm_operator(LU,extstruct,extvariables,i);
arg = extstruct(i).arg{1};
if min(size(arg))==1 & length(extstruct(i).arg)>1
    p = extstruct(i).arg{2};
    % Maximize each element conservatively, and use the norm of
    % that worst-case vector
    U = max(abs(LU(getvariables(arg),:)),[],2);
    if ~all(isinf(U))
        xmax = abs(getbase(arg))*[1;U];
        LU(extvariables(i),2) = min([norm(xmax,p) LU(extvariables(i),2)]);
    end
end
% norms are positive...
LU(extvariables(i),1) = max([0 LU(extvariables(i),1)]);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ know, x ] = p28_sol ( n )

%*****************************************************************************80
%
%% P28_SOL returns the solution for problem 28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = zeros ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ X, Y ] = plot_fun ( )

%*****************************************************************************80
%
%% PLOT_FUN demonstrates MATLAB's SPMD command for parallel programming.
%
%  Discussion:
%
%    Each worker computes X and Y data for a portion of a sine curve.
%
%    The client patches the composite data into numeric arrays and
%    returns that as the function output.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real X(N,100), Y(N,100), the data, stored as N vectors.
%

%
%  Set a list of colors that will be cycled through with each
%  plot command.
%
  spmd
    a = 2 * pi * ( labindex - 1 ) / numlabs;
    b = 2 * pi *   labindex       / numlabs;
    x = linspace ( a, b, 100 );
    y = sin ( x );
    fprintf ( 1, '  Lab %d works on [%f,%f].\n', labindex, a, b );
  end
%
%  Copy the composite data from the workers, and
%  convert it into numeric arrays.
%
%  (Only client data can be returned as output arguments!)
%
  n = matlabpool ( 'size' );

  X = [];
  Y = [];

  for i = 1 : n
    X = [ X; x{i} ];
    Y = [ Y; y{i} ];
  end

  return
end


text:
function [ state ] =movePTPArcXY_AC(t,theta,c,vel)
%% This function is used for moving the endeffector on an arc in the XY plane, for the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] =movePTPArcXY_AC(t,theta,c,vel)

%% About:
% This function is used to move the end-effector on an arc in the XY plane,

%% Arreguments:
% t: is the TCP/IP connection
% theta: is the arc angle, in radians
% c: the XY coordinates of the center of the circle, it is 1x2 vector.
% vel : is a double, defines the motion velocity mm/sec.

% Copyright, Mohammad SAFEEA, 9th of May 2017

    k=[0;0;1];
    pos=getEEFPos( t );
    c=colVec(c);
    c1=[c;pos{3}];
    state=movePTPArc_AC(t,theta,c1,k,vel);

end

function [ y ] = colVec( v)
% Convert a vector to a column vector:
    if(size(v,2)==1)
        y=v;
    else
        y=v';
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [rpav]=realized_preaveraged_bipower_variation(price,time,timeType,samplingType,samplingInterval,options)
% Estimated integrated variance using Preaveraged Bipower Variation
%
% USAGE:
%   [RPBV] = realized_preaveraged_bipower_variation(PRICE)
%   [RPBV] = realized_preaveraged_bipower_variation(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL)
%   [RPBV] = realized_preaveraged_bipower_variation(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,OPTIONS)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - [OPTIONAL] m by 1 vector of times where TIME(i) corresponds to PRICE(i).
%   TIMETYPE         - [OPTIONAL] String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS.mmm, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - [OPTIONAL] String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVAL ticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL  - [OPTIONAL] Scalar integer or n by 1 vector whose meaning depends on the
%                         selected SAMPLINGTYPE
%   OPTIONS           - [OPTIONAL] Preaveraged Realized Variance option structure initialized by calling
%                         realized_options('Preaveraging'). See help realized_options for a description of
%                         available options.
%
% OUTPUTS:
%   RPBV - Preaveraged bipower variation estimate
%
% COMMENTS:
%  Follows Christensen, Oomen and Podolski (2014) most closely, with the
%  noise variance estimator used in Hautsch and Podolski (2013)
%
% EXAMPLES:
%
%  See also PREALIZED_PREAVERAGED_VARIANCE, REALIZED_OPTIONS, REALIZED_KERNEL,
%  REALIZED_NOISE_ESTIMATE, REALIZED_VARIANCE, REALIZED_VARIANCE_OPTIMAL_SAMPLING,
%  REALIZED_RANGE, REALIZED_QUANTILE_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 2/27/2014

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 1
        m = length(price);
        time = linspace(9.5*3600,16*3600,m)';
        timeType = 'seconds';
        samplingType = 'businesstime';
        samplingInterval = 1;
        options = realized_options('preaveraging');
    case 5
        options = realized_options('preaveraging');
    case 6
        % Nothing
    otherwise
        error('One, five or six inputs required.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need method to estimate omega, and options.

% Inserted to protect against inputing integer times
time = double(time);

% 1. Filter the rpice
filteredPrice = realized_price_filter(price,time,timeType,samplingType,samplingInterval);
returns = diff(log(filteredPrice));
% 2. Compute rpreaveraged returns
theta = options.theta;
m = length(returns);
K = ceil(theta * sqrt(m));
g = @(x) min(x,1-x);
w = g((1:(K-1))/K);

preav_returns = nan(m-K+2,1);
for i=1:m-K+2
    preav_returns(i) = w*returns(i:i+K-2);
end

% 3. Compute constants
psi_1 = K*sum((g((1:K)/K) - g((0:(K-1))/K)).^2);
psi_2 = 1/K*sum((g((1:(K-1))/K)).^2);
% 4. Estimate noise
[noiseVariance, ~, ~, noiseEstimateOomen] = realized_noise_estimate(price, time, timeType, options);

omega = noiseEstimateOomen;
if omega<0
    omega = noiseVariance;
end

% 5. Compute Bipower Variation
mu_1 = sqrt(2/pi);
const1 = m/(m-2*K+2);
const2 = 1/(K*psi_2*mu_1^2);

bias = psi_1/(theta.^2 * psi_2) * omega^2;

lead = preav_returns(K:end);
lag = preav_returns(1:end-K+1);

rpav = const1 * const2 * sum(abs(lead.*lag)) - bias;

author: bashtage

repo: mfe-toolbox


text:
%PREX-MCPLOT   PRTools example on a multi-class classifier plot
help prex_mcplot
echo on
gridsize(100)
            % Generate twice normally distributed 2-class data in 2D
a = +gendath([20,20]);    % data only
b = +gendath([20,20]);    % data only
            % Shift the second data by a vector [5,5]
            % and combine it with the first dataset in to A
A = [a; b+5];
            % Generate 4-class labels
lab = genlab([20 20 20 20],[1 2 3 4]');
            % Construct a 4-class dataset A
A = prdataset(A,lab);
A = setname(A,'4-class dataset')
            % Plot this 4-class dataset
figure
            % Make a scatter-plot of the right size
scatterd(A,'.'); drawnow;
            % Compute normal densities based quadratic classifier
w = qdc(A);
            % Plot filled classification regions
plotc(w,'col'); drawnow;
hold on;
            % Redraw the scatter-plot
scatterd(A);
hold off
echo off

author: marianux

repo: ecg-kit

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


text:
clear; close all;
%% settings
folder = '291-aug/';  % 291 augment dataset
savepath = 'train_hdf5/train_x2.h5';  % save filename
size_input = 29; % 29 | 39
size_label = 57; % 57 | 77
scale = 2; % upsacling factor
stride = 29; % subimages sampling stride

%% initialization
data = zeros(size_input, size_input, 1, 1);
label = zeros(size_label, size_label, 1, 1);
bic = zeros(size_label, size_label, 1, 1);
count_input = 0;
count_label = 0;
count_bic = 0;

%% generate data
filepaths=dir(fullfile(folder,'*.bmp'));

for i = 1 :length(filepaths)

    image = imread(fullfile(folder, filepaths(i).name));
    if size(image,3) == 1
        continue;
    end

    image = rgb2ycbcr(image);
    image = im2double(image(:, :, 1)); % uint8 to double, ranges from [16/255, 235/255]

    im_label = modcrop(image, scale); % high resolution subimage
    im_input = imresize(im_label, 1/scale, 'bicubic'); % low resolution subimage
    im_bic=imresize(im_input,scale,'bicubic'); % interpolated low resolution subimage

    [hei,wid] = size(im_input); % LR subimage size
    [hei_label,wid_label]=size(im_label);% HR subimage size
    stride_label=stride*scale;
    % two-step cut
    % data
    for x = 1 : stride : hei - size_input + 1
        for y = 1 : stride : wid - size_input + 1
            subim_input = im_input(x : size_input + x - 1, y : size_input + y - 1);
            count_input = count_input + 1;
            data(:, :, 1, count_input) = subim_input;
        end

    end

    % label
    for x = 1 : stride_label : hei_label -size_label + 1
        for y = 1 : stride_label : wid_label -size_label +1
            locx=x+scale-1;
            locy=y+scale-1;
            subim_label=im_label(locx:size_label+locx-1,locy:size_label+locy-1);
            count_label=count_label+1;
            label(:, :, 1, count_label)=subim_label;
        end
    end

    % bicubic
    for x = 1 : stride_label : hei_label -size_label + 1
        for y = 1 : stride_label : wid_label -size_label +1
            locx=x+scale-1;
            locy=y+scale-1;
            subim_bic=im_bic(locx:size_label+locx-1,locy:size_label+locy-1);
            count_bic=count_bic+1;
            bic(:,:,1,count_bic)=subim_bic;
        end
    end

end

assert(count_label==count_input, 'Number of samples should be matched between data and labels');
count=count_label;

order = randperm(count);
data = data(:, :, 1, order);
label = label(:, :, 1, order);
bic = bic(:, :, 1, order);

%% writing to HDF5
chunksz = 64; % chunksize
created_flag = false;
totalct = 0;

for batchno = 1:floor(count/chunksz)
    last_read = (batchno-1)*chunksz;
    batchdata = data(:,:,1,last_read+1:last_read+chunksz);
    batchlabs = label(:,:,1,last_read+1:last_read+chunksz);
    batchbic = bic(:,:,1,last_read+1:last_read+chunksz);

    startloc = struct('dat',[1,1,1,totalct+1], 'lab', [1,1,1,totalct+1], 'bic', [1,1,1,totalct+1]);
    curr_dat_sz = store2hdf5(savepath, batchdata, batchlabs, batchbic, ~created_flag, startloc, chunksz);
    created_flag = true;
    totalct = curr_dat_sz(end);
end
h5disp(savepath);

repo: IDN-Caffe


text:
function sFileOut = out_fopen_bst(OutputFile, sFileIn, ChannelMat, EpochSize)
% OUT_FOPEN_BST: Saves the header of a new empty Brainstorm binary file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2014-2019; Martin Cousineau, 2018

% Get file comment
[fPath, fBase, fExt] = bst_fileparts(OutputFile);

% Create a new header structure
sFileOut = sFileIn;
sFileOut.filename  = OutputFile;
sFileOut.condition = '';
sFileOut.format    = 'BST-BIN';
sFileOut.byteorder = 'l';
sFileOut.comment   = fBase;
sFileOut.header    = struct();
sFileOut.header.device    = sFileOut.device;
sFileOut.header.sfreq     = sFileOut.prop.sfreq;
sFileOut.header.starttime = sFileOut.prop.times(1);
sFileOut.header.navg      = sFileOut.prop.nAvg;
% sFileOut.header.version   = 51;   % April 2019
sFileOut.header.version   = 52;   % March 2023
sFileOut.header.nsamples  = round((sFileOut.prop.times(2) - sFileOut.prop.times(1)) .* sFileOut.prop.sfreq) + 1;
sFileOut.header.epochsize = EpochSize;
sFileOut.header.nchannels = length(ChannelMat.Channel);
% Force the destination compensation level
sFileOut.prop.currCtfComp = sFileOut.prop.destCtfComp;
sFileOut.header.ctfcomp   = sFileOut.prop.destCtfComp;
% Open file
fid = fopen(OutputFile, 'w+', sFileOut.byteorder);
if (fid == -1)
    error(['Could not open output file: "' OutputFile '".']);
end

% ===== FORMAT HEADER =====
fwrite(fid, 'BSTBIN', 'char');                               % CHAR(6)    : Format
fwrite(fid, sFileOut.header.version, 'uint8');               % UINT8(1)   : Version of the format
fwrite(fid, str_zeros(sFileOut.header.device, 40), 'char');  % CHAR(40)   : Device used for recording
fwrite(fid, sFileOut.header.sfreq, 'float32');               % UINT32(1)  : Sampling frequency
fwrite(fid, sFileOut.header.starttime, 'float32');           % FLOAT32(1) : Start time
fwrite(fid, sFileOut.header.navg, 'uint32');                 % UINT32(1)  : Number of files averaged
if ~isempty(sFileOut.header.ctfcomp)
    fwrite(fid, sFileOut.header.ctfcomp, 'uint8');           % UINT8(1)   : CTF compensation status (0,1,2,3)
else
    fwrite(fid, 0, 'uint8');
end
fwrite(fid, sFileOut.header.nsamples, 'uint32');             % UINT32(1)  : Total number of samples
fwrite(fid, sFileOut.header.epochsize, 'uint32');            % UINT32(1)  : Number of samples per epoch
fwrite(fid, sFileOut.header.nchannels, 'uint32');            % UINT32(1)  : Number of channels

% ===== CHANNEL LOCATIONS =====
fwrite(fid, str_zeros(ChannelMat.Comment, 40), 'char');                 % CHAR(40)   : Channel file comment
for i = 1:length(ChannelMat.Channel)
    fwrite(fid, str_zeros(ChannelMat.Channel(i).Name, 20), 'char');     % CHAR(20)   : Channel name
    fwrite(fid, str_zeros(ChannelMat.Channel(i).Type, 20), 'char');     % CHAR(20)   : Channel type
    fwrite(fid, str_zeros(ChannelMat.Channel(i).Comment, 40), 'char');  % CHAR(40)   : Channel comment
    % Loc
    fwrite(fid, size(ChannelMat.Channel(i).Loc, 2), 'uint8');           % UINT8(1)   : Number of location points
    if ~isempty(ChannelMat.Channel(i).Loc)
        fwrite(fid, ChannelMat.Channel(i).Loc, 'float32');              % FLOAT32(N) : Locations [3xN]
    end
    % Orient
    fwrite(fid, size(ChannelMat.Channel(i).Orient, 2), 'uint8');        % UINT8(1)   : Number of orientations points
    if ~isempty(ChannelMat.Channel(i).Orient)
        fwrite(fid, ChannelMat.Channel(i).Orient, 'float32');           % FLOAT32(N) : Orientations [3xN]
    end
    % Weight
    fwrite(fid, length(ChannelMat.Channel(i).Weight), 'uint8');         % UINT8(1)   : Number of weights
    if ~isempty(ChannelMat.Channel(i).Weight)
        fwrite(fid, ChannelMat.Channel(i).Weight, 'float32');           % FLOAT32(N) : Weights [1xN]
    end
    % Channel flag
    fwrite(fid, sFileIn.channelflag(i), 'int8');                        % INT8(1)    : Channel flag
end

% ===== CTF =====
% CTF compensation matrix
fwrite(fid, size(ChannelMat.MegRefCoef), 'uint32');                     % UINT32(2)  : Dimensions of the MegRefCoef matrix
fwrite(fid, ChannelMat.MegRefCoef, 'float32');                          % FLOAT32(N) : MegRefCoef matrix

% ===== SSP PROJECTORS =====
fwrite(fid, length(ChannelMat.Projector), 'uint32');                      % UINT32(1)  : Number of projectors
for i = 1:length(ChannelMat.Projector)
    fwrite(fid, str_zeros(ChannelMat.Projector(i).Comment, 40), 'char');  % CHAR(40)   : Projector comment
    fwrite(fid, size(ChannelMat.Projector(i).Components), 'uint32');      % UINT32(2)  : Dimensions of the Components matrix
    fwrite(fid, ChannelMat.Projector(i).Components, 'float32');           % FLOAT32(N) : Components matrix
    fwrite(fid, size(ChannelMat.Projector(i).CompMask), 'uint32');        % UINT32(2)  : Dimensions of the CompMask matrix
    fwrite(fid, ChannelMat.Projector(i).CompMask, 'float32');             % FLOAT32(N) : CompMask matrix
    fwrite(fid, size(ChannelMat.Projector(i).SingVal), 'uint32');         % UINT32(2)  : Dimensions of the SingVal matrix
    if ~isempty(ChannelMat.Projector(i).SingVal)
        fwrite(fid, ChannelMat.Projector(i).SingVal, 'float32');          % FLOAT32(N) : SingVal matrix
    end
    fwrite(fid, ChannelMat.Projector(i).Status, 'int8');                % INT8(1)    : Status
end

% ===== HEAD POINTS =====
if ~isempty(ChannelMat.HeadPoints)
    fwrite(fid, length(ChannelMat.HeadPoints.Label), 'uint32');              % UINT32(1)  : Number of head points
    for i = 1:length(ChannelMat.HeadPoints.Label)
        fwrite(fid, ChannelMat.HeadPoints.Loc(:,i), 'float32');              % FLOAT32(3) : (X,Y,Z) positions in meters
        fwrite(fid, str_zeros(ChannelMat.HeadPoints.Label{i}, 10), 'char');  % CHAR(10)   : Point label
        fwrite(fid, str_zeros(ChannelMat.HeadPoints.Type{i}, 10), 'char');   % CHAR(10)   : Point type
    end
else
    fwrite(fid, 0, 'uint32');
end

% ===== FIDUCIALS =====
% Write fiducials positions (SCS)
if isfield(ChannelMat, 'SCS') && ~isempty(ChannelMat.SCS) ...
&& isfield(ChannelMat.SCS, 'NAS') && (length(ChannelMat.SCS.NAS) == 3) ...
&& isfield(ChannelMat.SCS, 'LPA') && (length(ChannelMat.SCS.LPA) == 3) ...
&& isfield(ChannelMat.SCS, 'RPA') && (length(ChannelMat.SCS.RPA) == 3) ...
&& isfield(ChannelMat.SCS, 'R')   && isequal(size(ChannelMat.SCS.R), [3,3]) ...
&& isfield(ChannelMat.SCS, 'T')   && (length(ChannelMat.SCS.T) == 3) ...
&& isfield(ChannelMat.SCS, 'Origin') && (length(ChannelMat.SCS.Origin) == 3)
    fwrite(fid, 1, 'int8');                                                  % INT8(1)    : Are the fiducials saved in the file?
    fwrite(fid, ChannelMat.SCS.NAS, 'float32');                              % FLOAT32(3) : NAS (X,Y,Z) positions in meters
    fwrite(fid, ChannelMat.SCS.LPA, 'float32');                              % FLOAT32(3) : LPA (X,Y,Z) positions in meters
    fwrite(fid, ChannelMat.SCS.RPA, 'float32');                              % FLOAT32(3) : RPA (X,Y,Z) positions in meters
    fwrite(fid, ChannelMat.SCS.R, 'float32');                                % FLOAT32(9) : Rotation matrix [3x3]
    fwrite(fid, ChannelMat.SCS.T, 'float32');                                % FLOAT32(3) : Translation vector [3x1]
    fwrite(fid, ChannelMat.SCS.Origin, 'float32');                           % FLOAT32(3) : Origin (X,Y,Z) positions in meters
else
    fwrite(fid, 0, 'int8');                                                  % INT8(1)    : Are the fiducials saved in the file?
end

% ===== EVENTS =====
fwrite(fid, length(sFileOut.events), 'uint32');                              % UINT32(1)  : Number of event categories
for iEvt = 1:length(sFileOut.events)
    isExtended = (size(sFileOut.events(iEvt).times,1) == 2);
    evtLabel = sFileOut.events(iEvt).label;
    labelLength = length(evtLabel);
    nOcc = size(sFileOut.events(iEvt).times,2);
    fwrite(fid, labelLength, 'uint8');                                       % UINT8(1)   : Length of event name (1 to 255)
    fwrite(fid, str_zeros(evtLabel, labelLength), 'char');                   % CHAR(??)   : Event name
    fwrite(fid, sFileOut.events(iEvt).color, 'float32');                     % FLOAT32(3) : Event color
    fwrite(fid, isExtended, 'int8');                                         % INT8(1)    : Event type (0=regular, 1=extended)
    fwrite(fid, nOcc, 'uint32');                                             % UINT32(1)  : Number of occurrences
    % If there are event occurrences
    if ~isempty(sFileOut.events(iEvt).times)
        % Write latencies
        fwrite(fid, sFileOut.events(iEvt).times, 'float32');                 % FLOAT32(2*N) : Time in seconds
        % Write list of channels associated to each event
        isChannels = ~isempty(sFileOut.events(iEvt).channels);
        fwrite(fid, uint8(isChannels), 'uint8');                             % UINT8(1) : 1 if channels list is present, 0 otherwise
        if isChannels
            for iOcc = 1:nOcc
                nChannels = length(sFileOut.events(iEvt).channels{iOcc});
                fwrite(fid, nChannels, 'uint16');                            % UINT16(1) : Number of channels associated to this event
                for iChan = 1:nChannels
                    chLabel = sFileOut.events(iEvt).channels{iOcc}{iChan};
                    labelLength = length(chLabel);
                    fwrite(fid, labelLength, 'uint8');                       % UINT8(1) : Length of channel name (1 to 255)
                    fwrite(fid, str_zeros(chLabel, labelLength), 'char');    % CHAR(??) : Channel name
                end
            end
        end
        % Read list of notes associated to each event
        isNotes = ~isempty(sFileOut.events(iEvt).notes);
        fwrite(fid, uint8(isNotes), 'uint8');                                % UINT8(1) : 1 if notes list is present, 0 otherwise
        if isNotes
            for iOcc = 1:nOcc
                noteLabel = sFileOut.events(iEvt).notes{iOcc};
                labelLength = length(noteLabel);
                fwrite(fid, labelLength, 'uint16');                          % UINT16(1) : Length of note text
                fwrite(fid, str_zeros(noteLabel, labelLength), 'char');      % CHAR(??) : Note text
            end
        end
    end
end

% Save total header size (to allow fast skipping when reading)
sFileOut.header.hdrsize = ftell(fid);
% Close file
fclose(fid);

end

%% ===== HELPER FUNCTIONS =====
function sout = str_zeros(sin, N)
    sout = char(zeros(1,N));
    if (length(sin) <= N)
        sout(1:length(sin)) = sin;
    else
        sout = sin(1:N);
    end
end

author: brainstorm-tools


text:
function pathname = filesepStandard(pathname0, options)
%
% Usage:
%    pathname = filesepStandard(pathname, options)
%
% Takes a pathname as argument and replaces any non-standard file/folder
% separators with standard ones, that is '/'. It also gets rid of redundant
% file seps
%
% Example:
%
%   >> pathname = 'C:\dir1\\\dir2\\dir3\test1/\test2/'
%   >> pathname = filesepStandard(pathname)
%
%   pathname =
%
%   C:/dir1/dir2/dir3/test1/test2/
%
%

pathname = '';

if nargin==0
    return
end
if isempty(pathname0)
    return
end
if ~ischar(pathname0) && ~iscell(pathname0)
    return
end

if ~exist('options', 'var')
    options = '';
end

if ischar(pathname0)
    pathname0 = removeExtraDots(pathname0);

    % Do basic error check to see if path exists; if not wexist without
    % doing anything
    if ~optionExists(options,'nameonly')
        if ~ispathvalid(pathname0)
            return
        end
    end

    % Change all path file separators to standard forward slash
    idxs = [];
    k = find(pathname0=='\' | pathname0=='/');
    for ii = 1:length(k)
        if (ii>1) && (k(ii) == k(ii-1)+1)
            idxs = [idxs, k(ii)]; %#ok<AGROW>
            continue;
        end
        pathname0(k(ii)) = '/';
    end

    % Remove any extraneous file separators
    pathname0(idxs) = '';

    % Change path to full path if option requesting it exists
    if optionExists(options,'full') || optionExists(options,'fullpath') || optionExists(options,'absolute')
        if ispathvalid(pathname0)
            pathname0 = fullpath(pathname0);
        end
    end

    % Add traling separator only for directory path names
    if (isdir_private(pathname0) || optionExists(options, 'dir')) && ~optionExists(options, 'file')
        if pathname0(end) ~= '/'
            pathname0(end+1) = '/';
        end
    elseif (isfile_private(pathname0) || optionExists(options, 'file')) && ~optionExists(options, 'dir')
        if pathname0(end) == '/'
            pathname0(end) = '';
        end
    else
        if pathname0(end) == '/'
            pathname0(end) = '';
        end
    end

    if optionExists(options, 'filesepwide')
        pathname0 = filesepwide(pathname0);
    end

    % Change path to full path if option requesting it exists
    pathname = pathname0;
    return;
end
pathname = pathname0;

% Call filesepStandard recursively for all path names in cell array
for ii = 1:length(pathname)
    pathname{ii} = filesepStandard(pathname{ii}, options);
end

% ---------------------------------------------------------------------------
function pname = removeExtraDots(pname)
k = cell(4,3);

% Case 1:
k{1,1} = strfind(pname, '/./');
k{2,1} = strfind(pname, '/.\');
k{3,1} = strfind(pname, '\.\');
k{4,1} = strfind(pname, '\./');
for ii = 1:length(k(:,1))
    for jj = length(k{ii,1}):-1:1
        pname(k{ii,1}(jj)+1:k{ii,1}(jj)+2) = '';
    end
end

% Case 2:
k{1,2} = strfind(pname, '/.');
k{2,2} = strfind(pname, '\.');
for ii = 1:length(k(:,2))
    if ~isempty(k{ii,2})
        if k{ii,2}+1<length(pname)
            continue
        end
        pname(k{ii,2}+1) = '';
    end
end

% --------------------------------------------------------------
function pathname = filesepwide(pathname0)
c = str2cell(pathname0, '/');
pathname = '';
for ii = 1:length(c)
    if isempty(pathname)
        if pathname0(1)=='/'
            pathname = ['/ ', c{ii}];
        else
            pathname = c{ii};
        end
    else
        pathname = sprintf('%s / %s', pathname, c{ii});
    end
end

author: BUNPC


text:
function [ x, seed ] = quasigeometric_sample ( a, b, seed )

%*****************************************************************************80
%
%% QUASIGEOMETRIC_SAMPLE samples the Quasigeometric PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the probability of 0 successes.
%    0.0 <= A <= 1.0.
%
%    Input, real B, the depreciation constant.
%    0.0 <= B < 1.0.
%
%    Input, integer SEED, a seed for the random
%    number generator.
%
%    Output, integer X, a sample of the PDF.
%
%    Output, integer SEED, a seed for the random
%    number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = quasigeometric_cdf_inv ( cdf, a, b );

  return
end


text:
% FORM_FACTOR Computer form factors for radiosity computation
%
% [F,V,BC,N] = form_factor(P,E)
%
% Inputs:
%   P  #P by 2 list of 2D vertices
%   E  #E by 2 list of edge indices into P
% Outputs:
%   F  #E by #E matrix of integrated form factor values
%   V  #E by #E matrix of visibility values
%   BC  #E by 2 list of edge barycenters
%   N  #E by 2 list of length-weighted normals

author: alecjacobson

repo: gptoolbox


text:
function [data] = normalise(data,u,std)

[N]=size(data,2);

data=(data-u*ones(1,N)) ./ (std*ones(1,N));
data=data';

author: alexandrebarachant


text:
function [y,yz,yzi] = autoGen_cst_swingFootHeight(q1,q2,q4,q5,l1,l2,l4,l5,stepLength,stepHeight)
%AUTOGEN_CST_SWINGFOOTHEIGHT
%    [Y,YZ,YZI] = AUTOGEN_CST_SWINGFOOTHEIGHT(Q1,Q2,Q4,Q5,L1,L2,L4,L5,STEPLENGTH,STEPHEIGHT)

%    This function was generated by the Symbolic Math Toolbox version 6.3.
%    25-Oct-2015 18:36:31

t3 = sin(q1);
t4 = l1.*t3;
t7 = sin(q2);
t8 = l2.*t7;
t9 = sin(q4);
t10 = l4.*t9;
t11 = sin(q5);
t12 = l5.*t11;
t2 = t4+t8-t10-t12;
t5 = 1.0./stepLength.^2;
t6 = cos(q1);
t13 = cos(q2);
t14 = cos(q4);
t15 = cos(q5);
y = -l1.*t6-l2.*t13+l4.*t14+l5.*t15-stepHeight.*(t2.^2.*t5-1.0);
if nargout > 1
    yz = [t4-l1.*stepHeight.*t2.*t5.*t6.*2.0;t8-l2.*stepHeight.*t2.*t5.*t13.*2.0;-t10+l4.*stepHeight.*t2.*t5.*t14.*2.0;-t12+l5.*stepHeight.*t2.*t5.*t15.*2.0];
end
if nargout > 2
    yzi = [2.0;3.0;5.0;6.0];
end

author: MatthewPeterKelly

repo: OptimTraj