text:
function [conf, args] = vl_argparse(conf, args, varargin)
% VL_ARGPARSE  Parse list of parameter-value pairs
%   CONF = VL_ARGPARSE(CONF, ARGS) updates the structure CONF based on
%   the specified parameter-value pairs ARGS={PAR1, VAL1, ... PARN,
%   VALN}. The function produces an error if an unknown parameter name
%   is passed in.
%
%   [CONF, ARGS] = VL_ARGPARSE(CONF, ARGS) copies any parameter in
%   ARGS that does not match CONF back to ARGS instead of producing an
%   error.
%
%   Example::
%     The function can be used to parse a list of arguments
%     passed to a MATLAB functions:
%
%       function myFunction(x,y,z,varargin)
%       conf.parameterName = defaultValue ;
%       conf = vl_argparse(conf, varargin)
%
%     If only a subset of the options should be parsed, for example
%     because the other options are interpreted by a subroutine, then
%     use the form
%
%      [conf, varargin] = vl_argparse(conf, varargin)
%
%     that copies back to VARARGIN any unknown parameter.
%
%   See also: VL_OVERRIDE(), VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% 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).

if ~isstruct(conf), error('CONF must be a structure') ; end

if length(varargin) > 0, args = {args, varargin{:}} ; end

remainingArgs = {} ;
names = fieldnames(conf) ;

if mod(length(args),2) == 1
  error('Parameter-value pair expected (missing value?).') ;
end

for ai = 1:2:length(args)
  paramName = args{ai} ;
  if ~ischar(paramName)
    error('The name of the parameter number %d is not a string.', (ai-1)/2+1) ;
  end
  value = args{ai+1} ;
  if isfield(conf,paramName)
    conf.(paramName) = value ;
  else
    % try case-insensitive
    i = find(strcmpi(paramName, names)) ;
    if isempty(i)
      if nargout < 2
        error('Unknown parameter ''%s''.', paramName) ;
      else
        remainingArgs(end+1:end+2) = args(ai:ai+1) ;
      end
    else
      conf.(names{i}) = value ;
    end
  end
end

args = remainingArgs ;

author: AlfredXiangWu


text:
 function G = Gtomo2_strip(sg, ig, varargin)
%function G = Gtomo2_strip(sg, ig, options)
%
% Generate a 2D system matrix for tomographic image reconstruction,
% based on a square pixel basis and strip-integral detector model,
% i.e., a rectangular detector response.
% Works for parallel-beam, arc-fan, and flat-fan geometries.
%
% Closely matches Aspire "system 13" and "system 14".
% The main limitation here is that it uses Matlab "sparse" data type,
% so you will run out of memory for large image or sinogram sizes.
%
% in:
%	sg		sino_geom()
%	ig		image_geom()
%
% options:
%	'strip_width'	if 0, then line integrals (default: sg.d)
%	'single'	if 1, then double(single()) the values, cf aspire
%	'gam_max'	max acceptance angle of detector [degrees] (def: 90)
%
% out:
%	G [nb*na,np]	Gsparse object, where np = sum(ig.mask(:))
%
% Copyright 2005-8-16, Jeff Fessler, The University of Michigan

if nargin == 1 && streq(sg, 'test'), Gtomo2_strip_test, return, end
if nargin < 2, help(mfilename), error(mfilename), end

% defaults
arg.gam_max = 90; % degrees
arg.strip_width = [];
arg.single = false;
arg.chat = false;
arg.sse = false;
arg = vararg_pair(arg, varargin);

if isempty(arg.strip_width), arg.strip_width = sg.d; end

if streq(sg.type, 'par')
  if (arg.sse)

    nthread = 1;
    mask = uint8(ig.mask);

    [A B] = jmh_sse_mex('jmh,sse,matrix',int32(nthread), ...
                        single(sg.orbit_start*pi/180), single(sg.orbit*pi/180), ...
                        int32(sg.nb), single(sg.dr), single(sg.offset_r), ...
                        single(arg.strip_width), int32(ig.nx), single(ig.dx), ...
                        single(ig.offset_x), int32(ig.ny), single(ig.dy), ...
                        single(ig.offset_y), mask,int32(sg.na));

    M = ceil((ig.dx * sqrt(2) + arg.strip_width) / sg.dr);

%     jj = [];

%     % Bad idea
%     for y=1:ig.ny
%       for x=1:ig.nx
%         for a=1:sg.na
%           for m=1:M
%             jj = [jj,(y-1)*ig.nx + x];
%           end
%         end
%       end
%     end
    jj = reshape(repmat(1:ig.nx*ig.ny,sg.na*M,1),[1,ig.nx*ig.ny*sg.na*M]);

    G = sparse(double(B+1),double(jj'),double(A),sg.na*sg.nb,ig.nx*ig.ny,length(A));

  else
      G = Gtomo2_strip_par(sg.nb, sg.na, sg.dr, sg.offset_r, ...
                           sg.orbit, sg.orbit_start, arg.strip_width, ...
                           ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y, ...
                           ig.mask, arg.single);
  end
elseif streq(sg.type, 'fan')
	G = Gtomo2_strip_fan(sg.nb, sg.na, sg.ds, sg.offset_s, ...
		sg.ar, arg.strip_width, ...
		deg2rad(arg.gam_max), ...
		sg.dsd, sg.dso, sg.dfs, ...
		sg.source_offset, ...
		ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y, ...
		ig.mask, arg.single, arg.chat);

else
	error('unknown geometry: %s', sg.type)
end

G = Gsparse(G, 'mask', ig.mask, 'idim', [ig.nx ig.ny], 'odim', [sg.nb sg.na]);

%
% = Gtomo2_strip_par()
%
function G = Gtomo2_strip_par(nb, na, dr, offset_r, ...
	orbit, orbit_start, strip_width, ...
	nx, ny, dx, dy, offset_x, offset_y, mask, is_single);

% pixel centers
wx = (nx-1)/2 + offset_x;
wy = (ny-1)/2 + offset_y;
wb = (nb-1)/2 + offset_r;

x = dx * ([0:nx-1] - wx);
y = dy * ([0:ny-1] - wy); % caution, may not match aspire if offset_y != 0
[x y] = ndgrid(x, y);
x = x(mask(:));
y = y(mask(:));
np = length(x);		% sum(mask(:)) - total # of support pixels

angle = deg2rad(orbit_start + [0:na-1]'/na * orbit);
cang = cos(angle);
sang = sin(angle);
tau = cang * x' + sang * y'; % [na,np] projected pixel center
tau = tau / dr;
tau = tau + wb;

d_max = (abs(cang) + abs(sang)) / 2;
tau_max = (dx * d_max + strip_width/2) / dr;
ib_min = 1 + floor(tau - repmat(tau_max, 1, np));
M = ceil((dx * sqrt(2) + strip_width) / dr);
jj = find(mask(:))'; % all-column G
jj = repmat(jj, na, 1);
jj = col(jj); % so that na=1 case works
list.ii = [];
list.jj = [];
list.ss = [];
for mm=0:M-1
	ticker(mfilename, mm+1, M)
	ib = ib_min + mm;
	good = (ib >= 0) & (ib < nb);
	good = good(:);
	val = square_strip_int(dr * (ib - tau), ...
		repmat(angle, 1, np), 'dx', dx, 'sw', strip_width);
	val = val(:);
	ii = col(ib + repmat([0:na-1]'*nb, 1, np)); % sinogram index
	list.ii = [list.ii; 1+ii(good)];
	list.jj = [list.jj; jj(good)];
	list.ss = [list.ss; val(good)];
end
if is_single
	list.ss = double(single(list.ss)); % stupid matlab insists on double
end
G = sparse(list.ii, list.jj, list.ss, nb*na, nx*ny, length(list.ss));

%
% = Gtomo2_strip_fan()
%
function G = Gtomo2_strip_fan(nb, na, ds, offset_s, ...
	beta, strip_width, ...
	gam_max, ... % maximum angle gamma accepted by detector [radians]
	dsd, dso, dfs, roff, ...
	nx, ny, dx, dy, offset_x, offset_y, mask, is_single, chat);

% pixel centers
wx = (nx-1)/2 + offset_x;
wy = (ny-1)/2 + offset_y;
wb = (nb-1)/2 + offset_s;

x = dx * ([0:nx-1] - wx);
y = dy * ([0:ny-1] - wy); % caution, may not match aspire if offset_y != 0
[x y] = ndgrid(x, y);
x = x(mask(:)); % [np,1]
y = y(mask(:));
np = length(x); % sum(mask(:)) - total # of support pixels
na = length(beta);

cbet = cos(beta);
sbet = sin(beta);

% s0 is "s" value for center of each pixel
gam0 = atan2(cbet * x' + sbet * y' - roff, ...
	dso + sbet * x' - cbet * y'); % [na,np]
clear cbet sbet
if dfs == 0 % 3rd gen (arc)
	s0 = dsd * gam0;
elseif isinf(dfs)
	s0 = dsd * tan(gam0);
else
	error 'unsupported dfs'
end

ang0 = repmat(beta, [1 np]) + gam0;
mag0 = dso * cos(gam0) - roff * sin(gam0) ...
	+ repmat(x', [na 1]) .* sin(ang0) - repmat(y', [na 1]) .* cos(ang0);
gamgood = abs(gam0) < gam_max;
if chat, printm('gam0: %g %g', rad2deg(minmax(gam0(gamgood)))), end
clear gam0
if dfs == 0 % 3rd gen (arc)
	mag0 = mag0 / dsd;
elseif isinf(dfs)
	mag0 = mag0 / dsd ./ (1 + (s0 ./ dsd).^2);
else
	error 'dfs bug'
end
mag0min = min(mag0(gamgood));
if isempty(mag0min)
	G = sparse([], [], [], nb*na, nx*ny, 0);
	return
end

tau = s0 / ds + wb; % [na,np], unitless

M = ceil((dx/mag0min * sqrt(2) + strip_width) / ds); % conservative!
if chat, printm('M=%d, mag0min=%g', M, mag0min), end
if M > 100
	printm('Warn: M=%d too large? probably x-ray source is too oblique', M)
	printm('type "return" and hope for the best, but it may be *slow*')
	printm('recommend using "gam_max" to impose incidence angle constraint')
	keyboard
end

ib_min = 1 + floor(tau - M/2);
jj = find(mask(:))';	% all-column G
jj = repmat(jj, na, 1); % [na,np]
jj = col(jj); % so that na=1 case works
list.ii = [];
list.jj = [];
list.ss = [];
for mm=0:M-1
	ticker(mfilename, mm+1, M)
	ib = ib_min + mm;
	val = mag0 .* square_strip_int(ds * (ib - tau), ...
		ang0, 'dx', dx ./ mag0, ...
		'sw', strip_width);
	val = val(:); % for na=1 case
%	val = square_strip_int(ds * (ib - tau) .* mag0, ...
%		ang0, 'dx', dx, 'sw', strip_width .* mag0); % same!
	ii = col(ib + repmat([0:na-1]'*nb, 1, np)); % sinogram index
%	if chat, printm('%d of %d zeros', sum(val(:)==0), length(val(:))), end
	good = (ib >= 0) & (ib < nb);
	good = good & gamgood;
	good = good(:);
	good = good & (val > 0); % remove extra zeros due to conservative M
	list.ii = [list.ii; 1+ii(good)];
	list.jj = [list.jj; jj(good)];
	list.ss = [list.ss; val(good)];
end
if is_single
	list.ss = double(single(list.ss)); % stupid matlab insists on double
end
G = sparse(list.ii, list.jj, list.ss, nb*na, nx*ny, length(list.ss));

%
% test demo
%
function Gtomo2_strip_test
ig = image_geom('nx', 512, 'ny', 480', 'fov', 500, 'down', 8);
ig.mask = ig.circ > 0;
% todo: 3 cases to test: parallel, arc-fan, flat-fan
sg = sino_geom('par', 'nb', 600, 'na', 480, 'dr', 1.2, 'down', ig.down);
% sg = sino_geom('fan', 'nb', 888, 'na', 984, 'ds', 1.0, 'down', ig.down, ...
% 	'offset_s', 0*0.25, 'source_offset', 0*3.0, ...
% 	'dsd', 949, 'dod', 408, 'dfs', 0);

ell = [20 50 150 150 0 1];
ell = [];
[x ell] = ellipse_im(ig, ell, 'oversample', 2);
ya = ellipse_sino(sg, ell, 'oversample', 4);

G = Gtomo2_strip(sg, ig, 'chat', 0);

yd = G * x;
sino = sg.zeros; sino(sg.nb/2, 10) = 1;
b1 = G' * sino;
bu = G' * sg.ones;
max_percent_diff(min(bu(ig.mask)), max(bu(ig.mask)), 'backproject ones')

if im
	clf, im pl 2 2
	im(1, x, 'test image')
	im(2, ya, 'sinogram ya'), cbar
	im(4, yd, 'sinogram yd'), cbar
	im(3, yd-ya, 'yd-ya'), cbar
	if 0
		im(1, ig.mask, 'support mask')
		im(2, b1, 'backproject 1 ray'), cbar
		im(3, bu, 'backproject ones'), cbar
	end
end

% verify consistency with Gtomo2_wtmex (aspire)
if 1 & has_aspire %& streq(sg.type, 'par')
prompt
	Gw = Gtomo2_wtmex(sg, ig, 'pairs', {'strip_width', sg.d});
	yw = Gw * x;
	ys = G * x;
	max_percent_diff(yw, ys, 'sino Gtomo2_wtmex vs Gtomo2_strip')
	if im
		im pl 2 3, im(1, x), cbar
		im(1, ya, 'ya'), cbar
		im(2, ys, 'ys'), cbar
		im(3, yw, 'yw'), cbar
		im(4, ys-yw, 'ys-yw'), cbar
		xlabelf('%g%%', 100*nrms(ys(:), yw(:)))
		im(5, ys-ya, 'ys-ya'), cbar
		xlabelf('%g%%', 100*nrms(ys(:), ya(:)))
		im(6, yw-ya, 'yw-ya'), cbar
		xlabelf('%g%%', 100*nrms(yw(:), ya(:)))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef TestConvertScaleAbs
    %TestConvertScaleAbs

    methods (Static)
        function test_1
            A = rand([10 10], 'single')*200 - 100;
            B = cv.convertScaleAbs(A, 'Alpha',5, 'Beta',3);
            validateattributes(B, {'uint8'}, {'size',size(A)});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [soln,eqn,info] = StokesCRP0(node,elem,bdFlag,pde,option)
%% STOKESCR Stokes equation: CR elements.
%
%   [u,p] = STOKESPCR(node,elem,bdFlag,pde) use Crouzeix and Raviart
%   nonconforming elements to approximate velocity u and piecewise constant
%   to approximate pressure p, repectively.
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0 in \Omega,
%   with
%       Dirichlet boundary condition        u = g_D  on \Gamma_D,
%       Neumann boundary condition du/dn - np = g_N  on \Gamma_N.
%
%  Created by Ming Wang at July, 2012, with discussion of Lin Zhong.
%
% See also Stokes, StokesP2P1, PoissonCR.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

%% Construct Data Structure
[elem2dof,edge] = dofedge(elem);
NE = size(edge,1); NT = size(elem,1); Nu = NE; Np = NT; N = size(node,1);

t = cputime;
%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix for Laplace operator
A = sparse(Nu,Nu);
for i = 1:3
    for j = i:3
        % local to global index map
        ii = double(elem2dof(:,i));
        jj = double(elem2dof(:,j));
        % local stiffness matrix
        Aij = 4*dot(Dlambda(:,:,i),Dlambda(:,:,j),2).*area;
        if (j==i)
            A = A + sparse(ii,jj,Aij,Nu,Nu);
        else
            A = A + sparse([ii,jj],[jj,ii],[Aij; Aij],Nu,Nu);
        end
    end
end
clear Aij
A = blkdiag(A,A);

%% Assemble matrix for divergence operator
d1 = -2.*Dlambda(:,:,1).*[area,area];
d2 = -2.*Dlambda(:,:,2).*[area,area];
d3 = -2.*Dlambda(:,:,3).*[area,area];
Dx = sparse(repmat((1:Np)',3,1),double(elem2dof(:)),...
            [d1(:,1);d2(:,1);d3(:,1)],Np,Nu);
Dy = sparse(repmat((1:Np)',3,1),double(elem2dof(:)),...
            [d1(:,2);d2(:,2);d3(:,2)],Np,Nu);
B = -[Dx Dy];
clear d1 d2 d3 B1 B2

%% Assemble right hand side
f1 = zeros(Nu,1);
f2 = zeros(Nu,1);
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    mid1 = (node(elem(:,2),:) + node(elem(:,3),:))/2;
    mid2 = (node(elem(:,3),:) + node(elem(:,1),:))/2;
    mid3 = (node(elem(:,1),:) + node(elem(:,2),:))/2;
    ft1 = repmat(area,1,2).*pde.f(mid1)/3;
    ft2 = repmat(area,1,2).*pde.f(mid2)/3;
    ft3 = repmat(area,1,2).*pde.f(mid3)/3;
    f1 = accumarray(elem2dof(:),[ft1(:,1);ft2(:,1);ft3(:,1)],[Nu 1]);
    f2 = accumarray(elem2dof(:),[ft1(:,2);ft2(:,2);ft3(:,2)],[Nu 1]);
end

[AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesCR;

%% Record assembeling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(ufreeDof), return; end
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if length(f)+length(g) <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % Multigrid-type  solver for large size systems
        option.solver = 'asmg';
    end
end
solver = option.solver;

%% Solver
switch solver
    case 'direct'
        t = cputime;
        bigA = [AD, BD'; ...
                BD, sparse(Np,Np)];
        bigF = [f; g];
        bigu = [u; p];
        bigFreeDof = [ufreeDof; 2*Nu+pDof];
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                        u(ufreeDof),p,elem,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                          u,p,node,elem,bdFlag,ufreeDof,option);
end

%% Post-process
if length(pDof) ~= Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(p.*area)/sum(area);
    p = p - c;
end

%% Output
soln = struct('u',u,'p',p);
eqn = struct('A',AD,'B',BD,'Lap',A,'f',f,'g',g,...
             'edge',edge,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesCR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesCR
        %% Boundary condition of Stokes equation: CR elements

        %% Initial set up
        f = [f1; f2];    % set in Neumann boundary condition
        g = zeros(Np,1);
        u = zeros(2*Nu,1);
        p = zeros(Np,1);
        ufreeDof = (1:Nu)';
        pDof = (1:Np)';

        if ~exist('bdFlag','var'), bdFlag = []; end
        if ~isfield(pde,'g_D'), pde.g_D = []; end
        if ~isfield(pde,'g_N'), pde.g_N = []; end
        if ~isfield(pde,'g_R'), pde.g_R = []; end

        %% Part 1: Find Dirichlet dof and modify the matrix
        % Find Dirichlet boundary dof: fixedDof and pDof
        isFixedDof = false(Nu,1);
        if ~isempty(bdFlag)       % case: bdFlag is not empty
            isDirichlet(elem2dof(bdFlag(:)==1)) = true;
            isFixedDof(isDirichlet) = true;% dof on D-edges
            fixedDof = find(isFixedDof);
            ufreeDof = find(~isFixedDof);
        end
        if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
            s = accumarray(elem2dof(:), 1, [Nu 1]);
            isFixedDof = (s==1);
            fixedDof = find(isFixedDof);
            ufreeDof = find(~isFixedDof);
        end
        if isempty(fixedDof) % pure Neumann boundary condition
            % pde.g_N could be empty which is homogenous Neumann boundary condition
            fixedDof = 1;
            ufreeDof = 2:Nu;    % eliminate the kernel by enforcing u(1) = 0;
        end

        % Modify the matrix
        % Build Dirichlet boundary condition into the matrix AD by enforcing
        % AD(fixedDof,fixedDof)=I, AD(fixedDof,ufreeDof)=0, AD(ufreeDof,fixedDof)=0.
        % BD(:,fixedDof) = 0 and thus BD'(fixedDof,:) = 0.
        bdidx = zeros(2*Nu,1);
        bdidx(fixedDof) = 1;
        bdidx(Nu+fixedDof) = 1;
        Tbd = spdiags(bdidx,0,2*Nu,2*Nu);
        T = spdiags(1-bdidx,0,2*Nu,2*Nu);
        AD = T*A*T + Tbd;
        BD = B*T;

        %% Part 2: Find boundary edges and modify the right hand side f and g
        % Find boundary edges: Neumann and Robin
        Neumann = []; Robin = []; %#ok<*NASGU>
        if ~isempty(bdFlag)
            isNeumann(elem2dof((bdFlag(:)==2)|(bdFlag(:) == 3))) = true;
            isRobin(elem2dof(bdFlag(:)==3)) = true;
            Neumannidx = find(isNeumann);
            Neumann   = edge(isNeumann,:);
            Robin     = edge(isRobin,:);
        end
        if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
            % no bdFlag, only pde.g_N or pde.g_R is given in the input
            [tempvar,Neumann] = findboundary(elem);
            if ~isempty(pde.g_R)
                Robin = Neumann;
            end
        end

        % Neumann boundary condition
        if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
            [lambda,w] = quadpts1(3);
            nQuad = size(lambda,1);
            % edge bases
            bdphi = 2*(lambda(:,1)+lambda(:,2))-1;
            % length of edge
            ve = node(Neumann(:,1),:) - node(Neumann(:,2),:);
            edgeLength = sqrt(sum(ve.^2,2));
            % update RHS
            for pp = 1:nQuad
                pxy = lambda(pp,1)*node(Neumann(:,1),:)+lambda(pp,2)*node(Neumann(:,2),:);
                gp = pde.g_N(pxy);
                f1(Neumannidx) = f1(Neumannidx) + w(pp)*edgeLength.*gp(:,1).*bdphi(pp); % interior bubble
                f2(Neumannidx) = f2(Neumannidx) + w(pp)*edgeLength.*gp(:,2).*bdphi(pp); % interior bubble
            end
        end
        f = [f1; f2];
        % The case non-empty Neumann but g_N=[] corresponds to the zero flux
        % boundary condition on Neumann edges and no modification is needed.

        % Dirichlet boundary conditions
        if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
            u1 = zeros(Nu,1);
            u2 = zeros(Nu,1);
            bdEdgeMid = (node(edge(fixedDof,1),:)+node(edge(fixedDof,2),:))/2;
            uD = pde.g_D(bdEdgeMid);         % bd values at middle points of edges
            u1(fixedDof) = uD(:,1);
            u2(fixedDof) = uD(:,2);
            u = [u1; u2]; % Dirichlet bd condition is built into u
            f = f - A*u;  % bring affect of nonhomgenous Dirichlet bd condition to
            g = g - B*u;  % the right hand side
            g = g - mean(g); % impose the compatible condition
            f(fixedDof) = u1(fixedDof);
            f(fixedDof+Nu) = u2(fixedDof);
            u = [u1; u2]; % Dirichlet bd condition is built into u
        end
        % The case non-empty Dirichlet but g_D=[] corresponds to the zero Dirichlet
        % boundary condition and no modification is needed.

        % modfiy pressure dof for pure Dirichlet
        if isempty(Neumann)
            pDof = (1:Np-1)';
        end

        ufreeDof = [ufreeDof; Nu+ufreeDof];
    end % end of function getbdStokesCR
end

repo: ifem


text:
function robot = parse_urdf(file)
% Loading file from urdf
% The function is tailored for UR robots: serial robots with 6 DOF
% Modify it if you have  a robot with different degrees of freedom
robot = xml2struct(file);
no_dof = 6;

% Extracting parameters of the robot
for i = 1:no_dof
    % axis of rotation of a joint i in coordinate system of joint i
    axis_of_rot = str2num(robot.robot.joint{i}.axis.Attributes.xyz)';
    % mass of link (i+1) because joint i rotates link (i+1) as the numbering of
    % links starts from base link that it not moving
    link_mass = str2double(robot.robot.link{i+1}.inertial.mass.Attributes.value);
    % poistion of the com in frame attached to link
    com_pos = str2num(robot.robot.link{i+1}.inertial.origin.Attributes.xyz)';
    com_vec2mat = vec2skewSymMat(com_pos);
    % inertial parameters of the link expressed in coordinate system attached
    % the center of mass.
    ixx = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.ixx);
    ixy = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.ixy);
    ixz = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.ixz);
    iyy = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.iyy);
    iyz = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.iyz);
    izz = str2double(robot.robot.link{i+1}.inertial.inertia.Attributes.izz);
    % the inertia tensor wrt the frame oriented as the body frame and with the
    % origin in the COM
    link_inertia = [ixx, ixy, ixz; ixy, iyy iyz; ixz, iyz, izz];
    % manipulator regressor
    robot.m(i) = link_mass;
    robot.k(:,i) = axis_of_rot;
    robot.r_com(:,i) = com_pos;
    robot.I(:,:,i) = link_inertia;
    robot.h(:,i) = link_mass*com_pos;
    robot.I_vec(:,i) = inertiaMatrix2Vector(link_inertia-...
                            link_mass*com_vec2mat*com_vec2mat);
    robot.pi(:,i) = [robot.I_vec(:,i); robot.h(:,i); robot.m(i)];
end

author: shamilmamedov


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 1;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function  [rankFR, rankFRV, rankFRvanilla, rankFRVvanilla, model] = checkRankFR(model, printLevel)
% Calculates the rank of `[F R]` and `[F; R]`, when restricted to certain
% rows and columns
%
% USAGE:
%
%    [rankFR, rankFRV, rankFRvanilla, rankFRVvanilla, model] = checkRankFR(model, printLevel)
%
% INPUTS:
%    model:             model structure
%    printLevel:        verbose level
%
% OUTPUTS:
%    rankFR:            rank of [`F R`], when using only `FRrows`
%    rankFRV:           rank of [`F; R`], when using only `FRVcols`
%    rankFRvanilla:     rank of [`F R`], when using all rows
%    rankFRVvanilla:    rank of [`F; R`], when using all cols
%    model:             structure with fields:
%
%                         * .FRrows - `m x 1` boolean of rows of [`F R`] that are nonzero,
%                           unique upto positive scaling and part of the
%                           maximal conservation vector
%                         * .FRVcols - `n x 1` boolean of cols of [`F; R`] that are nonzero,
%                           unique upto positive scaling and part of the
%                           maximal conservation vector
%                         * .FRirows - `m x 1` boolean of rows of [`F R`] that are independent
%                         * .FRdrows - `m x 1` boolean of rows of [`F R`] that are dependent
%                         * .FRwrows - `m x 1` boolean of independent rows of [`F R`] that
%                           have dependent rows amongst `model.FRdrows`
%                         * .FRVdcols - `n x 1` boolean of cols of [`F; R`] that are dependent
%                         * model.SConsistentMetBool - `m x 1` boolean vector indicating metabolites involved
%                           in the maximal consistent vector
%                         * .SConsistentRxnBool - `n x 1` boolean vector indicating metabolites involved
%                           in the maximal consistent vector
%                         * .FRnonZeroBool - `m x 1` boolean vector indicating metabolites involved
%                           in at least one internal reaction
%                         * .FRuniqueBool - `m x 1` boolean vector indicating metabolites with
%                           reaction stoichiometry unique upto scalar multiplication
%                         * .SIntRxnBool - `n x 1` boolean vector indicating the non exchange reactions
%                         * .FRVnonZeroBool - `n x 1` boolean vector indicating non exchange reactions
%                           with at least one metabolite involved
%                         * .FRVuniqueBool - `n x 1` boolean vector indicating non exchange reactions
%                           with stoichiometry unique upto scalar multiplication
%                         * .connectedRowsFRBool - `m x 1` boolean vector indicating metabolites in connected rows of [`F R`]
%                         * .connectedRowsFRVBool - `n x 1` boolean vector indicating complexes in connected columns of [`F; R`]
%                         * .V - :math:`S V = 0`; :math:`1^T |V| > 1` for all flux consistent reactions

if ~exist('printLevel','var')
    printLevel=1;
end

%scaling of stoichiometric matrix may be necessary when badly scaled, off
%by default
scaling='none';
%scaling='gmscal';
%scaling='minimumCoefficient';
switch scaling
    case 'minimumCoefficient'
        %rescale stoichiometric coefficients
        A=abs(model.S);
        minCoefficient=min(min(A(model.S~=0)));
        model.S=model.S*(1000/minCoefficient);
        model.lb=model.lb*(1000/minCoefficient);
        model.ub=model.ub*(1000/minCoefficient);
    case 'gmscal'
        %geometric mean scaling
        [cscale,rscale] = gmscal(model.S,0,0.9);

        % model.S =diag(rscale)*model.S*diag(cscale);
        % model.b =diag(rscale)*model.b;
        % model.lb=diag(cscale)*model.lb;
        % model.ub=diag(cscale)*model.ub;
        % model.c =diag(cscale)*model.c;

        % only scale columns
        model.S =model.S*diag(cscale);
        model.b =model.b;
        model.lb=diag(cscale)*model.lb;
        model.ub=diag(cscale)*model.ub;
        model.c =diag(cscale)*model.c;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%vanilla forward and reverse half stoichiometric matrices
F       = -model.S;
F(F<0)  =    0;
R       =  model.S;
R(R<0)  =    0;

%vanilla ranks
[rankFRvanilla,~,~]      = getRankLUSOL([F R]);
[rankFRVvanilla,~,~]     = getRankLUSOL([F;R]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nMet,nRxn]=size(model.S);

if ~isfield(model,'SIntRxnBool')  || ~isfield(model,'SIntMetBool')
    %finds the reactions in the model which export/import from the model
    %boundary i.e. mass unbalanced reactions
    %e.g. Exchange reactions
    %     Demand reactions
    %     Sink reactions
    model = findSExRxnInd(model);
end
model.SIntRxnBool_findSExRxnInd=model.SIntRxnBool;

%mass and charge balance
if isfield(model,'metFormulas')
    [massImbalance,imBalancedMass,imBalancedCharge,imBalancedRxnBool,Elements,missingFormulaeBool,balancedMetBool]...
        = checkMassChargeBalance(model,printLevel);
    model.balancedRxnBool=~imBalancedRxnBool;
    model.balancedMetBool=balancedMetBool;
    model.Elements=Elements;
    model.missingFormulaeBool=missingFormulaeBool;

    if nnz(model.missingFormulaeBool)<size(model.S,1)
        %assumes that all mass imbalanced reations are exchange reactions
        model.SIntRxnBool = model.SIntRxnBool & model.balancedRxnBool;
        model.SIntMetBool = model.SIntMetBool & model.balancedMetBool;
    end
end

%stoichiometric consistency
if ~isfield(model,'SConsistentMetBool') || ~isfield(model,'SConsistentRxnBool')
    %finds the metabolites in the model that are stoichiometrically
    %consistent

    % method.interface='solveCobraLP';
    % method.solver='mosek';
    % method.param.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_DUAL_SIMPLEX';
    method.interface='solveCobraLP';
    method.solver='gurobi';
    method.param.Method=1;

    [inform,m,model]=checkStoichiometricConsistency(model,printLevel-1,method);

    %assumes that all stoichiometrically inconsistent reactions are exchange reactions
    model.SIntRxnBool = model.SIntRxnBool & model.SConsistentRxnBool;
    model.SIntMetBool = model.SIntMetBool & model.SConsistentMetBool;
end

if printLevel>0
    fprintf('%u%s\n',nnz(~model.SIntRxnBool),' exchange reactions (imbalanced elementally, or involves a stoichiometrically inconsistent molecular species)')
    for i=1:nRxn
        if ~model.SIntRxnBool(i) && model.SIntRxnBool_findSExRxnInd(i)
            fprintf('%s%s\n',model.rxns{i},': deemed an exchange reaction')
        end
    end
    fprintf('---------\n')
end

if any(~model.SIntMetBool) & printLevel>0
    fprintf('%u%s\n',nnz(~model.SIntMetBool),' rows of S only involved in exchange reactions')
    for i=1:nMet
        if ~model.SIntMetBool(i)
            fprintf('%s%s\n',model.mets{i}',': deemed a molecular species only involved in exchange reactions')
        end
    end
    fprintf('---------\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%metBool1=(model.SConsistentMetBool | ~model.SIntMetBool); %incorrect to
%include mets that are only exchanged
metBool1=model.SConsistentMetBool;
rxnBool1=(model.SConsistentRxnBool | ~model.SIntRxnBool);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find rows that are not all zero when a subset of reactions omitted
A1=[F(:,rxnBool1) R(:,rxnBool1)];
model.FRnonZeroRowBool1 = any(A1,2);

%find cols that are not all zero when a subset of metabolites omitted
A1=[F(metBool1,:); R(metBool1,:)];
model.FRnonZeroColBool1 = any(A1,1)';

%only report for the consistent rows
if any(~model.FRnonZeroRowBool1 & metBool1) & printLevel>0
    fprintf('%u%s\n',nnz(~model.FRnonZeroRowBool1 & metBool1),' zero rows of [F,R]')
    for i=1:nMet
        if ~model.FRnonZeroRowBool1(i) && metBool1(i)
            fprintf('%s%s\n',model.mets{i},': is a zero row of [F,R]')
        end
    end
    fprintf('\n')
end

%only report for the consistent cols
if any(~model.FRnonZeroColBool1 & rxnBool1) && 0
    fprintf('%u%s\n',nnz(~model.FRnonZeroColBool1 & rxnBool1),' zero cols of consistent [F;R]')
    for i=1:nRxn
        if ~model.FRnonZeroColBool1(i) && rxnBool1(i)
            fprintf('%s%s\n',model.rxns{i},': is a zero col of consistent [F;R]')
        end
    end
    fprintf('\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=[F,R];%depends on the direction of reaction

%detect the rows of A that are identical upto scalar multiplication
%divide each row by the sum of each row.
sumA2            = sum(A,2);
sumA2(sumA2==0) = 1;
normalA2          = diag(1./sumA2)*A;

%get unique rows, but do not change the order
% [C,IA,IC] = unique(A,'rows') also returns index vectors IA and IC such
% that C = A(IA,:) and A = C(IC,:).
[uniqueRowsA,IA,IC] = unique(normalA2,'rows','stable');
model.FRuniqueRowBool=zeros(nMet,1);
model.FRuniqueRowBool(IA)=1;

if any(~model.FRuniqueRowBool & metBool1) & printLevel>0
    fprintf('%u%s\n',nnz(~model.FRuniqueRowBool & metBool1),' non-unique rows of stoich consistent [F,R]')
    for i=1:nMet
        if ~model.FRuniqueRowBool(i) && metBool1(i)
            fprintf('%s%s\n',model.mets{i},': not unique row of stoich consistent [F,R]')
        end
    end
    fprintf('\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=F+R;%invariant to direction of reaction

%detect the cols of A that are identical upto scalar multiplication
%divide each col by the sum of each row.
sumA1          = sum(A,1);
sumA1(sumA1==0)  = 1;
normalA1          = A*diag(1./sumA1);

%get unique cols, but do not change the order
% [C,IA,IC] = unique(A,'rows') also returns index vectors IA and IC such
% that C = A(IA,:) and A = C(IC,:).
[uniqueColsA,IA,IC] = unique(normalA1','rows','stable');
model.FRuniqueColBool=zeros(nRxn,1);
model.FRuniqueColBool(IA)=1;

if any(~model.FRuniqueColBool & rxnBool1) & printLevel>0
    fprintf('%u%s\n',nnz(~model.FRuniqueColBool & rxnBool1),' non-unique cols of stoich consistent [F;R]')
    for i=1:nRxn
        if ~model.FRuniqueColBool(i) && rxnBool1(i)
            fprintf('%s%s\n',model.rxns{i},': not unique col of stoich consistent [F;R]')
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if printLevel>0
    fprintf('\n%s\n','Diagnostics on size of boolean vectors for rows:')
    fprintf('%s%u\n','SIntMetBool:                 ',nnz(model.SIntMetBool))
    fprintf('%s%u\n','SConsistentMetBool:          ',nnz(model.SConsistentMetBool))
    fprintf('%s%u\n','FRnonZeroRowBool1:               ',nnz(model.FRnonZeroRowBool1))
    fprintf('%s%u\n','FRuniqueRowBool:               ',nnz(model.FRuniqueRowBool))

    fprintf('\n%s\n','Diagnostics on size of boolean vectors for cols:')
    fprintf('%s%u\n','SIntRxnBool:                 ',nnz(model.SIntRxnBool))
    fprintf('%s%u\n','SConsistentRxnBool:          ',nnz(model.SConsistentRxnBool))
    fprintf('%s%u\n','FRnonzeroColBool:               ',nnz(model.FRnonZeroColBool1))
    fprintf('%s%u\n','FRuniqueColBool:               ',nnz(model.FRuniqueColBool))
    fprintf('\n')
end

if ~isempty(find(strcmp(model.mets(~model.SConsistentMetBool),'h[c]')))
    warning('h[c] is one of the inconsistent molecules')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
metBool2=model.SConsistentMetBool & model.FRuniqueRowBool & model.FRnonZeroRowBool1;

% rxnBool2=(model.SConsistentRxnBool | ~model.SIntRxnBool);
rxnBool2=(model.SConsistentRxnBool | ~model.SIntRxnBool) & model.FRuniqueColBool;

%the uniqueness check has to be done before flux consistency since some
%models have a forward and backward reaction in separately, but these are
%the same reaction when both considered reversible.

%from the unique exchange reactions, and the subset of the unique non-exchange
%reactions that form a matrix with stoichiometrically consistent rows,
%check for reactions that are also flux consistent, with all reactions
%reversible
modelRev.S = model.S(metBool2,rxnBool2);

%infinite bounds may be necessary when stoichiometric matrix is badly
%scaled, in that case the minimum flux required (epsilon) is set to 1, the
%default lower and upper bounds are {-1000, 1000}, with epsilon getCobraSolverParams('LP', 'feasTol')*100,
%as that is two orders of magnitude higher than the typical feasibility
%tolerance for volation of mass balance constraints.
fluxConsistencyBounds='aThousand';
switch fluxConsistencyBounds
    case 'aThousand'
        %use all reactions reversible
        modelRev.lb=-1000*ones(size(modelRev.S,2),1);
        modelRev.ub= 1000*ones(size(modelRev.S,2),1);
        epsilon = getCobraSolverParams('LP', 'feasTol')*100;
    case 'infinity'
        %use all reactions reversible
        modelRev.lb=-inf*ones(size(modelRev.S,2),1);
        modelRev.ub= inf*ones(size(modelRev.S,2),1);
        epsilon = 1;
    case 'reconstruction'
        %use reconstruction reactions
        modelRev.lb=model.lb(rxnBool2);
        modelRev.ub=model.ub(rxnBool2);
        epsilon = getCobraSolverParams('LP', 'feasTol')*100;
end
modelRev.mets=model.mets(metBool2);
modelRev.rxns=model.rxns(rxnBool2);

%fast consistency check code from Nikos Vlassis et al
modeFlag=1;
[indFluxConsist,~,V0]=fastcc(modelRev,epsilon,printLevel-1,modeFlag);
modelRev.fluxConsistentRxnBool=false(size(modelRev.S,2),1);
modelRev.fluxConsistentRxnBool(indFluxConsist)=1;

%build the vector the same size as the original S
model.fluxConsistentRxnBool=false(nRxn,1);
%assign the boolean state to the original size model
model.fluxConsistentRxnBool(rxnBool2)=modelRev.fluxConsistentRxnBool;

%pad out V0 to correspond to the original size of S
model.V=sparse(nRxn,size(V0,2));
model.V(rxnBool2,:)=V0;

%check to make sure booleans are correct
if nnz(model.SConsistentRxnBool & ~model.SIntRxnBool)~=0
    error('No exchange reaction should be stoichiometrically consistent')
end

%check to make sure that V in nullspace of S.
tmp=norm(full(model.S(metBool2,model.SConsistentRxnBool & model.SIntRxnBool)*model.V(model.SConsistentRxnBool & model.SIntRxnBool,:)...
       + model.S(metBool2,~model.SIntRxnBool)*model.V(~model.SIntRxnBool,:)));
%two extra digits spare
if tmp>100*epsilon
    disp(tmp)
    error('Not flux consistent')
end
%create a stoichiometric matrix of perpetireactions only for
%stoichiometrically and flux consistent part
model.P=sparse(nMet,size(V0,2));
model.P(metBool2,:)=model.S(metBool2,~model.SIntRxnBool)*model.V(~model.SIntRxnBool,:);

% model.fluxConsistentRxnBool=true(nRxn,1);

%rows corresponding to flux consistent reactions
model.fluxConsistentMetBool = sum(model.S(:,model.fluxConsistentRxnBool)~=0,2)~=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%eliminate the metabolites and reactions that are stoichiometrically
%inconsistent or flux inconsistent from further consideration, but keep the
%flux consistent exchange reactions, also eliminate scalar multiples
metBool3=model.SConsistentMetBool & model.fluxConsistentMetBool & model.FRuniqueRowBool & model.FRnonZeroRowBool1;
rxnBool3=(model.SConsistentRxnBool | ~model.SIntRxnBool)  & model.FRuniqueColBool & model.fluxConsistentRxnBool;

%find rows that are not all zero when a subset of reactions omitted
A3=[F(:,rxnBool3) R(:,rxnBool3)];
model.FRnonZeroRowBool = any(A3,2);

%find cols that are not all zero when a subset of metabolites omitted
A3=[F(metBool3,:); R(metBool3,:)];
model.FRnonZeroColBool = any(A3,1)';

%only report for the latest subset of rows
if any(~model.FRnonZeroRowBool & metBool3) & printLevel>0
    fprintf('%u%s\n',nnz(~model.FRnonZeroRowBool & metBool3),' zero rows of [F,R]')
    for i=1:nMet
        if ~model.FRnonZeroRowBool(i) && metBool3(i)
            fprintf('%s%s\n',model.mets{i},': is a zero row of [F,R]')
        end
    end
    fprintf('\n')
end

%only report for the latest subset of cols
if any(~model.FRnonZeroColBool & rxnBool3) && 0
    fprintf('%u%s\n',nnz(~model.FRnonZeroColBool & rxnBool3),' zero cols of consistent [F;R]')
    for i=1:nRxn
        if ~model.FRnonZeroColBool(i) && rxnBool3(i)
            fprintf('%s%s\n',model.rxns{i},': is a zero col of consistent [F;R]')
        end
    end
    fprintf('\n')
end

%pause(eps)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A4=F(:,rxnBool3)+R(:,rxnBool3);%invariant to direction of reaction
%
% %detect the cols of A that are identical upto scalar multiplication
% %divide each col by the sum of each row.
% sumA4             = sum(A4,1);
% sumA4(sumA4==0)   = 1;
% normalA4          = A4*diag(1./sumA4);
%
% %get unique cols, but do not change the order
% % [C,IA,IC] = unique(A,'rows') also returns index vectors IA and IC such
% % that C = A(IA,:) and A = C(IC,:).
% [uniqueColsA,IA,IC] = unique(normalA4','rows','stable');
% model.FRuniqueColBool2=zeros(nRxn,1);
% model.FRuniqueColBool2(IA)=1;
%
% if any(~model.FRuniqueColBool2 & rxnBool3)
%     fprintf('%u%s\n',nnz(~model.FRuniqueColBool2 & rxnBool3),' non-unique cols of consistent [F;R]')
%     for i=1:nRxn
%         if ~model.FRuniqueColBool(i) && rxnBool1(i)
%             fprintf('%s%s\n',model.rxns{i},': not unique col of consistent [F;R]')
%         end
%     end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The largest connected component - not necessary, flux consistency takes
% care of this connected rows of [F,R] and columns of [F;R]
% [model.connectedRowsFRBool,model.connectedColsFRVBool]=connectedFR(F,R);
% [~,maxInd]=max(sum(model.connectedRowsFRBool,1));
% model.largestConnectedRowsFRBool=model.connectedRowsFRBool(:,maxInd);
% [~,maxInd]=max(sum(model.connectedColsFRVBool,1));
% model.largestConnectedColsFRVBool=model.connectedColsFRVBool(:,maxInd);

%bypass test for connectedness
model.largestConnectedRowsFRBool=true(nMet,1);
model.largestConnectedColsFRVBool=true(nRxn,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%diagnostics
if printLevel>0
    fprintf('\n%s\n','Diagnostics on size of boolean vectors for rows:')
    fprintf('%s%u\n','SConsistentMetBool:         ',nnz(model.SConsistentMetBool))
    fprintf('%s%u\n','fluxConsistentMetBool:      ',nnz(model.fluxConsistentMetBool))
    fprintf('%s%u\n','FRuniqueRowBool:               ',nnz(model.FRuniqueRowBool))
    fprintf('%s%u\n','FRnonZeroRowBool1:              ',nnz(model.FRnonZeroRowBool1))
    fprintf('%s%u\n','FRnonZeroRowBool:              ',nnz(model.FRnonZeroRowBool))
   % fprintf('%s%u\n','largestConnectedRowsFRBool: ',nnz(model.largestConnectedRowsFRBool))

    fprintf('\n%s\n','Diagnostics on size of boolean vectors for cols:')
    fprintf('%s%u\n','SIntRxnBool:                 ',nnz(model.SIntRxnBool))
    fprintf('%s%u\n','SConsistentRxnBool:          ',nnz(model.SConsistentRxnBool))
    fprintf('%s%u\n','fluxConsistentRxnBool:       ',nnz(model.fluxConsistentRxnBool))
    fprintf('%s%u\n','FRuniqueColBool:               ',nnz(model.FRuniqueColBool))
    fprintf('%s%u\n','FRnonZeroColBool1:              ',nnz(model.FRnonZeroColBool1))
    fprintf('%s%u\n','FRnonZeroColBool:              ',nnz(model.FRnonZeroColBool))
    %fprintf('%s%u\n','largestConnectedColsFRVBool: ',nnz(model.largestConnectedColsFRVBool))
    fprintf('\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%only use rows that are nonzero, unique upto positive scaling, part of
%the maximal conservation vector, and part of the largest component
selection='c';
switch selection
    case 'a'
        model.FRrows = model.SConsistentMetBool...
            & model.fluxConsistentMetBool...
            & model.FRnonZeroRowBool...
            & model.FRuniqueRowBool;
        model.FRVcols = (model.SConsistentRxnBool | ~model.SIntRxnBool)...
            & model.fluxConsistentRxnBool...
            & model.FRuniqueColBool...
            & model.FRnonZeroColBool;
    case 'b'
        %exchange reactions not considered when rank([F,R]) measured
        model.FRrows = model.SConsistentMetBool...
            & model.fluxConsistentMetBool...
            & model.FRnonZeroRowBool...
            & model.FRuniqueRowBool;
        model.FRVcols = model.SConsistentRxnBool...
            & model.fluxConsistentRxnBool...
            & model.FRuniqueColBool...
            & model.FRnonZeroColBool;
    case 'c'
        model.FRrows = model.SConsistentMetBool...
            & model.fluxConsistentMetBool...
            & model.FRnonZeroRowBool1...
            & model.FRnonZeroRowBool...
            & model.FRuniqueRowBool;
        model.FRVcols = (model.SConsistentRxnBool | ~model.SIntRxnBool)...
            & model.fluxConsistentRxnBool...
            & model.FRuniqueColBool...
            & model.FRnonZeroColBool;
    case 'd'
        model.FRrows = model.SConsistentMetBool...
            & model.fluxConsistentMetBool...
            & model.FRnonZeroRowBool1...
            & model.FRnonZeroRowBool...
            & model.FRuniqueRowBool;
        model.FRVcols = (model.SConsistentRxnBool | ~model.SIntRxnBool)...
            & model.fluxConsistentRxnBool;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%omit both the rows and the columns as specified above
Fr = F(model.FRrows,model.FRVcols);
Rr = R(model.FRrows,model.FRVcols);
Fc = F(model.FRrows,model.FRVcols);
Rc = R(model.FRrows,model.FRVcols);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% method.interface='solveCobraLP';
% method.solver='mosek';
% method.param.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_DUAL_SIMPLEX';

method.interface='solveCobraLP';
method.solver='gurobi5';
method.param.Method=1;

modelCheck.S = Rr-Fr;
modelCheck.SIntRxnBool=model.SIntRxnBool;
modelCheck.SIntMetBool=model.SIntMetBool(model.FRrows);
[informC,mC,~]=checkStoichiometricConsistency(modelCheck,0,method);
%save check on consistency
model.mC=mC;
if informC~=1
    disp(length(mC)-nnz(mC))
    error('(R-F) not stoichiometrically consistent')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [A,B,C]=bilinearDecomposition(model.S(model.FRrows,model.FRVcols));
% %forward and reverse half stoichiometric matrices
% Fb        =   -B;
% Fb(Fb<0)  =    0;
% Rb        =    B;
% Rb(Rb<0)  =    0;
% Fb(Fb~=0) =    1;
% Rb(Rb~=0) =    1;
% adj = inc2adj([Fb Rb]);
%
% % Test whether a graph is bipartite, if yes, return the two vertex sets
% % Inputs: graph in the form of adjancency list (neighbor list, see adj2adjL.m)
% % Outputs: True/False (boolean), empty set (if False) or two sets of vertices
% tic
% [isit,partA,partB]=isbipartite(adj2adjL(adj));
% toc
% %pause(eps)

%check rank of bilinear version
[A,B,C]=bilinearDecomposition(Fr-Rr);
%forward and reverse half stoichiometric matrices
Frb        =   -B;
Frb(Frb<0)  =    0;
Rrb        =    B;
Rrb(Rrb<0)  =    0;

%rank of [F R]
[rankBilinearFrRr,bilinearFrRrp,bilinearFrRrq] = getRankLUSOL([Frb Rrb]);
model.Frb=Frb;
model.Rrb=Rrb;
model.rankBilinearFrRr=rankBilinearFrRr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%rank of [F R]
[rankFR,FRp,FRq] = getRankLUSOL([Fr Rr]);

%indices of rows that are dependent
ind=find(model.FRrows);
iR=ind(FRp(1:rankFR));
dR=ind(FRp(rankFR+1:length(FRp)));
model.FRdrows=false(nMet,1);
model.FRdrows(dR)=1;
model.FRirows=false(nMet,1);
model.FRirows(iR)=1;
model.FRq=FRq;
model.FRp=FRp;
model.Fr=Fr;
model.Rr=Rr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rank of [F;R]
[rankFRV,FRVp,FRVq]  = getRankLUSOL([Fc;Rc]);

%boolean of non-exchange columns that are dependent
vertInd=find(model.FRVcols);
iC=vertInd(FRVq(1:rankFRV));
dC=vertInd(FRVq(rankFRV+1:length(FRVq)));
model.FRVdcols=false(nRxn,1);
model.FRVdcols(dC)=1;
model.FRVicols=false(nRxn,1);
model.FRVicols(iC)=1;
model.FRVq=FRVq;
model.FRVp=FRVp;
model.Fc=Fc;
model.Rc=Rc;

%sanity check
if nnz(model.FRrows) ~= nnz(model.FRrows & model.FRdrows)
    error('')
end
if nnz(model.FRVcols) ~= nnz(model.FRVcols & model.FRVdcols)
    error('')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%examine the dependencies between rows and columns if any
model.FRrowRankDeficiency=nnz(model.FRrows)-rankFR;
if model.FRrowRankDeficiency>0
    T  = [Fr Rr];
    %code from Michael
    C1 = T(FRp(1:rankFR),:);
    C2 = T(FRp(rankFR+1:length(FRp)),:);
    W0 = C1' \ C2'; % solves LS problems min ||C1'*W_k - C2_k'||

    model.FRW=sparse(model.FRrowRankDeficiency,nMet);
    model.FRW(:,iR)=W0';%L
    model.FRW(:,dR)=-speye(model.FRrowRankDeficiency);

    %indices of independent rows that other rows are dependent on
    wR=find(sum(abs(model.FRW),1)>0)';
    model.FRwrows=false(nMet,1);
    model.FRwrows(wR)=1;

    %sanity checks
    disp(norm(C1'*W0 - C2',inf)) %should be zero
    disp(norm([W0',-speye(model.FRrowRankDeficiency)]*[C1;C2],inf)) %should be zero

    W = sparse(model.FRrowRankDeficiency,size([Fr Rr],1));
    W(:,FRp(1:rankFR))=W0';
    W(:,FRp(rankFR+1:length(FRp)))=-speye(model.FRrowRankDeficiency);
    disp(norm(W*[Fr Rr],inf)) %should be zero

    model.FRW=sparse(model.FRrowRankDeficiency,nMet);
    model.FRW(:,model.FRrows)=W;
    disp(norm(model.FRW*[F(:,model.FRVcols), R(:,model.FRVcols)],inf)) %should be zero
    %%%%% N B %%%%
    disp(norm(model.FRW*[F,R],inf)) % will not in general be zero as other columns of [F R] could contribute

    % nnz(W)
    % figure; spy(W)
    % figure; spy(abs(W) > 1e-3)  % say
end

model.FRcolRankDeficiency=nnz(model.FRVcols)-rankFRV;
if model.FRcolRankDeficiency>0
    VT  = [Fc;Rc];
    %code from Michael
    VC1  = VT(:,model.FRVq(1:rankFRV)); % independent columns
    VC2  = VT(:,model.FRVq(rankFRV+1:size(VT,2))); % dependent columns
    VW0  = VC1 \ VC2; % solves LS problems min ||VC1*WV_k - VC2_k||

    model.FRVW=sparse(nRxn,model.FRcolRankDeficiency);
    model.FRVW(iC,:)=VW0; % independent columns
    model.FRVW(dC,:)=-speye(model.FRcolRankDeficiency); % dependent columns

    %indices of independent cols that other cols are dependent on
    wC=find(sum(abs(model.FRVW),2)>0);
    model.FRVwcols=false(nRxn,1);
    model.FRVwcols(wC)=1;

    if printLevel>2
        %sanity check
        disp(norm([F(model.FRrows,:); R(model.FRrows,:)]*model.FRVW,inf))%should be zero
        %%%%% N B %%%%
        disp(norm([F;R]*model.FRVW,inf)) % will not in general be zero as other rows [F;R] could contribute
    end

    % nnz(VW0)
    % figure; spy(VW0)
    % figure; spy(abs(VW0) > 1e-3)  % say
end

author: opencobra

repo: cobratoolbox


text:
function prob_test0251 ( )

%*****************************************************************************80
%
%% TEST0251 tests BUFFON_LAPLACE_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0251\n' );
  fprintf ( 1, '  BUFFON_LAPLACE_PDF evaluates the Buffon-Laplace PDF,\n' );
  fprintf ( 1, '  the probability that, on a grid of cells of width A\n' );
  fprintf ( 1, '  and height B, a needle of length L, dropped at random,\n' );
  fprintf ( 1, '  will cross at least one grid line.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A         B         L        PDF\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5
    a = i;
    for j = 1 : 5
      b = j;

      for k = 0 : 5
        l = k * min ( a, b ) / 5.0;
        pdf = buffon_laplace_pdf ( a, b, l );
        fprintf ( 1, '  %8.4f  %8.4f  %8.4f  %14f\n', a, b, l, pdf );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [zt,model] = latentlssvm(varargin)
% Calculate the latent variables of the LS-SVM classifier at the given test data
%
% >> Zt = latentlssvm({X,Y,'classifier',gam,sig2,kernel}, {alpha,b}, Xt)
% >> Zt = latentlssvm({X,Y,'classifier',gam,sig2,kernel}, Xt)
% >> [Zt, model] = latentlssvm(model, Xt)
%
% The latent variables of a binary classifier are the continuous
% simulated values of the test data which are used to make the
% final classifications. The classification of a testpoint depends
% on whether the latent value exceeds the model's threshold (b). If
% appropriate, the model is trained by the standard procedure (trainlssvm) first.
%
% As an application example: crossvalidation can be based on the latent variables:
%
% >> cost = crossvalidate(model, X, Y, 10, 'mse', 'mean', 'original', 'trainlssvm', 'latentlssvm')
%
%
% Full syntax
%
%     1. Using the functional interface:
%
% >> Zt = latentlssvm({X,Y,type,gam,sig2,kernel,preprocess}, Xt)
%
%       Outputs
%         Zt            : Nt x m matrix with predicted latent simulated outputs
%       Inputs
%         X             : N x d matrix with the inputs of the training data
%         Y             : N x 1 vector with the outputs of the training data
%         type          : 'classifier' ('c')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         Xt            : Nt x d matrix with the inputs of the test data
%
%
%     2. Using the object oriented interface:
%
% >> [Zt, model] = latentlssvm(model, Xt)
%
%       Outputs
%         Zt       : Nt x m matrix with continuous latent simulated outputs
%         model(*) : Trained object oriented representation of the LS-SVM model
%       Inputs
%         model    : Object oriented representation of the LS-SVM model
%         Xt       : Nt x d matrix with the inputs of the test data
%
% See also:
%   trainlssvm, simlssvm

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

model = varargin{1};
if iscell(model),
  model = initlssvm(model{:});
end

if model.type(1)~='c',
  error('Only usefull for classification tasks...');
end
[~, zt, model] = simlssvm(varargin{:});

author: peterhcharlton

repo: RRest


text:
function mtrPDB2TrackVis(volFile, trkExFile, pdbFile, trkFile)
% TrackVis file is in mm coordinates, but without center shift
% PDB file is loaded in AcPc coordinates.

% Load vol file to get AcPc xform
vol = niftiRead(volFile);
xformFrom = vol.qto_ijk;
xformTo = abs(vol.qto_xyz);
xformTo(1:3,4)=0;
xformFromAcPc = xformTo*xformFrom;

% Get mm dimensions of image
dim_mm = vol.dim(:).*abs(diag(vol.qto_xyz(1:3,1:3)));

% Load trks example header file
trks.header = read_trk_hdr(trkExFile);

% Load PDB file
fg = mtrImportFibers(pdbFile);

% Get the fiber coords into mm (without AcPc shift)
fg = dtiXformFiberCoords(fg,xformFromAcPc);

% Get into trk coordinates
trks.fiber = {};
for ll=1:length(fg.fibers)
    trks.fiber{ll}.num_points = size(fg.fibers{ll},2);
    trks.fiber{ll}.points = fg.fibers{ll}';
end

trks.header.n_count = length(fg.fibers);
write_trk(trks,trkFile);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = colorspace(Conversion,varargin)
%   COLORSPACE  Convert a color image between color representations.
%   B = COLORSPACE(S,A) converts the color representation of image A
%   where S is a string specifying the conversion.  S tells the
%   source and destination color spaces, S = 'dest<-src', or
%   alternatively, S = 'src->dest'.  Supported color spaces are
%
%     'RGB'              R'G'B' Red Green Blue (ITU-R BT.709 gamma-corrected)
%     'YPbPr'            Luma (ITU-R BT.601) + Chroma
%     'YCbCr'/'YCC'      Luma + Chroma ("digitized" version of Y'PbPr)
%     'YUV'              NTSC PAL Y'UV Luma + Chroma
%     'YIQ'              NTSC Y'IQ Luma + Chroma
%     'YDbDr'            SECAM Y'DbDr Luma + Chroma
%     'JPEGYCbCr'        JPEG-Y'CbCr Luma + Chroma
%     'HSV'/'HSB'        Hue Saturation Value/Brightness
%     'HSL'/'HLS'/'HSI'  Hue Saturation Luminance/Intensity
%     'XYZ'              CIE XYZ
%     'Lab'              CIE L*a*b* (CIELAB)
%     'Luv'              CIE L*u*v* (CIELUV)
%     'Lch'              CIE L*ch (CIELCH)
%
%  All conversions assume 2 degree observer and D65 illuminant.  Color
%  space names are case insensitive.  When R'G'B' is the source or
%  destination, it can be omitted. For example 'yuv<-' is short for
%  'yuv<-rgb'.
%
%  MATLAB uses two standard data formats for R'G'B': double data with
%  intensities in the range 0 to 1, and uint8 data with integer-valued
%  intensities from 0 to 255.  As MATLAB's native datatype, double data is
%  the natural choice, and the R'G'B' format used by colorspace.  However,
%  for memory and computational performance, some functions also operate
%  with uint8 R'G'B'.  Given uint8 R'G'B' color data, colorspace will
%  first cast it to double R'G'B' before processing.
%
%  If A is an Mx3 array, like a colormap, B will also have size Mx3.
%
%  [B1,B2,B3] = COLORSPACE(S,A) specifies separate output channels.
%  COLORSPACE(S,A1,A2,A3) specifies separate input channels.

%  Pascal Getreuer 2005-2006

%%% Input parsing %%%
if nargin < 2, error('Not enough input arguments.'); end
[SrcSpace,DestSpace] = parse(Conversion);

if nargin == 2
   Image = varargin{1};
elseif nargin >= 3
   Image = cat(3,varargin{:});
else
   error('Invalid number of input arguments.');
end

FlipDims = (size(Image,3) == 1);

if FlipDims, Image = permute(Image,[1,3,2]); end
if ~isa(Image,'double'), Image = double(Image)/255; end
if size(Image,3) ~= 3, error('Invalid input size.'); end

SrcT = gettransform(SrcSpace);
DestT = gettransform(DestSpace);

if ~ischar(SrcT) & ~ischar(DestT)
   % Both source and destination transforms are affine, so they
   % can be composed into one affine operation
   T = [DestT(:,1:3)*SrcT(:,1:3),DestT(:,1:3)*SrcT(:,4)+DestT(:,4)];
   Temp = zeros(size(Image));
   Temp(:,:,1) = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   Temp(:,:,2) = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   Temp(:,:,3) = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   Image = Temp;
elseif ~ischar(DestT)
   Image = rgb(Image,SrcSpace);
   Temp = zeros(size(Image));
   Temp(:,:,1) = DestT(1)*Image(:,:,1) + DestT(4)*Image(:,:,2) + DestT(7)*Image(:,:,3) + DestT(10);
   Temp(:,:,2) = DestT(2)*Image(:,:,1) + DestT(5)*Image(:,:,2) + DestT(8)*Image(:,:,3) + DestT(11);
   Temp(:,:,3) = DestT(3)*Image(:,:,1) + DestT(6)*Image(:,:,2) + DestT(9)*Image(:,:,3) + DestT(12);
   Image = Temp;
else
   Image = feval(DestT,Image,SrcSpace);
end

%%% Output format %%%
if nargout > 1
   varargout = {Image(:,:,1),Image(:,:,2),Image(:,:,3)};
else
   if FlipDims, Image = permute(Image,[1,3,2]); end
   varargout = {Image};
end

return;

function [SrcSpace,DestSpace] = parse(Str)
% Parse conversion argument

if isstr(Str)
   Str = lower(strrep(strrep(Str,'-',''),' ',''));
   k = find(Str == '>');

   if length(k) == 1         % Interpret the form 'src->dest'
      SrcSpace = Str(1:k-1);
      DestSpace = Str(k+1:end);
   else
      k = find(Str == '<');

      if length(k) == 1      % Interpret the form 'dest<-src'
         DestSpace = Str(1:k-1);
         SrcSpace = Str(k+1:end);
      else
         error(['Invalid conversion, ''',Str,'''.']);
      end
   end

   SrcSpace = alias(SrcSpace);
   DestSpace = alias(DestSpace);
else
   SrcSpace = 1;             % No source pre-transform
   DestSpace = Conversion;
   if any(size(Conversion) ~= 3), error('Transformation matrix must be 3x3.'); end
end
return;

function Space = alias(Space)
Space = strrep(Space,'cie','');

if isempty(Space)
   Space = 'rgb';
end

switch Space
case {'ycbcr','ycc'}
   Space = 'ycbcr';
case {'hsv','hsb'}
   Space = 'hsv';
case {'hsl','hsi','hls'}
   Space = 'hsl';
case {'rgb','yuv','yiq','ydbdr','ycbcr','jpegycbcr','xyz','lab','luv','lch'}
   return;
end
return;

function T = gettransform(Space)
% Get a colorspace transform: either a matrix describing an affine transform,
% or a string referring to a conversion subroutine
switch Space
case 'ypbpr'
   T = [0.299,0.587,0.114,0;-0.1687367,-0.331264,0.5,0;0.5,-0.418688,-0.081312,0];
case 'yuv'
   % R'G'B' to NTSC/PAL YUV
   % Wikipedia: http://en.wikipedia.org/wiki/YUV
   T = [0.299,0.587,0.114,0;-0.147,-0.289,0.436,0;0.615,-0.515,-0.100,0];
case 'ydbdr'
   % R'G'B' to SECAM YDbDr
   % Wikipedia: http://en.wikipedia.org/wiki/YDbDr
   T = [0.299,0.587,0.114,0;-0.450,-0.883,1.333,0;-1.333,1.116,0.217,0];
case 'yiq'
   % R'G'B' in [0,1] to NTSC YIQ in [0,1];[-0.595716,0.595716];[-0.522591,0.522591];
   % Wikipedia: http://en.wikipedia.org/wiki/YIQ
   T = [0.299,0.587,0.114,0;0.595716,-0.274453,-0.321263,0;0.211456,-0.522591,0.311135,0];
case 'ycbcr'
   % R'G'B' (range [0,1]) to ITU-R BRT.601 (CCIR 601) Y'CbCr
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   % Poynton, Equation 3, scaling of R'G'B to Y'PbPr conversion
   T = [65.481,128.553,24.966,16;-37.797,-74.203,112.0,128;112.0,-93.786,-18.214,128];
case 'jpegycbcr'
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   T = [0.299,0.587,0.114,0;-0.168736,-0.331264,0.5,0.5;0.5,-0.418688,-0.081312,0.5]*255;
case {'rgb','xyz','hsv','hsl','lab','luv','lch'}
   T = Space;
otherwise
   error(['Unknown color space, ''',Space,'''.']);
end
return;

function Image = rgb(Image,SrcSpace)
% Convert to Rec. 709 R'G'B' from 'SrcSpace'
switch SrcSpace
case 'rgb'
   return;
case 'hsv'
   % Convert HSV to R'G'B'
   Image = huetorgb((1 - Image(:,:,2)).*Image(:,:,3),Image(:,:,3),Image(:,:,1));
case 'hsl'
   % Convert HSL to R'G'B'
   L = Image(:,:,3);
   Delta = Image(:,:,2).*min(L,1-L);
   Image = huetorgb(L-Delta,L+Delta,Image(:,:,1));
case {'xyz','lab','luv','lch'}
   % Convert to CIE XYZ
   Image = xyz(Image,SrcSpace);
   % Convert XYZ to RGB
   T = [3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3);  % R
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3);  % G
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3);  % B
   % Desaturate and rescale to constrain resulting RGB values to [0,1]
   AddWhite = -min(min(min(R,G),B),0);
   Scale = max(max(max(R,G),B)+AddWhite,1);
   R = (R + AddWhite)./Scale;
   G = (G + AddWhite)./Scale;
   B = (B + AddWhite)./Scale;
   % Apply gamma correction to convert RGB to Rec. 709 R'G'B'
   Image(:,:,1) = gammacorrection(R);  % R'
   Image(:,:,2) = gammacorrection(G);  % G'
   Image(:,:,3) = gammacorrection(B);  % B'
otherwise  % Conversion is through an affine transform
   T = gettransform(SrcSpace);
   temp = inv(T(:,1:3));
   T = [temp,-temp*T(:,4)];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   AddWhite = -min(min(min(R,G),B),0);
   Scale = max(max(max(R,G),B)+AddWhite,1);
   R = (R + AddWhite)./Scale;
   G = (G + AddWhite)./Scale;
   B = (B + AddWhite)./Scale;
   Image(:,:,1) = R;
   Image(:,:,2) = G;
   Image(:,:,3) = B;
end

% Clip to [0,1]
Image = min(max(Image,0),1);
return;

function Image = xyz(Image,SrcSpace)
% Convert to CIE XYZ from 'SrcSpace'
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'xyz'
   return;
case 'luv'
   % Convert CIE L*uv to XYZ
   WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   L = Image(:,:,1);
   Y = (L + 16)/116;
   Y = invf(Y)*WhitePoint(2);
   U = Image(:,:,2)./(13*L + 1e-6*(L==0)) + WhitePointU;
   V = Image(:,:,3)./(13*L + 1e-6*(L==0)) + WhitePointV;
   Image(:,:,1) = -(9*Y.*U)./((U-4).*V - U.*V);                  % X
   Image(:,:,2) = Y;                                             % Y
   Image(:,:,3) = (9*Y - (15*V.*Y) - (V.*Image(:,:,1)))./(3*V);  % Z
case {'lab','lch'}
   Image = lab(Image,SrcSpace);
   % Convert CIE L*ab to XYZ
   fY = (Image(:,:,1) + 16)/116;
   fX = fY + Image(:,:,2)/500;
   fZ = fY - Image(:,:,3)/200;
   Image(:,:,1) = WhitePoint(1)*invf(fX);  % X
   Image(:,:,2) = WhitePoint(2)*invf(fY);  % Y
   Image(:,:,3) = WhitePoint(3)*invf(fZ);  % Z
otherwise   % Convert from some gamma-corrected space
   % Convert to Rec. 701 R'G'B'
   Image = rgb(Image,SrcSpace);
   % Undo gamma correction
   R = invgammacorrection(Image(:,:,1));
   G = invgammacorrection(Image(:,:,2));
   B = invgammacorrection(Image(:,:,3));
   % Convert RGB to XYZ
   T = inv([3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311]);
   Image(:,:,1) = T(1)*R + T(4)*G + T(7)*B;  % X
   Image(:,:,2) = T(2)*R + T(5)*G + T(8)*B;  % Y
   Image(:,:,3) = T(3)*R + T(6)*G + T(9)*B;  % Z
end
return;

function Image = hsv(Image,SrcSpace)
% Convert to HSV
Image = rgb(Image,SrcSpace);
V = max(Image,[],3);
S = (V - min(Image,[],3))./(V + (V == 0));
Image(:,:,1) = rgbtohue(Image);
Image(:,:,2) = S;
Image(:,:,3) = V;
return;

function Image = hsl(Image,SrcSpace)
% Convert to HSL
switch SrcSpace
case 'hsv'
   % Convert HSV to HSL
   MaxVal = Image(:,:,3);
   MinVal = (1 - Image(:,:,2)).*MaxVal;
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   Image(:,:,2) = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,3) = L;
otherwise
   Image = rgb(Image,SrcSpace);  % Convert to Rec. 701 R'G'B'
   % Convert R'G'B' to HSL
   MinVal = min(Image,[],3);
   MaxVal = max(Image,[],3);
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   S = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,1) = rgbtohue(Image);
   Image(:,:,2) = S;
   Image(:,:,3) = L;
end
return;

function Image = lab(Image,SrcSpace)
% Convert to CIE L*a*b* (CIELAB)
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'lab'
   return;
case 'lch'
   % Convert CIE L*CH to CIE L*ab
   C = Image(:,:,2);
   Image(:,:,2) = cos(Image(:,:,3)*pi/180).*C;  % a*
   Image(:,:,3) = sin(Image(:,:,3)*pi/180).*C;  % b*
otherwise
   Image = xyz(Image,SrcSpace);  % Convert to XYZ
   % Convert XYZ to CIE L*a*b*
   X = Image(:,:,1)/WhitePoint(1);
   Y = Image(:,:,2)/WhitePoint(2);
   Z = Image(:,:,3)/WhitePoint(3);
   fX = f(X);
   fY = f(Y);
   fZ = f(Z);
   Image(:,:,1) = 116*fY - 16;    % L*
   Image(:,:,2) = 500*(fX - fY);  % a*
   Image(:,:,3) = 200*(fY - fZ);  % b*
end
return;

function Image = luv(Image,SrcSpace)
% Convert to CIE L*u*v* (CIELUV)
WhitePoint = [0.950456,1,1.088754];
WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));

Image = xyz(Image,SrcSpace); % Convert to XYZ
U = (4*Image(:,:,1))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
V = (9*Image(:,:,2))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
Y = Image(:,:,2)/WhitePoint(2);
L = 116*f(Y) - 16;
Image(:,:,1) = L;                        % L*
Image(:,:,2) = 13*L.*(U - WhitePointU);  % u*
Image(:,:,3) = 13*L.*(V - WhitePointV);  % v*
return;

function Image = lch(Image,SrcSpace)
% Convert to CIE L*ch
Image = lab(Image,SrcSpace);  % Convert to CIE L*ab
H = atan2(Image(:,:,3),Image(:,:,2));
H = H*180/pi + 360*(H < 0);
Image(:,:,2) = sqrt(Image(:,:,2).^2 + Image(:,:,3).^2);  % C
Image(:,:,3) = H;                                        % H
return;

function Image = huetorgb(m0,m2,H)
% Convert HSV or HSL hue to RGB
N = size(H);
H = min(max(H(:),0),360)/60;
m0 = m0(:);
m2 = m2(:);
F = H - round(H/2)*2;
M = [m0, m0 + (m2-m0).*abs(F), m2];
Num = length(m0);
j = [2 1 0;1 2 0;0 2 1;0 1 2;1 0 2;2 0 1;2 1 0]*Num;
k = floor(H) + 1;
Image = reshape([M(j(k,1)+(1:Num).'),M(j(k,2)+(1:Num).'),M(j(k,3)+(1:Num).')],[N,3]);
return;

function H = rgbtohue(Image)
% Convert RGB to HSV or HSL hue
[M,i] = sort(Image,3);
i = i(:,:,3);
Delta = M(:,:,3) - M(:,:,1);
Delta = Delta + (Delta == 0);
R = Image(:,:,1);
G = Image(:,:,2);
B = Image(:,:,3);
H = zeros(size(R));
k = (i == 1);
H(k) = (G(k) - B(k))./Delta(k);
k = (i == 2);
H(k) = 2 + (B(k) - R(k))./Delta(k);
k = (i == 3);
H(k) = 4 + (R(k) - G(k))./Delta(k);
H = 60*H + 360*(H < 0);
H(Delta == 0) = nan;
return;

function Rp = gammacorrection(R)
Rp = real(1.099*R.^0.45 - 0.099);
i = (R < 0.018);
Rp(i) = 4.5138*R(i);
return;

function R = invgammacorrection(Rp)
R = real(((Rp + 0.099)/1.099).^(1/0.45));
i = (R < 0.018);
R(i) = Rp(i)/4.5138;
return;

function fY = f(Y)
fY = real(Y.^(1/3));
i = (Y < 0.008856);
fY(i) = Y(i)*(841/108) + (4/29);
return;

function Y = invf(fY)
Y = fY.^3;
i = (Y < 0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function organizeExperiments_Radiomics_Old(pathFeatures,pathLR,pathCR,outcomes,timeToEvent,nPatient,nonTextName,fSetNames,paramSEP,baselineSEP,freedomSEP,textType,textName)

startpath = pwd;

permTrain = randperm(nPatient.HGJ + nPatient.CHUS); permTest = randperm(nPatient.HMR + nPatient.CHUM); % Random permutation of instances.

nNonText = numel(nonTextName);
nFset = numel(fSetNames);
nameOutcomes = fieldnames(outcomes.HGJ); nOutcomes = numel(nameOutcomes);
nTextType = numel(textType);

training.nonText = struct; training.text = struct; training.outcomes = struct; training.timeToEvent = struct;
testing.nonText = struct; testing.text = struct; testing.outcomes = struct; testing.timeToEvent = struct;

% Organizing non-texture features
cd(pathFeatures)
train1 = load('nonText_HGJ_GTVtot'); train1 = struct2cell(train1); train1 = train1{1};
train2 = load('nonText_CHUS_GTVtot'); train2 = struct2cell(train2); train2 = train2{1};
test1 = load('nonText_HMR_GTVtot'); test1 = struct2cell(test1); test1 = test1{1};
test2 = load('nonText_CHUM_GTVtot'); test2 = struct2cell(test2); test2 = test2{1};
for i = 1:nNonText
    training.nonText.(nonTextName{i}).Data = [train1.(nonTextName{i}).Data;train2.(nonTextName{i}).Data];
    testing.nonText.(nonTextName{i}).Data = [test1.(nonTextName{i}).Data;test2.(nonTextName{i}).Data];
    training.nonText.(nonTextName{i}).Data = training.nonText.(nonTextName{i}).Data(permTrain);
    testing.nonText.(nonTextName{i}).Data = testing.nonText.(nonTextName{i}).Data(permTest);
end

% Organizing texture features
for f = 1:nFset
    training.text.(fSetNames{f}) = struct; testing.text.(fSetNames{f}) = struct;
    if strcmp(fSetNames{f},'PET') % Using only PET
        train1 = load('text_HGJ_PT_GTVtot'); train1 = struct2cell(train1); train1 = train1{1};
        train2 = load('text_CHUS_PT_GTVtot'); train2 = struct2cell(train2); train2 = train2{1};
        test1 = load('text_HMR_PT_GTVtot'); test1 = struct2cell(test1); test1 = test1{1};
        test2 = load('text_CHUM_PT_GTVtot'); test2 = struct2cell(test2); test2 = test2{1};
        train1 = {train1}; train2 = {train2};
        test1 = {test1}; test2 = {test2};
        training.text.(fSetNames{f}).param = paramSEP;
        training.text.(fSetNames{f}).baseline = baselineSEP;
        training.text.(fSetNames{f}).freedom = freedomSEP;
        training.text.(fSetNames{f}).paramName = {'Scale','Quant.algo','Ng'};
        textCellName = {'PT'};
    elseif strcmp(fSetNames{f},'CT') % Using only CT
        train1 = load('text_HGJ_CT_GTVtot'); train1 = struct2cell(train1); train1 = train1{1};
        train2 = load('text_CHUS_CT_GTVtot'); train2 = struct2cell(train2); train2 = train2{1};
        test1 = load('text_HMR_CT_GTVtot'); test1 = struct2cell(test1); test1 = test1{1};
        test2 = load('text_CHUM_CT_GTVtot'); test2 = struct2cell(test2); test2 = test2{1};
        train1 = {train1}; train2 = {train2};
        test1 = {test1}; test2 = {test2};
        training.text.(fSetNames{f}).param = paramSEP;
        training.text.(fSetNames{f}).baseline = baselineSEP;
        training.text.(fSetNames{f}).freedom = freedomSEP;
        training.text.(fSetNames{f}).paramName = {'Scale','Quant.algo','Ng'};
        textCellName = {'CT'};
    end

    for i = 1:numel(train1)
        sizeParam = size(train1{1});
        training.text.(fSetNames{f}).(textCellName{i}) = cell(sizeParam);
        testing.text.(fSetNames{f}).(textCellName{i}) = cell(sizeParam);
        for n = 1:numel(train1{1})
            for type = 1:nTextType
                for text = 1:numel(textName{type})
                    training.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = [train1{i}{n}.(textType{type}).(textName{type}{text}).Data;train2{i}{n}.(textType{type}).(textName{type}{text}).Data];
                    testing.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = [test1{i}{n}.(textType{type}).(textName{type}{text}).Data;test2{i}{n}.(textType{type}).(textName{type}{text}).Data];
                    training.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = training.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data(permTrain);
                    testing.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = testing.text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data(permTest);
                end
            end
        end
    end
end

% Organizing outcomes
for o = 1:nOutcomes
    training.outcomes.(nameOutcomes{o}) = [outcomes.HGJ.(nameOutcomes{o});outcomes.CHUS.(nameOutcomes{o})];
    testing.outcomes.(nameOutcomes{o}) = [outcomes.HMR.(nameOutcomes{o});outcomes.CHUM.(nameOutcomes{o})];
    training.timeToEvent.(nameOutcomes{o}) = [timeToEvent.HGJ.(nameOutcomes{o});timeToEvent.CHUS.(nameOutcomes{o})];
    testing.timeToEvent.(nameOutcomes{o}) = [timeToEvent.HMR.(nameOutcomes{o});timeToEvent.CHUM.(nameOutcomes{o})];
    training.outcomes.(nameOutcomes{o}) = training.outcomes.(nameOutcomes{o})(permTrain);
    testing.outcomes.(nameOutcomes{o}) = testing.outcomes.(nameOutcomes{o})(permTest);
    training.timeToEvent.(nameOutcomes{o}) = training.timeToEvent.(nameOutcomes{o})(permTrain);
    testing.timeToEvent.(nameOutcomes{o}) = testing.timeToEvent.(nameOutcomes{o})(permTest);
end

% Re-arranging by outcomes (same features for all outcomes here)
tempTraining = training; tempTesting = testing; clear training testing
for o = 1:nOutcomes
    training.(nameOutcomes{o}).nonText = tempTraining.nonText;
    training.(nameOutcomes{o}).text = tempTraining.text;
    training.(nameOutcomes{o}).outcome = tempTraining.outcomes.(nameOutcomes{o});
    training.(nameOutcomes{o}).timeToEvent = tempTraining.timeToEvent.(nameOutcomes{o});
    testing.(nameOutcomes{o}).nonText = tempTesting.nonText;
    testing.(nameOutcomes{o}).text = tempTesting.text;
    testing.(nameOutcomes{o}).outcome = tempTesting.outcomes.(nameOutcomes{o});
    testing.(nameOutcomes{o}).timeToEvent = tempTesting.timeToEvent.(nameOutcomes{o});
end

cd(pathLR), save('training','training'), save('testing','testing'), save('permTrain','permTrain'), save('permTest','permTest')
cd(pathCR), save('training','training'), save('testing','testing'), save('permTrain','permTrain'), save('permTest','permTest')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function diagnostic = bisection(varargin)
%BISECTION Solve simple quasi-convex MAXIMIZATION problem by bisection
%
%   DIAGNOSTIC = BISECTION(F,h,options,tolerance)
%
%    min   t
%    subject to
%            F(x,t) >=0
%
%   NOTES
%    It is assumed that the problem is quasi-convex in the scalar simple
%    variable t.
%
%    Lower and upper bounds are automatically derived
%    (you are adviced not to add simple bounds on the t-variable)
%
%    Default bisection tolerance 1e-5.
%
%    The algorithm simply solves a series of feasibility problems where the
%    parameter t is fixed, and hones in on an optimal value using bisection
%
%    It is recommended to explicitly set a solver. Otherwise YALMIP will
%    have to try to figure out a suitable solver for the feasibility
%    problems

solvertime = tic;
Constraints = varargin{1};
Objective = varargin{2};
if isequal(getbase(Objective),[0 -1])
    % User wants to maximize something, so we can reuse old code format
    varargin{2} = -Objective;
    options = varargin{3};
    options.bisection.switchedsign = 0;
    options.solver = options.bisection.solver;
    varargin{3} = options;
    diagnostic = bisection_core(varargin{:});
elseif isequal(getbase(Objective),[0 1])
    % User wants to minimize something. Rewrite as old max code
    options = varargin{3};
    options.bisection.switchedsign = 1;
    options.solver = options.bisection.solver;
    varargin{3} = options;
    Constraints = replace(Constraints, Objective, -Objective);
    varargin{1} = Constraints;
    diagnostic = bisection_core(varargin{:});
end
diagnostic.yalmiptime = toc(solvertime)-diagnostic.solvertime;
diagnostic.info = yalmiperror(diagnostic.problem,'bisection');

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ess,m] = spm_mci_ess (x,p)
% Compute Effective Sample Size
% FORMAT [ess,m] = spm_mci_ess (x,p)
%
% x      Univariate time series
% p      Maximum lag for autocovariance estimation
%
% ess    Effective Sample Size
% m      Number of lags used in ESS estimate
%
% This routine is based on the Initial Positive Sequence estimate
% proposed in C. Geyer (1992) Practical Markov Chain Monte Carlo,
% Statistical Science, 7(4):473-511.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_ess.m 6697 2016-01-27 14:57:28Z spm $

N=length(x);

try, pmax=p; catch, pmax=min(ceil(N/10),256); end

for i=1:pmax,
    y(:,i)=x(pmax-i+1:end-i);
end
C=cov(y);
c=C(1,:);
gamma=c(2:end);
gamma0=c(1);
r=gamma/gamma0;

G=[];
for j=1:floor(pmax/2)-1,
    % Sum of adjacent pairs of autocovariances
    G(j)=gamma(2*j)+gamma(2*j+1);
end

if ~isempty(G)
    % Find minimum j such that all G's up to j are positive
    Gneg=find(G<0);
    if isempty(Gneg)
        m=length(G);
    else
        m1=min(Gneg);
        m=m1-1;
    end
else
    m=0;
end

ess=N/(1+2*sum(r(1:2*m)));

% figure;
% plot(c);
% title('Autocovariance');
%
% figure
% plot(G);
% title('Sum of adjacent covariances');

author: spm


text:
function prod = tensor_array_innerprod( X1, X2 )

prod = 0;
N = length(X1);
for i = 1:N
    prod = prod + innerprod(X1{i},X2{i});
end

end

author: andrewssobral

repo: lrslibrary


text:
function exp_map_iteration(vt_id, transV, yawRotV, heightV, xGc, yGc, zGc, curYawHdc, curHeight)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    %% define some variables

    %%% some variables of em

    % define a variable of experiences
    global EXPERIENCES;

    % current experience id
    global CUR_EXP_ID;

    % the number of experiences
    global NUM_EXPS;

    % the min change of em between previous e with current e
    global MIN_DELTA_EM;

    % the experience delta threshold
    global DELTA_EXP_GC_HDC_THRESHOLD;

    %%% some variables of vt

    % define a variable of visual template
    global VT;

    % the previous visual template id
    global PREV_VT_ID;

    %%% some variables of 3D GC

    % The x, y, z dimension of 3D Grid Cells Model (3D CAN)
    global GC_X_DIM;
    global GC_Y_DIM;
    global GC_Z_DIM;

    %%% some variables of 3D HDC

    % The dimension of yaw in yaw_height_hdc network
    global YAW_HEIGHT_HDC_Y_DIM;

    % The dimension of height in yaw_height_hdc network
    global YAW_HEIGHT_HDC_H_DIM;

    %%% computing delta x,y,z, yaw, height

    % integrate the delta x, y, z, yaw, height
    % Rad radian  accumulative
    global ACCUM_DELTA_X;
    global ACCUM_DELTA_Y;
    global ACCUM_DELTA_Z;
    global ACCUM_DELTA_YAW;   % accum_delta_facing
%     global ACCUM_DELTA_HEIGHT;

    ACCUM_DELTA_YAW = clip_radian_180(ACCUM_DELTA_YAW + yawRotV);
%     ACCUM_DELTA_HEIGHT = mod(ACCUM_DELTA_HEIGHT + heightV, YAW_HEIGHT_HDC_H_DIM);
%
    ACCUM_DELTA_X = ACCUM_DELTA_X + transV * cos(ACCUM_DELTA_YAW);
    ACCUM_DELTA_Y = ACCUM_DELTA_Y + transV * sin(ACCUM_DELTA_YAW);
    ACCUM_DELTA_Z = ACCUM_DELTA_Z + heightV;

    % trajectory of delta of em
    global DELTA_EM;

    minDeltaX = get_min_delta(EXPERIENCES(CUR_EXP_ID).x_gc, xGc, GC_X_DIM);
    minDeltaY = get_min_delta(EXPERIENCES(CUR_EXP_ID).y_gc, yGc, GC_Y_DIM);
    minDeltaZ = get_min_delta(EXPERIENCES(CUR_EXP_ID).z_gc, zGc, GC_Z_DIM);

    minDeltaYaw = get_min_delta(EXPERIENCES(CUR_EXP_ID).yaw_hdc, curYawHdc, YAW_HEIGHT_HDC_Y_DIM);
    minDeltaHeight = get_min_delta(EXPERIENCES(CUR_EXP_ID).height_hdc, curHeight, YAW_HEIGHT_HDC_H_DIM);

    minDeltaYawReversed = get_min_delta(EXPERIENCES(CUR_EXP_ID).yaw_hdc, (YAW_HEIGHT_HDC_Y_DIM /2) - curYawHdc, YAW_HEIGHT_HDC_Y_DIM);
    minDeltaYaw = min(minDeltaYaw, minDeltaYawReversed);

%     minDeltaHeightReversed = get_min_delta(EXPERIENCES(CUR_EXP_ID).height_hdc, (YAW_HEIGHT_HDC_H_DIM /2) - curHeight, YAW_HEIGHT_HDC_H_DIM);
%     minDeltaHeight = min(minDeltaHeight, minDeltaHeightReversed);

    delta_em = sqrt((minDeltaX).^2 + (minDeltaY).^2 + (minDeltaZ).^2 + (minDeltaYaw).^2 + (minDeltaHeight).^2);
    DELTA_EM = [DELTA_EM delta_em];

    % if the visual template is new or the 3d grid cells (x,y,z) and head direction cells (yaw, height)
    % has changed enough create a new experience
    if VT(vt_id).numExp == 0 || delta_em > DELTA_EXP_GC_HDC_THRESHOLD

        NUM_EXPS = NUM_EXPS + 1;
        create_new_exp(CUR_EXP_ID, NUM_EXPS, vt_id, xGc, yGc, zGc, curYawHdc, curHeight);

        PREV_EXP_ID = CUR_EXP_ID;
        CUR_EXP_ID = NUM_EXPS;

        ACCUM_DELTA_X = 0;
        ACCUM_DELTA_Y = 0;
        ACCUM_DELTA_Z = 0;

        ACCUM_DELTA_YAW = EXPERIENCES(CUR_EXP_ID).yaw_exp_rad;
%         ACCUM_DELTA_HEIGHT = EXPERIENCES(CUR_EXP_ID).height_hdc;

        % if the visual template has changed (but isn't new) search for the matching experience
    elseif vt_id ~= PREV_VT_ID

        % find the experience associated with the current visual template and that is under the
        % threshold distance to the centre of grid cell and head direction cell activity
        % if multiple experiences are under the threshold then don't match (to reduce hash collisions)
        matched_exp_id = 0;
        matched_exp_count = 0;

        for search_id = 1:VT(vt_id).numExp

            minDeltaYaw = get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).yaw_hdc, curYawHdc, YAW_HEIGHT_HDC_Y_DIM);
%             minDeltaYawReversed = get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).yaw_hdc, (YAW_HEIGHT_HDC_Y_DIM /2) - curYawHdc, YAW_HEIGHT_HDC_Y_DIM);
%             minDeltaYaw = min(minDeltaYaw, minDeltaYawReversed);

            minDeltaHeight = get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).height_hdc, curHeight, YAW_HEIGHT_HDC_Y_DIM);
%             minDeltaHeightReversed = get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).height_hdc, (YAW_HEIGHT_HDC_Y_DIM /2) - curHeight, YAW_HEIGHT_HDC_H_DIM);
%             minDeltaHeight = min(minDeltaHeight, minDeltaHeightReversed);

            delta_em(search_id) = sqrt(get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).x_gc, xGc, GC_X_DIM)^2 ...
            + get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).y_gc, yGc, GC_Y_DIM)^2 ...
            + get_min_delta(EXPERIENCES(VT(vt_id).EXPERIENCES(search_id).id).z_gc, yGc, GC_Z_DIM)^2 ...
            + minDeltaYaw^2 + minDeltaHeight^2) ;

            if delta_em(search_id) < DELTA_EXP_GC_HDC_THRESHOLD
               matched_exp_count = matched_exp_count + 1;
            end
        end

        if matched_exp_count > 1
            % this means we aren't sure about which experience is a match due to hash table collision
            % instead of a false posivitive which may create blunder links in
            % the experience map keep the previous experience

        else
            [min_delta, min_delta_id] = min(delta_em);

            MIN_DELTA_EM = [MIN_DELTA_EM; min_delta];
            if min_delta < DELTA_EXP_GC_HDC_THRESHOLD

                matched_exp_id = VT(vt_id).EXPERIENCES(min_delta_id).id;

                % see if the previous experience already has a link to the current experience
                link_exists = 0;
                for link_id = 1 : EXPERIENCES(CUR_EXP_ID).numlinks
                    if EXPERIENCES(CUR_EXP_ID).links(link_id).exp_id == matched_exp_id
                        link_exists = 1;
                        break;
                    end
                end

                % if we didn't find a link then create the link between current
                % experience and the experience for the current visual template
                if link_exists == 0
                    EXPERIENCES(CUR_EXP_ID).numlinks = EXPERIENCES(CUR_EXP_ID).numlinks + 1;
                    EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).exp_id = matched_exp_id;
                    %  EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).d_xy = sqrt(ACCUM_DELTA_X^2 + ACCUM_DELTA_Y^2 + ACCUM_DELTA_Z^2);
                    EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).d_xy = sqrt(ACCUM_DELTA_X^2 + ACCUM_DELTA_Y^2);
                    EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).d_z = ACCUM_DELTA_Z;

                    EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).heading_yaw_exp_rad = ...
                        get_signed_delta_radian(EXPERIENCES(CUR_EXP_ID).yaw_exp_rad, atan2(ACCUM_DELTA_Y, ACCUM_DELTA_X)); % heading is the delta angle between current pose of exp and previous pose of exp

                    EXPERIENCES(CUR_EXP_ID).links(EXPERIENCES(CUR_EXP_ID).numlinks).facing_yaw_exp_rad = ... % facing is the direction of each exp
                        get_signed_delta_radian(EXPERIENCES(CUR_EXP_ID).yaw_exp_rad, ACCUM_DELTA_YAW);

                 end

            end

            % if there wasn't an experience with the current visual template and grid cell (x y z) and head direction cell (yaw, height)
            % then create a new experience
            if matched_exp_id == 0
                NUM_EXPS = NUM_EXPS + 1;
                create_new_exp(CUR_EXP_ID, NUM_EXPS, vt_id, xGc, yGc, zGc, curYawHdc, curHeight);
                matched_exp_id = NUM_EXPS;
            end

            PREV_EXP_ID = CUR_EXP_ID;
            CUR_EXP_ID = matched_exp_id;

            ACCUM_DELTA_X = 0;
            ACCUM_DELTA_Y = 0;
            ACCUM_DELTA_Z = 0;
            ACCUM_DELTA_YAW = EXPERIENCES(CUR_EXP_ID).yaw_exp_rad;
%             ACCUM_DELTA_HEIGHT = EXPERIENCES(CUR_EXP_ID).height_hdc;
        end

    end

    global EXP_CORRECTION;
    global EXP_LOOPS;

    % do the experience map correction interatively for all the links in all the experiences
    for i = 1:EXP_LOOPS

        for exp_id = 1:NUM_EXPS

            for link_id=1:EXPERIENCES(exp_id).numlinks

                % experience 0 has a link to experience 1
                e0 = exp_id;
                e1 = EXPERIENCES(exp_id).links(link_id).exp_id;

                % work out where e0 thinks e1 (x,y) should be based on the stored link information
                lx = EXPERIENCES(e0).x_exp + EXPERIENCES(e0).links(link_id).d_xy * cos(EXPERIENCES(e0).yaw_exp_rad + EXPERIENCES(e0).links(link_id).heading_yaw_exp_rad);
                ly = EXPERIENCES(e0).y_exp + EXPERIENCES(e0).links(link_id).d_xy * sin(EXPERIENCES(e0).yaw_exp_rad + EXPERIENCES(e0).links(link_id).heading_yaw_exp_rad);
                lz = EXPERIENCES(e0).z_exp + EXPERIENCES(e0).links(link_id).d_z;  %

                % correct e0 and e1 (x,y) by equal but opposite amounts
                % a 0.5 correction parameter means that e0 and e1 will be fully
                % corrected based on e0's link information
                EXPERIENCES(e0).x_exp = EXPERIENCES(e0).x_exp + (EXPERIENCES(e1).x_exp - lx) * EXP_CORRECTION;
                EXPERIENCES(e0).y_exp = EXPERIENCES(e0).y_exp + (EXPERIENCES(e1).y_exp - ly) * EXP_CORRECTION;
                EXPERIENCES(e0).z_exp = EXPERIENCES(e0).z_exp + (EXPERIENCES(e1).z_exp - lz) * EXP_CORRECTION;

                EXPERIENCES(e1).x_exp = EXPERIENCES(e1).x_exp - (EXPERIENCES(e1).x_exp - lx) * EXP_CORRECTION;
                EXPERIENCES(e1).y_exp = EXPERIENCES(e1).y_exp - (EXPERIENCES(e1).y_exp - ly) * EXP_CORRECTION;
                EXPERIENCES(e1).z_exp = EXPERIENCES(e1).z_exp - (EXPERIENCES(e1).z_exp - lz) * EXP_CORRECTION;

                % determine the angle between where e0 thinks e1's facing
                % should be based on the link information
                TempDeltaYawFacing = get_signed_delta_radian((EXPERIENCES(e0).yaw_exp_rad + EXPERIENCES(e0).links(link_id).facing_yaw_exp_rad), EXPERIENCES(e1).yaw_exp_rad);

                % correct e0 and e1 facing by equal but opposite amounts
                % a 0.5 correction parameter means that e0 and e1 will be fully
                % corrected based on e0's link information
                EXPERIENCES(e0).yaw_exp_rad = clip_radian_180(EXPERIENCES(e0).yaw_exp_rad + TempDeltaYawFacing * EXP_CORRECTION);
                EXPERIENCES(e1).yaw_exp_rad = clip_radian_180(EXPERIENCES(e1).yaw_exp_rad - TempDeltaYawFacing * EXP_CORRECTION);

            end
        end

    end

    % keep a frame by frame history of which experience was currently active
    global EXP_HISTORY;
    EXP_HISTORY = [EXP_HISTORY; CUR_EXP_ID];
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [idx,weights] = findconvexcone(x0,xs)
%FINDCONVEXCONE selects up to 3 points from x0 with xs in their conic span
%
%   Usage: [idx,weights] = findconvexcone(x0,xs)
%
%   Input parameters:
%       x0          - point cloud on a sphere around the origin / m [nx3]
%       xs          - desired direction as point in R^3 / m [1x3]
%
%   Output parameters:
%       idx         - row indices of N points in x0 [Nx1]
%                     where N is 1,2 or 3
%       weights     - weights [Nx1]
%
%   FINDCONVEXCONE(x0,xs) returns 1,2 or 3 row indices into x0 and non-negative
%   weights w1, ..., w3 such that w1*x1 + w2*x2 + w3*x3 with
%   [x1; x2; x3] == x0(idx,:) composes the point inside the triangle spanned
%   by x1, x2, x3.
%
%   x1...x3 are selected from the convex hull in R3.
%   Various precautions are taken to make this well-behaved in most cases.
%
%   (If all x0 and xs have unit norm this is VBAP.)
%
%   See also: findnearestneighbour, test_interpolation_point_idx

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);

%% ===== Prepare Grid (see local functions below) ========================
% Normalise x0 and xs to lie on unit sphere as only direction is relevant
xs = xs./norm(xs,2);
radii = vector_norm(x0,2);
if abs(max(radii) - min(radii)) >1e-3
     warning('%s: Grid is apparently not a sphere.', upper(mfilename))
end
x0 = x0./repmat(radii,[1,size(x0,2)]);

% Rotate to principal axes to enable 2D arrays
[x0, xs] = rotate_to_principal_axes(x0, xs);

% Calculate dummy points to enable "partial" arrays
dummy_points = augment_bounding_box(x0);
if ~isempty(dummy_points)
    dummy_indices = (1:size(dummy_points,1)) + size(x0,1);
    x0 = [x0; dummy_points];
end

%% ===== Computation =====================================================
% Delaunay triangulation of convex hull
simplices = convhulln(x0);

% Find x0 with smallest angle to xs
[~,most_aligned_point] = ...
    max(vector_product(x0,repmat(xs,size(x0,1),1),2));

% The simplices at "most aligned point" are the most likely candidates,
% put them at the beginning of the list
mask = logical(sum(simplices==most_aligned_point,2));
simplices = [simplices(mask,:); simplices(~mask,:) ];

% One of these simplices contains xs
for n = 1:size(simplices,1);
    A = x0(simplices(n,:),:);
    weights = xs/A;
    weights(abs(weights)<1e-10) = 0;
    if all(weights >= 0) % non-negative weights == convex combination
        idx = simplices(n,:);
        break;
    end
end
assert(all(weights >= 0), '%s: Negative weights. Shall never happen.', ...
    upper(mfilename))

if ~isempty(dummy_points)
    % Remove possible dummies from selected points
    dummy_mask = (dummy_indices == idx);
    if any(dummy_mask)
        idx(dummy_mask) = [];
        if ~weights(dummy_mask)==0
            warning('%s: Requested point lies outside grid.', upper(mfilename))
        end
        weights(dummy_mask) = [];
    end
end

% Normalise weights
weights = weights/sum(weights);

[weights,order] = sort(weights.','descend');
idx = idx(order).';
end

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

function [x0, xs] = rotate_to_principal_axes(x0, xs, gamma)
%ROTATE_TO_PRINCIPAL_AXES rotates x0 and xs to x0's principal axes.
% If the ratio of second-to-smallest singular values is < gamma,
% the third dimension is discarded.
%
%   Input parameters:
%       x0          - point cloud in R^3
%       xs          - point in R^3
%       gamma       - scalar in 0 < gamma << 1
%
%   Output parameters:
%       x0          - point cloud in R^3 or R^2
%       xs          - point in R^3 or R^2
if nargin < 3
    gamma = 0.1; % inverse of aspect ratio of principal axes
end

[~,S,V] = svd(x0);
x0 = x0*V;
xs = xs*V;
S = diag(S);
if S(end)/S(end-1) < gamma
    x0(:,3) = [];
    xs(3) = [];
    warning('SFS:findconvexcone','%s: Grid is apparently two-dimensional. ', ...
        upper(mfilename));
end
end

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

function dummy_points = augment_bounding_box(x0)
%AUGMENT_BOUNDING_BOX yields dummy points such that the origin is
% contained in the cartesian bounding box of x0.
dummy_points = -diag(sign(max(x0)) + sign(min(x0)));
dummy_points(~any(dummy_points,2),:) = [];
end

author: sfstoolbox

repo: sfs-matlab


text:
function volume = sphere_imp_zone_volume_3d ( r, h1, h2 )

%*****************************************************************************80
%
%% SPHERE_IMP_ZONE_VOLUME_3D computes the volume of a spherical zone in 3D.
%
%  Discussion:
%
%    An implicit sphere in 3D satisfies the equation:
%
%      sum ( ( P(1:DIM_NUM) - CENTER(1:DIM_NUM) )**2 ) = R**2
%
%    Draw any radius of the sphere and note the point P where the radius
%    intersects the sphere.  Now choose two points on the radius line, a
%    distance H1 and H2 from the point P.  Consider all the points on or within
%    the sphere whose projection onto the radius lies between these two points.
%    These points constitute the spherical zone, which can also be considered
%    the difference of two spherical caps.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real H1, H2, the distances that define the
%    thickness of the zone.  H1 and H2 must be between 0 and 2 * R.
%
%    Output, real VOLUME, the volume of the spherical zone
%
  h11 = min ( h1, h2 );
  h11 = max ( h11, 0.0 );

  if ( 2.0 * r <= h11 )
    volume = 0.0;
    return
  end

  h22 = max ( h1, h2 );
  h22 = min ( h22, 2.0 * r );

  if ( h22 <= 0.0 )
    volume = 0.0;
    return
  end

  volume = ( 1.0 / 3.0 ) * pi * ( ...
      h22 * h22 * ( 3.0 * r - h22 ) ...
    - h11 * h11 * ( 3.0 * r - h11 ) );

  return
end