text:
 function proj = cuboid_proj(cg, params, varargin)
%function proj = cuboid_proj(cg, params, varargin)
%|
%| Compute set of 2d line-integral projection views of cuboids.
%| Works for both parallel-beam and cone-beam geometry.
%|
%| in
%|	cg			ct_geom()
%|	params [ne 9]		cuboid parameters:
%|			[x_center y_center z_center  x_radius y_radius z_radius
%|				xy_angle_degrees z_angle_degrees  amplitude]
%| options
%|	oversample		over-sampling factor for emulating "strips"
%|				(to account for finite detector size)
%|
%| out
%|	proj	[ns nt na]	projection views
%|
%| Yong Long, 2008-08-28, University of Michigan, adapted from ellipsoid_proj()
%| 2013-05-24 Jeff Fessler

if nargin == 1 && streq(cg, 'test'), cuboid_proj_test, return, end
if nargin < 2, ir_usage, end

arg.oversample = 1;
arg = vararg_pair(arg, varargin);

proj = cuboid_proj_do(params, cg.s, cg.t, cg.dt, cg.ar, cg.source_zs, ...
		cg.dso, cg.dod, cg.dfs, arg.oversample);

end % cuboid_proj()

% cuboid_proj_line1()
function [lxmin lxmax] = cuboid_proj_line1(rx, p1, e1)
tmp = (e1 == 0);
e1(tmp) = inf;
% bounds of l corresponding to rect(x/rx)
lxmin = (-rx/2 - p1) ./ e1;
lxmax = ( rx/2 - p1) ./ e1;
% re-arrange the bounds so that lxmin contains the minimum l values
% and lxmax contains the maximum l values
temp = lxmin;
lxmin = min(lxmin, lxmax);
lxmax = max(temp, lxmax);
% exclude points where e1=0 by setting lxmin = -Inf and lxmax = Inf
lxmin(tmp) = -inf;
lxmax(tmp) = inf;
end % cuboid_proj_line1()

% cuboid_proj_do()
function proj = cuboid_proj_do(params, ss, tt, dt, ...
		beta, ... % [radians]
		source_zs, dso, dod, dfs, oversample)

if size(params, 2) ~= 9, error '9 parameters per cuboid', end

if oversample > 1
	ds = ss(2) - ss(1);
%	dt = tt(2) - tt(1); % fails if nt=1 so pass dt instead
	if any(abs(diff(ss) / ds - 1) > 1e-6) ...
	|| any(abs(diff(tt) / dt - 1) > 1e-6)
		fail 'uniform spacing required for oversampling'
	end
	No = oversample;
	% determine new finer sampling positions
	ss = outer_sum([-(No-1):2:(No-1)]'/(2*No)*ds, ss(:)'); % [No ns]
	tt = outer_sum([-(No-1):2:(No-1)]'/(2*No)*dt, tt(:)'); % [No nt]
	proj = cuboid_proj_do(params, ss(:), tt(:), dt/No, ...
		beta, source_zs, dso, dod, dfs, 1);
	proj = downsample3(proj, [No No 1]);
return
end

% determine equivalent parallel-beam projection coordinates, at beta=0
ns = length(ss);
nt = length(tt);
[sss ttt] = ndgrid(ss, tt);

if isinf(dso) % parallel beam
	uu = sss;
	vv = ttt;
	azim0 = zeros(size(sss));
	polar = zeros(size(sss));

elseif isinf(dfs) % cone-beam with flat detector
	[uu vv azim0 polar] = ir_coord_cb_flat_to_par(sss, ttt, dso, dod);

elseif dfs == 0 % cone-beam with arc detector
	[uu vv azim0 polar] = ir_coord_cb_arc_to_par(sss, ttt, dso, dod);

else
	fail 'not done'
end

clear sss ttt

cpolar = cos(polar);
spolar = sin(polar);
proj = zeros(ns, nt, numel(beta));

% loop over cuboids
for ip = 1:size(params,1)
	par = params(ip,:);

	cx = par(1);	rx = par(4);
	cy = par(2);	ry = par(5);
	cz = par(3);	rz = par(6);
	xang = deg2rad(par(7)); % xy-plane rotation
	zang = deg2rad(par(8)); % z-plane rotation
	if zang, error 'z rotation not done', end
	val = par(9);

	for ib = 1:length(beta)
		az = beta(ib) + azim0;

		% shift property of 3D transform:
		cz_eff = cz - source_zs(ib); % center relative to source
		ushift = cx * cos(az) + cy * sin(az);
		vshift = (cx * sin(az) - cy * cos(az)) .* spolar + cz_eff * cpolar;

		az = az - xang;
		us = uu - ushift;
		vs = vv - vshift;
		p1 = us .* cos(az) + vs .* sin(az) .* spolar;
		p2 = us .* sin(az) - vs .* cos(az) .* spolar;
		p3 = vs .* cpolar;

		e1 = -sin(az) .* cpolar; % x = p1 + l*e1
		e2 = cos(az) .* cpolar; % y = p2 + l*e2
		e3 = spolar; % z = p3 + l*e3

		[lxmin lxmax] = cuboid_proj_line1(rx, p1, e1);
%		clear p1
		[lymin lymax] = cuboid_proj_line1(ry, p2, e2);
%		clear p2
		[lzmin lzmax] = cuboid_proj_line1(rz, p3, e3);
%		clear p3

		lmin = max(lxmin, lymin);
		lmin = max(lmin, lzmin); % lower bound for l

		lmax = min(lxmax, lymax);
		lmax = min(lmax, lzmax); % upper bound for l

		ll = max(lmax - lmin, 0); % intersection only when lmax > lmin

		% cases where e(k) = 0 (rays along axes)

		tmp = e1 == 0; % sin(az) = 0
		zero_e = (-rx/2 <= us(tmp)) & (us(tmp) <= rx/2);
		ll(tmp) = ll(tmp) .* zero_e;

		tmp = e2 == 0; % cos(az) = 0
		zero_e = (-ry/2 <= us(tmp)) & (us(tmp) <= ry/2);
		ll(tmp) = ll(tmp) .* zero_e;

		tmp = e3 == 0; % sin(polar) = 0
		zero_e = (-rz/2 <= vs(tmp)) & (vs(tmp) <= rz/2);
		ll(tmp) = ll(tmp) .* zero_e;

if 0 % old way.  todo cut.  and verify e1 and e2 cases above
		if (e3 == 0) % sin(polar) = 0; so line along z
			zero_e = (-rz/2 <= vs) & (vs <= rz/2);
			ll = ll .* zero_e;
		end

		if (e1 == 0) % sin(az) = 0;
			zero_e = (-rx/2 <= us) & (us <= rx/2);
			ll = ll .* zero_e;
		end

		if (e2 == 0) % cos(az) = 0;
			zero_e = (-ry/2 <= us) & (us <= ry/2);
			ll = ll .* zero_e;
		end
end

		proj(:,:,ib) = proj(:,:,ib) + val * ll;
		clear proj_i
	end
end

end % cuboid_proj_do()

% cuboid_proj_test()
% internal test routine
function cuboid_proj_test

if 1
	params = [30 20 5  200 100 20  45 0 10];
	fun_proj = @(cg) cuboid_proj(cg, params, 'oversample', 2); % analytical
	fun_im = @(ig) cuboid_im(ig, params, 'oversample', 2); %, 'checkfov', true);

	ir_proj3_compare1(fun_proj, fun_im, 'chat', 1);
return % todo
end

down = 30;
cg = ct_geom('fan', 'ns', round(888/down), 'nt', 64, ...
	'na', 18, ...
	'ds', 1.0*down, 'dt', 1.1, ...
	'down', 1, ... % only downsample s and beta
	'dsd', 949, 'dod', 408, 'dfs', 0, ... % 3rd gen CT
	'dsd', inf, ... % parallel-beam
	'pitch', 0.5, ...
	'offset_t', 0.0, ...
	'offset_s', 0.25); % quarter detector
%	'dsd', 949, 'dod', 408, 'dfs', inf, ... % flat detector

ig = image_geom('nx', 512/2^2, 'nz', 64, 'fov', 512, 'dz', 0.625);
%cg.plot3(ig);

params = [30 20 5  200 100 20  45 0 10];
x = cuboid_im(ig, params, 'oversample', 2);

im plc 2 2
im(1, x), cbar

ya = cuboid_proj(cg, params, 'oversample', 2);
im(2, cg.s, cg.t, ya), cbar
titlef('analytical cone-beam projections, dfs=%g', cg.dfs)
xlabel s, ylabel t

%im clf, im(cg.s, cg.ad, permute(ya, [1 3 2])), cbar
%xlabel s, ylabel '\beta'

A = Gcone(cg, ig);
yd = A * x;

im(3, cg.s, cg.t, yd, 'discrete projections'), cbar

im(4, cg.s, cg.t, ya-yd), cbar
titlef('analytical - discrete')

max_percent_diff(ya, yd)

end % cuboid_proj_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function view = computeProbMap(view,scanList)
%
% view = computeProbMap(view,[scanList])
%
% Cycles through tSeries, computing the -log10 P values
% for each voxel. This is a first pass at a P map - not taking into account
% the non-gaussian properties of the spatio-temporal noise distribution.
%
% scanList:
%   0 - do all scans
%   number or list of numbers - do only those scans
%   default - prompt user via selectScans dialog
%
% Based on computeMeanMap
% and  Bantettini el al, 1993,
% Processing Strategies for Time Course Data Sets
% MRM 30:161-173 (1993) pp 171
%

nScans = numScans(view);
logProb=cell(1,nScans);
if strcmp(view.mapName,'logProbMap')
    % If exists, initialize to existing map
    map=view.map;
else
    % Otherwise, initialize to empty cell array
    map = cell(1,nScans);
end

% (Re-)set scanList
if ~exist('scanList','var')
    scanList = selectScans(view);
elseif scanList == 0
    scanList = 1:nScans;
end
if isempty(scanList)
    error('Analysis aborted');
end

% Compute it
waitHandle = mrvWaitbar(0,'Computing log10 P values from the tSeries.  Please wait...');
ncScans = length(scanList);
for iScan = 1:ncScans
    scan = scanList(iScan);
    nFrames=numFrames(view,scan);

    logProb{scanList(iScan)}=-(log10(computeCoherenceSignificance(view.co{scanList(iScan)},nFrames)));

    mrvWaitbar(scan/ncScans)
end
close(waitHandle);

% Set parameter map
view = setParameterMap(view,logProb,'log10Prob');

% Save file
saveParameterMap(view);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ins=ins_init(opt,avp0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INS initilization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sample date
ins.nt = 1/opt.sample_rate;
ins.old_imud = zeros(1,6);

% initialize position,attitude and velocity
ins.att = avp0(1:3);
ins.vel = avp0(4:6);
ins.pos = avp0(7:9);
ins.acc = zeros(3,1);
ins.Cnb = att2Cnb(ins.att);

% initialize imu error
ins.bg = zeros(3,1);
ins.ba = zeros(3,1);
ins.Kg = eye(3);
ins.Ka = eye(3);
ins.tauG = [inf;inf;inf];
ins.tauA = [inf;inf;inf];

% initialize earth related parameters
ins.eth = earth_update(ins.pos,ins.vel);
ins.Mpv = [0, ins.eth.Mpv4, 0; ins.eth.Mpv2, 0, 0; 0, 0, 1];
ins.wib = zeros(3,1);
ins.fb  = zeros(3,1);
ins.fn  = -ins.eth.gn;
ins.web = zeros(3,1);

% initialize state, state covariance matrix, system noise matrix
% state transition matrix, lever arm
init_att_unc = opt.init_att_unc;
init_vel_unc = opt.init_vel_unc;
init_pos_unc = opt.init_pos_unc;
init_bg_unc  = repmat(opt.init_bg_unc,1,3);
init_ba_unc  = repmat(opt.init_ba_unc,1,3);
psd_gyro = repmat(opt.psd_gyro,1,3);
psd_acce = repmat(opt.psd_acce,1,3);
psd_bg   = repmat(opt.psd_bg,1,3);
psd_ba   = repmat(opt.psd_ba,1,3);

ins.x  = [ins.att;ins.vel;ins.pos;ins.bg;ins.ba];
ins.P  = diag([init_att_unc,init_vel_unc,init_pos_unc,init_bg_unc,init_ba_unc].^2);
ins.Q  = diag([psd_gyro, psd_acce, zeros(1,3), psd_bg, psd_ba])*ins.nt;
ins.Phi= update_trans_mat(ins);
ins.lever = opt.lever';

ins.xa = zeros(15,1);
ins.Pa = zeros(15,15);

return

repo: GINav


text:
function x = rotatearound(x0,theta,mu)

[n1,n2] = size(x0);
[m1,m2] = size(mu);
theta = -theta;
R = [cos(theta),  sin(theta),
     -sin(theta), cos(theta)];
didflip = false;
if n1 ~= 2 & n2 == 2,
  x0 = x0';
  n2 = n1;
  didflip = true;
end;
if m1 ~= 2 & m2 == 2,
  mu = mu';
  m2 = m1;
end;

if m2 ~= n2,
  mu = repmat(mu(:,1),[1,n2]);
end;

x = R * (x0 - mu) + mu;
if didflip,
  x = x';
end

author: kristinbranson

repo: JAABA


text:
function install_blendenpik
% install_blendenpik
%
% Compiles and install Blendenpik.
% Blenenpik is a library for fast solution of dense least squares problem.
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.

%% Get configuration
[s, hostname] = system('hostname');
hostname = hostname(1:end-1);
config_file = sprintf('config/config.%s', hostname);

% Generate config file?
if (exist(config_file) ~= 0)
    regen_config = input('Configuration file already exist. Regenerate? ', 's');
    if (strcmp(lower(regen_config), 'yes'))
        configure_blenenpik;
    end
else
    configure_blenenpik;
end

% Read config file
fid = fopen(config_file);
if (fid == -1)
    error(sprintf('Config file %s not found!', CONFIG_FILE));
end
while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end
    eval(tline);
end

WISDOM_FILE = sprintf('%s/dat/wisdom_file.%s.%d.dat', pwd, hostname, FFTW_PLANNING_LEVEL);
if (ispc)
    WISDOM_FILE(WISDOM_FILE == '\') = '/';
end

%% Write config.h
fid = fopen('config.h', 'w');
fprintf(fid, '#define FFTW_QUANT %d\n', FFTW_QUANT);
fprintf(fid, '#define FFTW_TIMES %d\n', FFTW_TIMES);
fprintf(fid, '#define FFTW_WISDOM_FILE "%s"\n', WISDOM_FILE);
if (FFTW_PLANNING_LEVEL > 0)
    wisdom_flag = 'FFTW_WISDOM_ONLY';
else
    wisdom_flag = 'FFTW_ESTIMATE';
end
fprintf(fid, '#define FFTW_WISDOM_FLAG %s\n', wisdom_flag);
if (BLAS_UNDERSCORE)
    fprintf(fid, '#define BLAS_UNDERSCORE\n');
end

if (USE_FFTW)
    fprintf(fid, '#define USE_FFTW\n');
end

if (USE_SWHT)
    fprintf(fid, '#define USE_SWHT\n');
end

fclose(fid);

%% Compile
MEX_DIR = pwd;

if (~exist(MEX_DIR, 'dir'))
    mkdir(MEX_DIR);
end

incs = sprintf('-I"%s"', pwd);
libs = '';

if (USE_FFTW)
    incs = sprintf('%s -I"%s"', incs, FFTW_DIR_INC);
    libs = sprintf('%s -L"%s" -l%s', libs, FFTW_DIR_LIB, FFTW_NAME);
end

if (USE_SWHT)
    incs = sprintf('%s -I"%s"', incs, SWHT_DIR_INC);
    libs = sprintf('%s -L"%s" -l%s', libs, SWHT_DIR_LIB, SWHT_NAME);
end

prefix0 = sprintf('mex -outdir ''%s'' %s -O %s', MEX_DIR, incs, libs);
prefix = sprintf('%s -L"%s" -l%s -L"%s" -l%s %s ', prefix0,  LAPACK_DIR, LAPACK_NAME, BLAS_DIR, BLAS_NAME, ADDITIONAL_FLAGS);

eval([prefix '-c ' 'fftw_r2r.c']);
FFTW_R2R_OBJ = sprintf('"%s/fftw_r2r%s"', MEX_DIR, OBJ_EXT);
eval([prefix 'mex_fftw_r2r.c ' FFTW_R2R_OBJ]);
eval([prefix 'fast_unitary_transform.c ' FFTW_R2R_OBJ]);
eval([prefix 'fast_unitary_transform_size.c']);
eval([prefix 'build_fftw_wisdom.c']);
eval([prefix 'wtime.c']);
eval([prefix 'mex_dtrcon.c']);
eval([prefix 'mex_dgeqrf.c']);
eval([prefix 'mex_dormqr.c']);
eval([prefix 'mex_dtrsm.c']);
eval([prefix 'mex_dlange.c']);
eval([prefix 'lapack_solve_ls.c']);
eval([prefix 'dense_overdetermined_lsqr.c']);
eval([prefix 'dense_full_overdetermined_lsqr.c']);
eval([prefix 'dense_underdetermined_lsqr.c']);
addpath(MEX_DIR);
addpath(pwd);
savepath;

if (~exist('dat', 'dir'))
    mkdir('dat');
end

if (USE_FFTW)
    if (exist(WISDOM_FILE, 'file'))
        build_wisdom_text = input('FFTW wisdom already exist. Rebuild? (yes/no) ', 's');
        build_wisdom = strcmp(lower(build_wisdom_text), 'yes');
    else
        build_wisdom = true;
    end
    if (build_wisdom)
        disp('Building FFTW''s wisdom...');
        build_fftw_wisdom(FFTW_PLANNING_LEVEL, WISDOM_FILE);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function isbn = isbn_fill ( isbn )

%*****************************************************************************80
%
%% ISBN_FILL fills in a missing digit in an ISBN code.
%
%  Example:
%
%    Input:
%
%      0-8493-9?40-9
%
%    Output:
%
%      0-8493-9640-9
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Book Industry Study Group,
%    The Evolution in Product Identification:
%    Sunrise 2005 and the ISBN-13,
%    http://www.bisg.org/docs/The_Evolution_in_Product_ID.pdf
%
%  Parameters:
%
%    Input, string ISBN, a partial ISBN code, with a
%    single digit replaced by the character '?', signifying
%    that that digit is missing.
%
%    Output, string ISBN, the completed ISBN code.
%
  lenc = s_len_trim ( isbn );

  i = 0;
  isbn_pos = -1;
  digit_pos = -1;
  num_digit = 0;

  while ( 1 )

    i = i + 1;

    if ( lenc < i )
      break
    end

    c = isbn(i);

    if ( ch_is_digit ( c ) )

      num_digit = num_digit + 1;
      digit(num_digit) = isbn_to_i4 ( c );

    elseif ( ( num_digit == 9 & isbn(i:i) == 'X' ) | ...
             ( num_digit == 9 & isbn(i:i) == 'x' ) )

      num_digit = num_digit + 1;
      digit(num_digit) = isbn_to_i4 ( c );

    elseif ( c == '?' )

      if ( isbn_pos == -1 )

        num_digit = num_digit + 1;
        digit(num_digit) = 0;
        digit_pos = num_digit;
        isbn_pos = i;

      else
        fprintf ( 1, '\n' );
        fprintf ( 1, 'ISBN_FILL - Fatal error!\n' );
        fprintf ( 1, '  Only one question mark is allowed!\n' );
        error ( 'ISBN_FILL - Fatal error!' );
      end

    end

    if ( 10 <= num_digit )
      break
    end

  end

  if ( num_digit ~= 10 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ISBN_FILL - Fatal error!\n' );
    fprintf ( 1, '  The input ISBN code did not have 10 digits.\n' );
    error ( 'ISBN_FILL - Fatal error!' );
  end

  if ( isbn_pos == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ISBN_FILL - Fatal error!\n' );
    fprintf ( 1, '  A question mark is required!\n' );
    error ( 'ISBN_FILL - Fatal error!' );
  end

  check = 0;
  for i = 1 : 10
    check = check + ( 11 - i ) * digit(i);
  end

  check = mod ( check, 11 );

  if ( check == 0 )

    k = 0;
%
%  Need to solve the modular equation:
%
%    A * X = B mod C
%
%  Below is a stupid way.  One day I will come back and fix this up.
%
  else

    for i = 1 : 10
      j = ( 11 - digit_pos ) * i + check;
      if ( mod ( j, 11 ) == 0 )
        k = i;
      end
    end

  end

  isbn(isbn_pos) = i4_to_isbn ( k );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = r83p_sl ( n, a_lu, b, job, work2, work3, work4 )

%*****************************************************************************80
%
%% R83P_SL solves a R83P system.
%
%  Discussion:
%
%    The R83P storage format stores a periodic tridiagonal matrix as
%    a 3 by N array, in which each row corresponds to a diagonal, and
%    column locations are preserved.  The matrix value
%    A(1,N) is stored as the array entry A(3,N), and the matrix value
%    A(N,1) is stored as the array entry A(1,1).
%
%    The linear system must have been factored by R83P_FA.
%
%  Example:
%
%    Here is how a R83P matrix of order 5 would be stored:
%
%      A51 A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54 A15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 3.
%
%    Input, real A_LU(3,N), the LU factors from R83P_FA.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Input, integer JOB, specifies the system to solve.
%    0, solve A * x = b.
%    nonzero, solve A' * x = b.
%
%    Input, real WORK2(N-1), WORK3(N-1), WORK4, factor data from R83P_FA.
%
%    Output, real X(N), the solution to the linear system.
%
  x(1:n) = b(1:n);

  if ( job == 0 )
%
%  Solve A1 * X1 = B1.
%
    x(1:n-1) = r83_np_sl ( n-1, a_lu, x, job );
%
%  X2 = B2 - A3 * X1
%
    x(n) = x(n) - a_lu(1,1) * x(1) - a_lu(3,n-1) * x(n-1);
%
%  Solve A4 * X2 = X2
%
    x(n) = x(n) / work4;
%
%  X1 := X1 - inverse ( A1 ) * A2 * X2.
%
    x(1:n-1) = x(1:n-1) - work2(1:n-1) * x(n);

  else
%
%  Solve A1' * X1 = B1.
%
    x(1:n-1) = r83_np_sl ( n-1, a_lu, x, job );
%
%  X2 := X2 - A2' * B1
%
    x(n) = x(n) - a_lu(3,n) * x(1) - a_lu(1,n) * x(n-1);
%
%  Solve A4 * X2 = X2.
%
    x(n) = x(n) / work4;
%
%  X1 := X1 - transpose ( inverse ( A1 ) * A3 ) * X2.
%
    x(1:n-1) = x(1:n-1) - work3(1:n-1) * x(n);

  end

  return
end


text:
function [A,b,Aeq,beq] = row2gen(Ain,rl,ru)
%ROW2GEN Convert Linear A with Row Bounds to Linear Inequality & Equality Constraints
%   [A,b,Aeq,beq] = row2gen(A,rl,ru)

%   Copyright (C) 2012 Jonathan Currie (I2C2)

if(isempty(Ain))
    A = []; b = []; Aeq = []; beq = [];
    return;
elseif(isempty(ru))
    error('You must supply Ain, rl and ru!');
end
if(isempty(rl))
    rl = -Inf(size(ru)); %default
end
%Transpose as Required
if(size(rl,2) > 1)
    rl = rl';
end
if(size(ru,2) > 1)
    ru = ru';
end

%Indices
eq = rl == ru; neq = ~eq;
ile = isfinite(ru) & neq;
ige = isfinite(rl) & neq;

%Ineq
A = [ Ain(ile,:);
     -Ain(ige,:)];
b = [ ru(ile);
     -rl(ige)];
%Eq
if(any(eq))
    Aeq = Ain(eq,:); beq = ru(eq);
else
    Aeq = []; beq = [];
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function hdr = ReadEfileHeader(EfileName)
% header = ReadEfileHeader(EfileName)
%
% Returns header info from Efile produced by G. Glover's auto-recon program.
% Output is structure having field names corresponding to Glover's parameters.
%
% Ress 9/01
% $Author: sayres $
% $Date: 2005/08/23 08:06:25 $
if nargin==0
    if exist('Raw/Pfiles','dir')
        startDir = 'Raw/Pfiles';
    else
        startDir = pwd;
    end
    [f p] = myUiGetfile(startDir,'*.7','Choose an E-file');
    EfileName = fullfile(p,f);
end

hdr = [];

fid = fopen(EfileName);
if fid == -1
    Alert(['Could not open file: ' EfileName])
    return
end

nLine = 0;
while 1
    line = fgetl(fid);
    if ~ischar(line), break, end
    ieq = findstr(line, '=');
    lhs = line(1:ieq-2);
    rhs = line(ieq+2:end);
    switch lhs
    case 'rev'
        % Revision
        hdr.rev = rhs;
    case 'date of scan'
        % Date
        hdr.date = rhs;
        % Make this more reasonable
        hdr.date=[hdr.date(1:6),hdr.date(8:end)];
    case 'time of scan'
        % Time
        hdr.time = rhs;
    case 'patient name'
        % Name
        hdr.name = rhs;

    case 'psd'
        % PSD
        hdr.psd = rhs;
    case 'coil'
        % Coil
        hdr.coil = rhs;
    case 'slquant'
        % Number of slices
        hdr.slquant = getNumberFromString(rhs);
    case 'num time frames'
        % Number of frames
        hdr.nframes = getNumberFromString(rhs);
    case 'numextra discards'
        % Number of discarded shots
        hdr.nextra = getNumberFromString(rhs);
    case 'nshot'
        % Number of interleaves
        hdr.nshots = getNumberFromString(rhs);
    case 'FOV'
        % Field of view (mm)
        hdr.FOV = getNumberFromString(rhs);
    case 'slice thick'
        % Slice thickness (mm)
        hdr.sliceThickness = getNumberFromString(rhs);
    case 'skip'
        % Slice spacing (skip)
        hdr.skip = getNumberFromString(rhs);
    case 'TR'
        % Echo time, TR (ms)
        hdr.TR = getNumberFromString(rhs);
    case 'TE'
        % Repetition time, TE (ms)
        hdr.TE = getNumberFromString(rhs);
    case 'time/frame'
        % Acquisition time (ms)
        hdr.tAcq = getNumberFromString(rhs);
    case 'equiv matrix size'
        % Equivalent matrix size
        hdr.equivMatSize = getNumberFromString(rhs);
    case 'imgsize'
        % Image size
        hdr.imgsize = getNumberFromString(rhs);
    case 'pixel size'
        % Pixel size (mm)
        hdr.pixel = getNumberFromString(rhs);
    case 'freq'
        % Frequency (MHz)
        hdr.freq = getNumberFromString(rhs)/1.e6;
    case 'R1, R2, TG (mps)'
        % Gains [R1, R2, TG]
        hdr.R1 = getNumberFromString(rhs, 1);
        hdr.R2 = getNumberFromString(rhs, 2);
        hdr.TG = getNumberFromString(rhs, 3);
    otherwise
    end
end


text:
%  Figure 10.24      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%   fig10_24.m is a script to generate Fig. 10.24,
%   the transient response of the LQR symmetric rootlocus compensator
%   of the satellite position control, non-colocated case WITH ESTIMATOR

% parameter values
m=[1, 0.1]; k=[0, 0.091] ; d=[0, 0.0036]; k1=[0, 0.4];
% call function
[f,g,h,j] = twomass(m,k,d);
s=[f, g;h, 0];r=[0*g;1]; n=s\r;nx=n(1:4);nu=n(5);
% call function
[f1,g,h,j] = twomass(m,k1,d);

% form G(s)G(-s) model
a=[f, 0*f;
-h'*h, -f'];
b=[g;0*g];
d=[0];
c=[0*h, g'];
hold off; clf
P=eig(a-b*c*0.1621);
pc=P(real(P<0)==1);
K=place(f,g,pc);
nbar=nu+K*nx;
% eig(f-g*K)
P=eig(a-b*c*3.056e7);
pe=P(real(P<0)==1);
L=place(f',h',pe)';
ac=f-g*K-L*h ;bc=L;cc=K;dc=0;
[Aol,Bol,Col,Dol]=series(ac,bc,cc,dc,f,g,h,j);
[acl,bcl,ccl,dcl]=feedback(f,g,h,j,ac,bc,cc,dc);
[acl1,bcl,ccl,dcl]=feedback(f1,g,h,j,ac,bc,cc,dc);
bcl= nbar*[g;g];
t=0:.25:30;
syscl=ss(acl,bcl,ccl,dcl);
step(syscl,t);
hold on;
gtext('nominal case')
syscl1=ss(acl1,bcl,ccl,dcl);
step(syscl1,t) ;
gtext('stiff-spring case');
title('Fig. 10.24 Closed-loop step response for the SRL design with an estimator')
%grid
nicegrid


text:
function polylogarithm_values_test ( )

%*****************************************************************************80
%
%% POLYLOGARITHM_VALUES_TEST demonstrates the use of POLYLOGARITHM_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYLOGARITHM_VALUES_TEST:\n' );
  fprintf ( 1, '  POLYLOGARITHM_VALUES stores values of\n' );
  fprintf ( 1, '  the polylogarithm.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       N           Z           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, z, fx ] = polylogarithm_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %12f  %24.16f\n',  n, z, fx );

  end

  return
end


text:
function [EventsFUN,EventsValue]=sdeeventsfun(solver,tspan,y0,options)
%SDEEVENTSFUN  Process events function arguments for SDE solvers and functions.
%
%   See also:
%       SDEARGUMENTS, SDEARGUMENTS_PROCESS, SDEOUTPUTFUN, SDERANDFUN, SDEGET,
%       FUNCTION_HANDLE

%   Andrew D. Horchler, horchler @ gmail . com, Created 5-2-13
%   Revision: 1.2, 5-2-13

EventsFUN = sdeget(options,'EventsFUN',[],'flag');
if ~isempty(EventsFUN)
    if ~isa(EventsFUN,'function_handle')
        error('SDETools:sdeeventsfun:EventsFUNNotAFunctionHandle',...
              'EventsFUN, if specified, must be a function handle.  See %s.',...
              solver);
    end

    % Check output of EventsFUN at initial condition and save value
    try
        [EventsValue,isterminal,direction] = EventsFUN(tspan(1),y0);
    catch err
        switch err.identifier
            case 'MATLAB:TooManyInputs'
                error('SDETools:sdeeventsfun:EventsFUNTooFewInputs',...
                      'EventsFUN must have at least two inputs.  See %s.',...
                      solver);
            case 'MATLAB:TooManyOutputs'
                error('SDETools:sdeeventsfun:EventsFUNNoOutput',...
                     ['The output of EventsFUN was not specified. EventsFUN '...
                      'must return three non-empty vectors.  See %s.'],solver);
            case 'MATLAB:unassignedOutputs'
                error('SDETools:sdeeventsfun:EventsFUNUnassignedOutput',...
                     ['The first output of EventsFUN was not assigned.  '...
                      'See %s.'],solver);
            case 'MATLAB:minrhs'
                error('SDETools:sdeeventsfun:EventsFUNTooManyInputs',...
                     ['EventsFUN requires one or more input arguments '...
                      '(parameters) that were not supplied.  See %s.'],solver);
            otherwise
                rethrow(err);
        end
    end
    if ~isvector(EventsValue) || isempty(EventsValue)...
            || ~all(isfinite(EventsValue))
        error('SDETools:sdeeventsfun:InvalidEventsValue',...
             ['The first output of EventsFUN, ''Value'', must be a '...
              'non-empty finite vector.  See %s.'],solver);
    end
    if ~isvector(isterminal)...
            || ~any(length(isterminal) == [length(EventsValue) 1])
        error('SDETools:sdeeventsfun:EventsIsterminalDimensionMismatch',...
             ['The second output of EventsFUN, ''IsTerminal'', must be a '...
              'scalar or a vector the same length as the first output.  '...
              'See %s.'],solver);
    end
    if ~all(isterminal == 0 | isterminal == 1)
        error('SDETools:sdeeventsfun:InvalidEventsIsterminal',...
             ['The elements of the second output of EventsFUN, '...
              '''IsTerminal'', must be equal to 0 (false) or 1 (true).  '...
              'See %s.'],solver);
    end
    if ~isempty(direction)
        if ~isvector(direction)...
                || ~any(length(direction) == [length(EventsValue) 1])
            error('SDETools:sdeeventsfun:EventsDirectionDimensionMismatch',...
                 ['If the third output of EventsFUN, ''Direction'', is not '...
                  'specified as the empty matrix, [], it must be a scalar '...
                  'or a vector the same length as first output.  See %s.'],...
                  solver);
        end
        if ~all(direction == 0 | direction == 1 | direction == -1)
            error('SDETools:sdeeventsfun:InvalidEventsDirection',....
                 ['The third output of EventsFUN, ''Direction'', must be '...
                  'equal to 0 (default), 1, or -1.  See %s.'],solver);
        end
    end
else
    EventsValue = [];
end

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
function test_bug2221

% MEM 1gb
% WALLTIME 00:15:00
% DEPENDENCY ft_selectdata ft_selectdata_new ft_multiplotER

% this test function intends to test the functionality of:
%  - ft_selectdata_new with respect to multiple input: does it correctly
%  identifies the intersection between the multiple input dimensions?
%  - ft_selectdata_new (new functionality): take a union rather than
%  intersection, this is to be built in as an optional switch
%  - ft_multiplotER: if input data contains not entirely overlapping time
%  axes, use selectdata for getting the union and plotting should then work

% as of september 2019, ft_selectdata_new/old have become obsolete, and the
% explicit call of ft_selectdata_new has been replaced by a call to
% ft_selectdata

% create some data first
freq1 = [];
freq1.label = {'chan01';'chan02'};
freq1.freq  = (1:10);
freq1.time  = (1:15);
freq1.powspctrm = randn(2,10,15);
freq1.dimord = 'chan_freq_time';

freq2 = freq1;
freq2.powspctrm = randn(2,8,12);
freq2.label = {'chan02';'chan03'};
freq2.freq  = (5:12);
freq2.time  = (3:14);

cfg                  = [];
cfg.select           = 'intersect';
[freq1out, freq2out] = ft_selectdata(cfg,freq1,freq2);
assert(isequal(rmfield(freq1out, {'cfg','powspctrm'}),rmfield(freq2out, {'cfg','powspctrm'})));

cfg.select        = 'union';
[freq1out, freq2out] = ft_selectdata(cfg,freq1,freq2);

%%%%%%%%%%%
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2221.mat');
load(filename);

data1 = data_S10_cond1;
data2 = data_S10_cond2;
cfg = [];
cfg.select  = 'intersect';

[data1out, data2out] = ft_selectdata(cfg, data1, data2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [eef_transform,J]=directKinematics_1(q,TefTool,dh)
%% Calculates the direct kinematics of a robotic manipulator

% Syntax:
% [eef_transform,J]=directKinematics(q,TefTool,dh)

% Arreguments:
% q: angles of the joints of the robot.
% TefTool: transfomr matrix from the tool frame to the flange frame of the
% robot.
% dh: is a structure with the Denavit-Hartenberg parameters of the robot.

% Return value:
% eef_transform: transfomration matrix from end-effector to tool
% J: jacobean at the tool center point (TCP) of the end-effector EEF.

% Copyright:
% Mohammad SAFEEA
% 16th-Aug-2017
% updated: 22nd-June-2018

%% DH PARAMETERS FOR THE ROBOT
alfa=dh.alfa;
d=dh.d;
a=dh.a;

%% Calculating the direct Kinematics
T=zeros(4,4,7);
i=1;
T(:,:,i)=getDHMatrix(alfa{i},q(i),d{i},a{i});
    for i=2:7
        T(:,:,i)=T(:,:,i-1)*getDHMatrix(alfa{i},q(i),d{i},a{i});
        T(:,:,i)=normalizeColumns(T(:,:,i));
    end
        T(:,:,7)=T(:,:,7)*TefTool;
        T(:,:,7)=normalizeColumns(T(:,:,7));

    pef=T(1:3,4,7);
    for i=1:7
        k=T(1:3,3,i);
        pij=pef-T(1:3,4,i);
        J(1:3,i)=cross(k,pij);
        J(4:6,i)=k;
    end
%% End effector transform
eef_transform=T(:,:,7);
end

function T=getDHMatrix(alfa,theta,d,a)
T=zeros(4,4);

calpha=cos(alfa);
sinalpha=sin(alfa);
coshteta=cos(theta);
sintheta=sin(theta);

T(1,1)=coshteta;
T(2,1)=sintheta*calpha;
T(3,1)=sintheta*sinalpha;

T(1,2)=-sintheta;
T(2,2)=coshteta*calpha;
T(3,2)=coshteta*sinalpha;

T(2,3)=-sinalpha;
T(3,3)=calpha;

T(1,4)=a;
T(2,4)=-sinalpha*d;
T(3,4)=calpha*d;
T(4,4)=1;

end

function normalizedT=normalizeColumns(T)
%% This function is used to normalize the columns of a rotation matrix with
% some numerical errors resulting from matrix multiplication problems
r=zeros(4,3); % corrected rotatio matrix, with zero badding row at the end
for j=1:3
    r(1:3,j)=T(1:3,j)/norm(T(1:3,j));
end
normalizedT=[r,T(:,4)];
end

repo: KST-Kuka-Sunrise-Toolbox


text:
dataset_path = '/home/zhixuan/data1/HUMBI/HUMBI_uploaded/Hand_81_140';
subject = 82;
frame = 1;
side = 'right'; % 'left' or 'right'
showAxis = false; % whether show 3D axis when plotting
showKpsIdx = false; % whether plot indices of keypoints
vis3D = true; % if true, visualize 3D mesh and kps
               % if false, visualize their reprojections on image
%%%%%%%%%%%% modify above lines accordingly %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
config;

% extract camera params
intrinsics_path = sprintf(intrinsics_path_format, subject);
extrinsics_path = sprintf(extrinsics_path_format, subject);
[M, C, ~, ~] = extractCameraParameters(intrinsics_path, extrinsics_path, num_cam);
[ground_plane_normal, ground_plane_center] = getGroundPlaneFromCameras(C(:, 2:2:66), cameras_to_plane_offset);

% read data
reconstruction_dir = sprintf(recon_dir_format, subject, frame);
[vertices, kps, mano_params] = read_hand_recon(reconstruction_dir, side);

% visualization
set(gcf, 'OuterPosition', get(0, 'Screensize')); % maximize the figure
if vis3D  % visualize 3D mesh and kps
    subplot(121); setup_vis; camlight('headlight');
    h = trimesh(tri_list, vertices(:,1), vertices(:,2), vertices(:,3), ...
        'FaceColor', [0 150 230]/255, 'EdgeColor', 'none', 'LineWidth', 0.1);
    h.FaceLighting = 'flat'; % or 'gouraud'
    subplot(122); setup_vis;
    drawHand(kps, showKpsIdx);
    sgtitle(sprintf('subject %d, frame %d', subject, frame));
else
    bboxes = readmatrix([sprintf(img_dir_format, subject, frame), '/list.txt']); % n x 5, bboxes on original image
    for i = 1:size(bboxes, 1)
        % compute reprojections
        cam = bboxes(i, 1);
        bbox = bboxes(i, 2:end); % 1 x 4, [xmin, xmax, ymin, ymax]
        scale_x = (bbox(2) - bbox(1) + 1) / img_size(2);
        scale_y = (bbox(4) - bbox(3) + 1) / img_size(1);
        % Notice: scale_x and scale_y may not be exactly the same but should be very close
        rvertices = reproject(vertices, M(:, :, cam+1)); % 3448 x 2
        rvertices = (rvertices - [bbox(1), bbox(3)]) ./ [scale_x, scale_y];
        rkps = reproject(kps, M(:, :, cam+1)); % 66 x 2
        rkps = (rkps - [bbox(1), bbox(3)]) ./ [scale_x, scale_y];

        % draw
        img_path = sprintf([img_dir_format, '/image%07d.png'], subject, frame, cam);
        img = imread(img_path);
        subplot(121); imshow(img); hold on
        scatter(rvertices(:,1), rvertices(:,2), 100, '.', 'MarkerEdgeColor', [0 150 230]/255, 'MarkerEdgeAlpha', 0.8);
        subplot(122); imshow(img); hold on
        drawHand(rkps, showKpsIdx);
        sgtitle(sprintf('subject %d, frame %d, cam %d (press "Enter" to go to next view)', subject, frame, cam));
        pause;
    end
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
function c = l1nn_cholesky ( n, h )

%*****************************************************************************80
%
%% L1NN_CHOLESKY computes the Cholesky factor of the 1D NN Laplacian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real C(N,N), the Cholesky factor.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1NN_CHOLESKY - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1NN_CHOLESKY - Fatal error!' );
  end

  c = zeros ( n, n );

  for i = 1 : n - 1
    c(i,i)   = + 1.0;
    c(i,i+1) = - 1.0;
  end

  c(1:n,1:n) = c(1:n,1:n) / h;

  return
end


text:
function [pred_boxes, scores, box_deltas_, anchors_, scores_] = proposal_im_detect(conf, caffe_net, im)
% [pred_boxes, scores, box_deltas_, anchors_, scores_] = proposal_im_detect(conf, im, net_idx)
% --------------------------------------------------------
% Faster R-CNN
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    im = single(im);
    [im_blob, im_scales] = get_image_blob(conf, im);
    im_size = size(im);
    scaled_im_size = round(im_size * im_scales);

    % permute data into caffe c++ memory, thus [num, channels, height, width]
    im_blob = im_blob(:, :, [3, 2, 1], :); % from rgb to brg
    im_blob = permute(im_blob, [2, 1, 3, 4]);
    im_blob = single(im_blob);

    net_inputs = {im_blob};

    % Reshape net's input blobs
    caffe_net.reshape_as_input(net_inputs);
    output_blobs = caffe_net.forward(net_inputs);

    % Apply bounding-box regression deltas
    box_deltas = output_blobs{1};
    featuremap_size = [size(box_deltas, 2), size(box_deltas, 1)];
    % permute from [width, height, channel] to [channel, height, width], where channel is the
        % fastest dimension
    box_deltas = permute(box_deltas, [3, 2, 1]);
    box_deltas = reshape(box_deltas, 4, [])';

    anchors = proposal_locate_anchors(conf, size(im), conf.test_scales, featuremap_size);
    pred_boxes = fast_rcnn_bbox_transform_inv(anchors, box_deltas);
      % scale back
    pred_boxes = bsxfun(@times, pred_boxes - 1, ...
        ([im_size(2), im_size(1), im_size(2), im_size(1)] - 1) ./ ([scaled_im_size(2), scaled_im_size(1), scaled_im_size(2), scaled_im_size(1)] - 1)) + 1;
    pred_boxes = clip_boxes(pred_boxes, size(im, 2), size(im, 1));

    assert(conf.test_binary == false);
    % use softmax estimated probabilities
    scores = output_blobs{2}(:, :, end);
    scores = reshape(scores, size(output_blobs{1}, 1), size(output_blobs{1}, 2), []);
    % permute from [width, height, channel] to [channel, height, width], where channel is the
        % fastest dimension
    scores = permute(scores, [3, 2, 1]);
    scores = scores(:);

    box_deltas_ = box_deltas;
    anchors_ = anchors;
    scores_ = scores;

    if conf.test_drop_boxes_runoff_image
        contained_in_image = is_contain_in_image(anchors, round(size(im) * im_scales));
        pred_boxes = pred_boxes(contained_in_image, :);
        scores = scores(contained_in_image, :);
    end

    % drop too small boxes
    [pred_boxes, scores] = filter_boxes(conf.test_min_box_size, pred_boxes, scores);

    % sort
    [scores, scores_ind] = sort(scores, 'descend');
    pred_boxes = pred_boxes(scores_ind, :);
end

function [data_blob, rois_blob, im_scale_factors] = get_blobs(conf, im, rois)
    [data_blob, im_scale_factors] = get_image_blob(conf, im);
    rois_blob = get_rois_blob(conf, rois, im_scale_factors);
end

function [blob, im_scales] = get_image_blob(conf, im)
    if length(conf.test_scales) == 1
        [blob, im_scales] = prep_im_for_blob(im, conf.image_means, conf.test_scales, conf.test_max_size);
    else
        [ims, im_scales] = arrayfun(@(x) prep_im_for_blob(im, conf.image_means, x, conf.test_max_size), conf.test_scales, 'UniformOutput', false);
        im_scales = cell2mat(im_scales);
        blob = im_list_to_blob(ims);
    end
end

function [rois_blob] = get_rois_blob(conf, im_rois, im_scale_factors)
    [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, im_scale_factors);
    rois_blob = single([levels, feat_rois]);
end

function [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, scales)
    im_rois = single(im_rois);

    if length(scales) > 1
        widths = im_rois(:, 3) - im_rois(:, 1) + 1;
        heights = im_rois(:, 4) - im_rois(:, 2) + 1;

        areas = widths .* heights;
        scaled_areas = bsxfun(@times, areas(:), scales(:)'.^2);
        levels = max(abs(scaled_areas - 224.^2), 2);
    else
        levels = ones(size(im_rois, 1), 1);
    end

    feat_rois = round(bsxfun(@times, im_rois-1, scales(levels)) / conf.feat_stride) + 1;
end

function [boxes, scores] = filter_boxes(min_box_size, boxes, scores)
    widths = boxes(:, 3) - boxes(:, 1) + 1;
    heights = boxes(:, 4) - boxes(:, 2) + 1;

    valid_ind = widths >= min_box_size & heights >= min_box_size;
    boxes = boxes(valid_ind, :);
    scores = scores(valid_ind, :);
end

function boxes = clip_boxes(boxes, im_width, im_height)
    % x1 >= 1 & <= im_width
    boxes(:, 1:4:end) = max(min(boxes(:, 1:4:end), im_width), 1);
    % y1 >= 1 & <= im_height
    boxes(:, 2:4:end) = max(min(boxes(:, 2:4:end), im_height), 1);
    % x2 >= 1 & <= im_width
    boxes(:, 3:4:end) = max(min(boxes(:, 3:4:end), im_width), 1);
    % y2 >= 1 & <= im_height
    boxes(:, 4:4:end) = max(min(boxes(:, 4:4:end), im_height), 1);
end

function contained = is_contain_in_image(boxes, im_size)
    contained = boxes >= 1 & bsxfun(@le, boxes, [im_size(2), im_size(1), im_size(2), im_size(1)]);

    contained = all(contained, 2);
end

author: ShaoqingRen


text:
function [lat,lon] = importPlates(filename, startRow, endRow)
%IMPORTPLATES Import numeric data from a text file as column vectors.
%   [LAT,LON] = IMPORTPLATES(FILENAME) Reads data from text file FILENAME for
%   the default selection.
%
%   [LAT,LON] = IMPORTPLATES(FILENAME, STARTROW, ENDROW) Reads data from rows
%   STARTROW through ENDROW of text file FILENAME.
%
% Example:
%   [lat,lon] = importPlates('All_boundaries.txt',2, 12377);
%
%    See also TEXTSCAN.

% Auto-generated by MATLAB on 2015/01/23 11:51:24

%% Initialize variables.
if nargin<=2
    startRow = 2;
    endRow = inf;
end

%% Read columns of data as strings:
% For more information, see the TEXTSCAN documentation.
formatSpec = '%9s%s%[^\n\r]';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', '', 'WhiteSpace', '', 'HeaderLines', startRow(1)-1, 'ReturnOnError', false);
for block=2:length(startRow)
    frewind(fileID);
    dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', '', 'WhiteSpace', '', 'HeaderLines', startRow(block)-1, 'ReturnOnError', false);
    for col=1:length(dataArray)
        dataArray{col} = [dataArray{col};dataArrayBlock{col}];
    end
end

%% Close the text file.
fclose(fileID);

%% Convert the contents of columns containing numeric strings to numbers.
% Replace non-numeric strings with NaN.
raw = repmat({''},length(dataArray{1}),length(dataArray)-1);
for col=1:length(dataArray)-1
    raw(1:length(dataArray{col}),col) = dataArray{col};
end
numericData = NaN(size(dataArray{1},1),size(dataArray,2));

for col=[1,2]
    % Converts strings in the input cell array to numbers. Replaced non-numeric
    % strings with NaN.
    rawData = dataArray{col};
    for row=1:size(rawData, 1);
        % Create a regular expression to detect and remove non-numeric prefixes and
        % suffixes.
        regexstr = '(?<prefix>.*?)(?<numbers>([-]*(\d+[\,]*)+[\.]{0,1}\d*[eEdD]{0,1}[-+]*\d*[i]{0,1})|([-]*(\d+[\,]*)*[\.]{1,1}\d+[eEdD]{0,1}[-+]*\d*[i]{0,1}))(?<suffix>.*)';
        try
            result = regexp(rawData{row}, regexstr, 'names');
            numbers = result.numbers;

            % Detected commas in non-thousand locations.
            invalidThousandsSeparator = false;
            if any(numbers==',');
                thousandsRegExp = '^\d+?(\,\d{3})*\.{0,1}\d*$';
                if isempty(regexp(thousandsRegExp, ',', 'once'));
                    numbers = NaN;
                    invalidThousandsSeparator = true;
                end
            end
            % Convert numeric strings to numbers.
            if ~invalidThousandsSeparator;
                numbers = textscan(strrep(numbers, ',', ''), '%f');
                numericData(row, col) = numbers{1};
                raw{row, col} = numbers{1};
            end
        catch me
        end
    end
end

%% Replace non-numeric cells with NaN
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),raw); % Find non-numeric cells
raw(R) = {NaN}; % Replace non-numeric cells

%% Allocate imported array to column variable names
lat = cell2mat(raw(:, 1));
lon = cell2mat(raw(:, 2));

author: geoscience-community-codes

repo: GISMO


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

%%
clear; close all; clc;

%% Syntax
% |[FLs,VLs,logicAttachFaces,indCurve,Eb]=localRefineMap(FL,VL,V_SLIL_Lunate,distGrow,numGrowStepsRefine,cParSmoothCurve);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function prob = LogisticRegressionTest(data, model)

prob = 1 ./ (1 + exp(- (model.w(1:end-1)' * data.feat) - model.w(end)));
% prob = 1 ./ (1 + exp(- (model.w' * data.feat)));

author: flyers

repo: drone-tracking


text:
[out, th_value] = threshold(rcc,'isodata',1);
siz = size(out);
% out=fliplr(out);
a=zeros(siz(1:2));
num=1;
for ind=1:siz(1)*siz(2)
    [i,j] = ind2sub(siz(1:2),ind);
    if a(ind) == 0
        indexlabel = intersect(find(out(:,:,ind-1)>0),find(a==0));
        if ~isempty(indexlabel)
            a([indexlabel; ind]) = num;

            num=num+1;
        end
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Example file for ba_interp3

mex -O ba_interp3.cpp
%% Interpolation of 3D volumes (e.g. distance transforms)

% Create Low Resolution Volume
sz=16;
[dist.x dist.y dist.z] = meshgrid(linspace(-1,1,sz), linspace(-1,1,sz), linspace(-1,1,sz));
dist.D = min(sqrt((dist.x - 0).^2 + (dist.y - 0.35).^2 + (dist.z - 0).^2) - 0.5, ...
                   sqrt((dist.x - 0).^2 + (dist.y + 0.35).^2 + (dist.z - 0).^2) - 0.5);

interpolators={@interp3, @ba_interp3};
methods = {'nearest', 'linear', 'cubic'};
sizes=[20:10:100];
samples=20;

for interpolator_i=1:numel(interpolators)
  interp_fun = interpolators{interpolator_i};
  figure(interpolator_i);

  subplot 221;
  cla
  hx = slice(dist.x, dist.y, dist.z, dist.D, 0.2, [], []);
  set(hx,'EdgeColor','none');
  p = patch(isosurface(dist.x,dist.y,dist.z,dist.D, 0));
  isonormals(dist.x,dist.y,dist.z,dist.D,p)
  set(p,'FaceColor','red','EdgeColor','none');
  daspect([1 1 1])
  view(3);
  axis([-1 1 -1 1 -1 1]);
  camlight
  lighting gouraud
  title('Low resolution distance image');

  %% Increase resolution by interpolation
  for size_i=1:numel(sizes)
    sz1 = sizes(size_i);
    [dist_h.x dist_h.y dist_h.z] = meshgrid(linspace(-1,1,sz1), linspace(-1,1,sz1), linspace(-1,1,sz1));
    for i=1:numel(methods)
      for sample=1:samples
        tic;
        dist_h.D = interp_fun(dist.x, dist.y, dist.z, dist.D, dist_h.x, dist_h.y, dist_h.z, methods{i});
        took(interpolator_i,size_i,i,sample)=toc;
      end

      subplot(2,2,1+i);
      cla
      hx=slice(dist_h.x, dist_h.y, dist_h.z, dist_h.D, 0.2, [], []);
      set(hx,'EdgeColor','none');
      p = patch(isosurface(dist_h.x,dist_h.y,dist_h.z,dist_h.D, 0));
      isonormals(dist_h.x,dist_h.y,dist_h.z,dist_h.D,p)
      set(p,'FaceColor','red','EdgeColor','none');
      daspect([1 1 1])
      view(3);
      axis([-1 1 -1 1 -1 1]);
      camlight
      lighting gouraud
      title({sprintf('Interpolated distance image (%s)', methods{i}), ...
             sprintf('Interpolation with %s took %.2fs', strrep(func2str(interp_fun), '_', '\_'), took(interpolator_i,size_i,i))});
      drawnow
    end
  end

end
%%
print -dpng ba_interp3.png
!trim-images ba_interp3.png -o cropped_; mv cropped_ba_interp3.png ba_interp3.png; convert ba_interp3.png -resize 500x500 ba_interp3_small.png
%%
figure(3);
clf
for i=1:numel(methods)
  subplot(numel(methods),1,i);
  errorbar(sizes.^3, mean(took(1,:,i,:),4), std(took(1,:,i,:),1,4), 'g.--', 'linewidth', 2, 'markersize', 20);
  hold on
  errorbar(sizes.^3, mean(took(2,:,i,:),4), std(took(2,:,i,:),1,4), 'b.--', 'linewidth', 2, 'markersize', 20);
  hold off
  time_steps = [0.1:0.1:0.9 1:1:9 10:5:50];
  axis([0 max(sizes.^3) 0 min(time_steps(time_steps>max(max(max(took(:,:,i,:))))))]);
  title({sprintf('Speed comparision: "%s" interpolation', methods{i}), sprintf('Average of %g runs', samples)});
  legend(cellfun(@(x) strrep(func2str(x), '_', '\_'), interpolators, 'uniformoutput', false), 'location', 'northwest');
  xlabel('Number of Interpolation points');
  ylabel('Elapsed time (s)');
  grid on
end
print -depsc ba_interp3_speed.eps
!convert ba_interp3_speed.eps ba_interp3_speed.png
!trim-images ba_interp3_speed.png -o cropped_; mv cropped_ba_interp3_speed.png ba_interp3_speed.png

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [ state ] = movePTPLineEefRelBase( t , Pos, relVel)
%% This function is used for moving the endeffector on a line, for the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] = movePTPLineEefRelBase( t , Pos, vel)

%% About:
% This function is used to perform linear relative movetion of the end-effector on a line,
% this performs only a linear motion, the orientation of end effector will not be changed.
% For more information refer to the example file (kuka0_move_lin_relative.m).

%% Arreguments:
% t: is the TCP/IP connection
% Pos: is the relative displacement, described in base frame, it is 1x3
% cell array, (unit millimeters)
% vel: is a double, represents the linear velocity of the motion (unit
% mm/sec).

% Copy right, Mohammad SAFEEA, 9th of May 2017

    theCommand=['jRelVel_',num2str(relVel),'_']; % set over ride.
    fprintf(t, theCommand);
    message=fgets(t);
    % The new position of end-effector, described in robot
    newPos={0,0,0,0,0,0};

    newPos{1}=Pos{1};
    newPos{2}=Pos{2};
    newPos{3}=Pos{3};

    sendEEfPositions( t ,newPos); % send destination joint positions.

    theCommand='doPTPinCSRelBase';
    fprintf(t, theCommand); % start the point to point motion.
    message=fgets(t);

    readingFlag=false;

    message='';
    state=false;
    while readingFlag==false
        message=fgets(t);
        state=checkAcknowledgment(message);
        if state
            readingFlag=true;
        end
    end

end


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

%% Syntax
% |[E_HEX20,V_HEX20,V_HEX20_cell,Fb_HEX20,Fb_HEX20_QUAD8]=hex8_hex20(E_HEX8,V_HEX8,V_HEX8_cell,Fb_HEX8);|

%% Description
% The |hex8_hex20| converts 4-node tetrahedral elements to 10-node
% tetrahedral elements.

%% Examples

%%
clear; close all; clc;

% Plot settings
fontSize=15;
faceAlpha=0.5;
edgeColor='k';
edgeWidth1=2;
edgeWidth2=1;
markerSize1=75;
markerSize2=10;

%% CONVERSION FROM HEX8 TO HEX20, EXAMPLE FOR A SINGLE HEXAHEDRON
% Creating a single 4-node HEXAHEDRON
[V8,~]=platonic_solid(2,1); %q indicates solid type, r is the radius
HEX8=[1:8];
[F8,~]=element2patch(HEX8,[],'hex8');

%%
% Converting to a single 10-node HEXAHEDRON
[HEX20,V20,~]=hex8_hex20(HEX8,V8,{});
[F20,~]=element2patch(HEX20,[],'hex20');

%%
% Plotting elements
cFigure;
subplot(1,2,1); hold on;
title('A linear HEXAHEDRON','FontSize',fontSize);

hp=gpatch(F8,V8,'gw','k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize1;

patchNormPlot(F8,V8,0.75); %Plotting face normals

for q=1:1:size(HEX8,2)
    text(V8(HEX8(1,q),1),V8(HEX8(1,q),2),V8(HEX8(1,q),3),['  ',num2str(q)],'FontSize',fontSize*2,'color','b');
end

axisGeom(gca,fontSize);
axis off;
camlight('headlight'); lighting flat;

subplot(1,2,2); hold on;
title('A quadratic HEXAHEDRON','FontSize',fontSize);

hp=gpatch(F20,V20,'rw','k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize1;
patchNormPlot(F20,V20,0.75); %Plotting face normals

for q=1:1:size(HEX20,2)
    text(V20(HEX20(1,q),1),V20(HEX20(1,q),2),V20(HEX20(1,q),3),[' ',num2str(q)],'FontSize',fontSize*2,'color','b');
end

axisGeom(gca,fontSize);
axis off;
camlight('headlight'); lighting flat;

drawnow;

%% CONVERSION FROM HEX8 TO HEX20, EXAMPLE FOR A HEXAHEDRON MESH WITH NODAL PARAMETERS
n=3;
for q=1:1:n
    [HEX8,V8]=subHex(HEX8,V8,1);
end

[F8,~]=element2patch(HEX8,[],'hex8');

%Create nodal result e.g. displacement and color
V4d=V8;
V4d(:,1)=V4d(:,1)+0.5.*sin(pi*V4d(:,3));
V4d(:,2)=V4d(:,2)+0.5.*cos(pi*V4d(:,3));
D4=V4d-V8;

C4=V8(:,1); %Color towards X

%%
% Converting to a single 10-node HEXAHEDRON
HEX8_cell={D4,C4};
[HEX20,V20,HEX20_cell]=hex8_hex20(HEX8,V8,HEX8_cell);
[F20,~]=element2patch(HEX20,[],'hex20');
D10=HEX20_cell{1};
C10=HEX20_cell{2};
V10d=V20+D10;

%%
% Plotting elements
hf=cFigure; % Open figure for plotting
subplot(2,2,1); hold on;
title('A linear HEXAHEDRON mesh','FontSize',fontSize);

hp=gpatch(F8,V8,'gw','k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize2;

axisGeom(gca,fontSize); axis off;
view([-50,12])
camlight('headlight'); lighting flat;

subplot(2,2,3); hold on;
title('Data on linear HEXAHEDRON mesh','FontSize',fontSize);

hp=gpatch(F8,V4d,C4,'k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize2;

axisGeom(gca,fontSize); axis off;
view([-50,12])
camlight('headlight'); lighting flat;

subplot(2,2,2); hold on;
title('A quadratic HEXAHEDRON mesh','FontSize',fontSize);

hp=gpatch(F20,V20,'rw','k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize2;

axisGeom(gca,fontSize); axis off;
view([-50,12])
camlight('headlight'); lighting flat;

subplot(2,2,4); hold on;
title('Mapped data on quadratic HEXAHEDRON mesh','FontSize',fontSize);

hp=gpatch(F20,V10d,C10,'k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize2;

axisGeom(gca,fontSize); axis off;
view([-50,12])
camlight('headlight'); lighting flat;

drawnow;

%% CONVERSION FROM HEX8 TO HEX20, EXAMPLE FOR KEEPING TRACK OF BOUNDARY FACES

[F8,~]=element2patch(HEX8,[],'hex8');
[indBoundary]=tesBoundary(F8);
Fb4=F8(indBoundary,:);

[HEX20,V20,HEX20_cell,Fb10]=hex8_hex20(HEX8,V8,{},Fb4);
[F20,~]=element2patch(HEX20,[],'hex20');

% [indBoundary]=tesBoundary(F20,V20);
% Fb20=F20(indBoundary,:);

%%
%

hf=cFigure; % Open figure for plotting
subplot(1,2,1); hold on;
title('A linear HEXAHEDRON mesh','FontSize',fontSize);

gpatch(F8,V8,'gw','k',faceAlpha);

gpatch(Fb4,V8,'rw','k',faceAlpha);
axisGeom(gca,fontSize); axis off;
camlight('headlight'); lighting flat;

subplot(1,2,2); hold on;
title('A quadratic HEXAHEDRON mesh','FontSize',fontSize);

hp=gpatch(Fb10,V20,'rw','k',faceAlpha);
hp.Marker='.';
hp.MarkerSize=markerSize2;

axisGeom(gca,fontSize); axis off;
camlight('headlight'); lighting flat;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function handles = dtiCmapInit(handles)
% Initialize color maps for mrDiffusion (dtiFiberUI)
%
% handles:  The dtiFiberUI handles
%
%   handles = dtiCmapInit(handles)
%
% Called in dtiFiberUI startup.  Not elsewhere.
%
% Brian (c) VISTASFOT Team, 2012

% I wonder if we need a cmap structure for the future?

handles.cmaps(1).name = 'autumn';
handles.cmaps(1).rgb = autumn(256);
handles.cmaps(2).name = 'winter';
handles.cmaps(2).rgb = winter(256);
handles.cmaps(3).name = 'spring';
handles.cmaps(3).rgb = spring(256);
handles.cmaps(4).name = 'summer';
handles.cmaps(4).rgb = summer(256);
handles.cmaps(5).name = 'hot';
handles.cmaps(5).rgb = hot(256);
handles.cmaps(6).name = 'cool';
handles.cmaps(6).rgb = cool(256);
handles.cmaps(7).name = 'gray';
handles.cmaps(7).rgb = gray(256);
handles.cmaps(8).name = 'red';
handles.cmaps(8).rgb = vividColormap(256,[0.15 0.85],'r');
handles.cmaps(9).name = 'green';
handles.cmaps(9).rgb = vividColormap(256,[0.15 0.85],'g');
handles.cmaps(10).name = 'blue';
handles.cmaps(10).rgb = vividColormap(256,[0.15 0.85],'b');
handles.cmaps(11).name = 'cyan';
handles.cmaps(11).rgb = vividColormap(256,[0.15 0.85],'c');
handles.cmaps(12).name = 'magenta';
handles.cmaps(12).rgb = vividColormap(256,[0.15 0.85],'m');
handles.cmaps(13).name = 'yellow';
handles.cmaps(13).rgb = vividColormap(256,[0.15 0.85],'y');

end


text:
function UFcollection_install (nlargefile)
%UFCOLLECTION_INSTALL install the UFcollection toolbox
%
% Example:
%   UFcollection_install
%
% See also UFget.

% Copyright 2006-2007, Timothy A. Davis

if (nargin < 1)
    % try with large-file I/O
    nlargefile = 0 ;
end

if (nlargefile)
    fprintf ('Trying to compile without large file support...\n') ;
    mex -I../UFconfig -DNLARGEFILE UFfull_write.c
else
    try
	mex -I../UFconfig UFfull_write.c
    catch
	fprintf ('Trying to compile without large file support...\n') ;
	mex -I../UFconfig -DNLARGEFILE UFfull_write.c
    end
end

addpath (pwd) ;
fprintf ('UFcollection toolbox successfully compiled.\n') ;


text:
function [matrix xmin xmax ymin ymax]=grd_read_v2(namefile)
% Function to read a GRD file
%                  (from Golden Software Surfer, ASCII format)
%
% [matrix xmin xmax ymin ymax]=grd_read_v2(name of file)
%
% Input:
%      nomarch = name of the file to be read, including ".grd" extension
% Output:
%      matrix =  matrix of the read data
%      xmin xmax ymin ymax = grid limits
%
% Coded by Alberto Avila Armella.
%          UPDATED & IMPROVED BY Jose Maria Garcia-Valdecasas

grdfile=fopen(namefile,'r');    % Open file
code=fgetl(grdfile);            % Reads surfer code 'DSAA'
% Grid dimensions (number of nodes)
aux=str2num(fgetl(grdfile)); nx=aux(1); ny=aux(2);
% X limits
aux=str2num(fgetl(grdfile)); xmin=aux(1); xmax=aux(2);
% Y limits
aux=str2num(fgetl(grdfile)); ymin=aux(1); ymax=aux(2);
% Z limits
aux=str2num(fgetl(grdfile)); zmin=aux(1); aux(2);
% Read matrix
[matrix,count] = fscanf(grdfile, '%f', [nx,ny]);
matrix=matrix';   % Trasposes matrix

fclose(grdfile);


text:
function [img1_2,img1_4,img1_8,img2_2,img2_4,img2_8,img2mask_2,img2mask_4,img2mask_8]=GPReduceAll(img1,img2,img2mask,levels,displayflag)
%
% [img1_2,img1_4,img1_8,img2_2,img2_4,img2_8]=GPReduceAll(img1,img2,levels)
%

global Gimg1 Gimg2 Gimg2mask Gimg1_2 Gimg1_4 Gimg1_8 Gimg2_2 Gimg2_4 Gimg2_8 Gimg2mask_2 Gimg2mask_4 Gimg2mask_8;

if ~isequal(img1,Gimg1) || isempty(Gimg1_2)
	disp('Down sampling image #1 for step 2...');
% 	if( size(img1,3) > size(img2,3) & size(img2,3)==16 )
% 		Gimg1_2 = GPReduce2(img1,16,displayflag);
% 	else
		Gimg1_2 = GPReduce(img1,displayflag);
% 	end
end

if ~isequal(img2,Gimg2) || isempty(Gimg2_2)
	disp('Down sampling image #2 for step 2...');
	Gimg2_2 = GPReduce(img2,displayflag);
end

if ~isequal(img2mask,Gimg2mask) || isempty(Gimg2mask_2)
	if isempty(img2mask)
		Gimg2mask_2 = [];
	else
		disp('Down sampling image #2 mask for step 2...');
		Gimg2mask_2 = GPReduce(img2mask,displayflag);
	end
end

if levels > 2
	if ~isequal(img1,Gimg1) || isempty(Gimg1_4)
		disp('Down sampling image #1 for step 3...');
% 		if( size(img1,3) > size(img2,3) & size(img2,3)==16 )
% 			Gimg1_4 = GPReduce2(Gimg1_2,8,displayflag);
% 		else
			Gimg1_4 = GPReduce(Gimg1_2,displayflag);
% 		end
	end

	if ~isequal(img2,Gimg2) || isempty(Gimg2_4)
		disp('Down sampling image #2 for step 3...');
		Gimg2_4 = GPReduce(Gimg2_2,displayflag);
	end

	if ~isequal(img2mask,Gimg2mask) || isempty(Gimg2mask_4)
		if isempty(img2mask)
			Gimg2mask_4 = [];
		else
			disp('Down sampling image #2 mask for step 3...');
			Gimg2mask_4 = GPReduce(Gimg2mask_2,displayflag);
		end
	end
end

if levels > 3
	if ~isequal(img1,Gimg1) || isempty(Gimg1_8)
		disp('Down sampling image #1 for step 4...');
% 		if( size(img1,3) > size(img2,3) & size(img2,3)==16 )
% 			Gimg1_8 = GPReduce2(Gimg1_4,4,displayflag);
% 		else
			Gimg1_8 = GPReduce(Gimg1_4,displayflag);
% 		end
	end

	if ~isequal(img2,Gimg2) || isempty(Gimg2_8)
		disp('Down sampling image #2 for step 4...');
		Gimg2_8 = GPReduce(Gimg2_4,displayflag);
	end

	if ~isequal(img2mask,Gimg2mask) || isempty(Gimg2mask_8)
		if isempty(img2mask)
			Gimg2mask_8 = [];
		else
			disp('Down sampling image #2 mask for step 4...');
			Gimg2mask_8 = GPReduce(Gimg2mask_4,displayflag);
		end
	end
end

Gimg1 = img1;
Gimg2 = img2;
Gimg2mask = img2mask;

img1_2 = Gimg1_2;
img1_4 = Gimg1_4;
img1_8 = Gimg1_8;
img2_2 = Gimg2_2;
img2_4 = Gimg2_4;
img2_8 = Gimg2_8;
img2mask_2 = Gimg2mask_2;
img2mask_4 = Gimg2mask_4;
img2mask_8 = Gimg2mask_8;

clear global Gimg1 Gimg2 Gimg2mask Gimg1_2 Gimg1_4 Gimg1_8 Gimg2_2 Gimg2_4 Gimg2_8 Gimg2mask_2 Gimg2mask_4 Gimg2mask_8;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function A = affineBuild(trans, rot, scale, skew)
%
% A = affineBuild(trans, rot, scale, skew)
%
% Builds an affine (orthogonal linear) transformation matrix from the four
% componets: translations, rotations scales and skews.
%
% trans, rot, scale, skew should all be 1x3 vectors specifying x,y,z.
% Rotations are in radians, and can be described as pitch (x-rotation),
% roll (y-rotation) and yaw (z-rotation).
%
% The transformations are applied in the following order:
%
% 1) translations
% 2) rotations
% 3) scaling
% 4) skews
%
% The form of the transform assumes a PRE-multiplication format:
% Y = A*X where X and Y are 4 x n arrays of n coordinates.
%
% HISTORY:
%   2004.03.12 RFD (bob@white.stanford.edu) shamelessly copied the core
%   algorithm from spm99's spm_matrix.m.
%   2005.07.08 ras (sayres at stanford edu) imported into mrVista 2.0
%   2005.07.28 ras enorces double class for inputs
if ~exist('skew','var') | isempty(skew), skew = [0 0 0];    end
if ~exist('scale','var') | isempty(scale), scale = [1 1 1]; end
if ~exist('rot','var') | isempty(rot), rot = [0 0 0];       end
if ~exist('trans','var') | isempty(trans), trans = [0 0 0]; end
if ~isa(trans,'double'), trans = double(trans);             end
if ~isa(rot,'double'), rot = double(rot);                   end
if ~isa(scale,'double'), scale = double(scale);             end
if ~isa(skew,'double'), skew = double(skew);                end

A  = eye(4);

A  = A*[1 	0 	0  trans(1);
        0 	1 	0  trans(2);
        0 	0 	1  trans(3);
        0 	0 	0  1];

A  = A*[1   0            0            0;
        0   cos(rot(1))  sin(rot(1))  0;
        0  -sin(rot(1))  cos(rot(1))  0;
        0   0    	      0           1];

A  = A*[cos(rot(2))  0  sin(rot(2))  0;
        0    	     1  0            0;
       -sin(rot(2))  0  cos(rot(2))  0;
        0            0  0            1];

A  = A*[cos(rot(3))  sin(rot(3))  0  0;
       -sin(rot(3))  cos(rot(3))  0  0;
        0            0            1  0;
        0     	     0    	      0  1];

A  = A*[scale(1) 0         0         0;
        0        scale(2)  0         0;
        0        0         scale(3)  0;
        0        0         0         1];

A  = A*[1  skew(1)  skew(2)  0;
        0  1        skew(3)  0;
        0  0        1        0;
        0  0        0        1];

return;


text:
function [Sim_Rate, Sim_Mse, Sim_Ber] = simulation_plot(Data, vec, Algorithms)
% function for plotting figures
% Data is a cell containing all the results
% vec is the x axis vector

[row,col] = size(Data);

for i = 1 : row
    for j = 1 : col
        data_one = Data{i,j};
        Sim_Rate(i,j) = data_one.Rate;
        Sim_Mse(i,j) = data_one.Mse;
        Sim_Ber(i,j) = data_one.Ber;
    end
end

plot_Rate;
plot_Mse;
plot_Ber;

    function plot_Rate
        figure;
        for m = 1: col
            plot(vec,Sim_Rate(:,m),'LineWidth',2);
            hold on
        end
        title('Rate');
        ylabel('Rate(bits/s/Hz');
        legend(Algorithms);
    end

    function plot_Mse
        figure;
        for m = 1: col
            semilogy(vec,Sim_Mse(:,m),'LineWidth',2);
            hold on
        end
        title('MSE');
        ylabel('MSE Value');
        legend(Algorithms);
    end

function plot_Ber
        figure;
        for m = 1: col
            semilogy(vec,Sim_Ber(:,m),'LineWidth',2);
            hold on
        end
        title('BER');
        ylabel('BER');
        legend(Algorithms);
    end
end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function r = quantity(a)

r = a.quantity;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function out = mrdivide(T,S)
% implements T / S
%

if isa(S,'double')

  T.M = T.M / S;
  out = T;

elseif isa(S,'tensor')

  T.M = T.M / S.M;
  out = T;

end

author: mtex-toolbox

repo: mtex

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


text:
function S = sparse2 (i,j,s,m,n,nzmax)					    %#ok
%SPARSE2 replacement for SPARSE
%
%   Example:
%   S = sparse2 (i,j,s,m,n,nzmax)
%
%   Identical to the MATLAB sparse function (just faster).
%   An additional feature is added that is not part of the MATLAB sparse
%   function, the Z matrix.  With an extra output,
%
%   [S Z] = sparse2 (i,j,s,m,n,nzmax)
%
%   the matrix Z is a binary real matrix whose nonzero pattern contains the
%   explicit zero entries that were dropped from S.  Z only contains entries
%   for the sparse2(i,j,s,...) usage.  [S Z]=sparse2(X) where X is full always
%   returns Z with nnz(Z) = 0, as does [S Z]=sparse2(m,n).  More precisely,
%   Z is the following matrix (where ... means the optional m, n, and nzmax
%   parameters).
%
%       S = sparse (i,j,s, ...)
%       Z = spones (sparse (i,j,1, ...)) - spones (S)
%
%   See also sparse.

%   Copyright 2006-2007, Timothy A. Davis, http://www.suitesparse.com

error ('sparse2 mexFunction not found') ;

author: ayonga

repo: frost-dev

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


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = 1 / m * X' * X;
[U, S, V] = svd(Sigma);

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

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
function sammi(model, parser, data, secondaries, options)
% Visualize the given model, set of reactions, and/or data using SAMMI.
% Documentation at: https://sammim.readthedocs.io/en/latest/index.html
%
% Citation: Schultz, A., & Akbani, R. (2019). SAMMI: A Semi-Automated
%     Tool for the Visualization of Metabolic Networks. Bioinformatics.
%
% USAGE:
% sammi(model,parser,data,secondaries,options)
%
% INPUT:
%   model: COBRA model to be visualized
%
% OPTIONAL INPUTS:
%   parser: How the model is to be parsed. There are four possible
%   options for this parameter. Default empty array.
%       *empty array: If this parameter is an empty arrray, all reaction
%       in the model will be loaded in a single map. Not advisable for
%       large maps.
%       *string: If this parameter is a characters array there are two
%       options. Either the parameter defines the path to a SAMMI map (JSON
%       file downloaded from a previous instance of SAMMI), in which case
%       the given map will be used, or the parameter defines a field in the
%       model struct, in which case this field will be used to parse the
%       model into subgraphs.
%       *cell array: If this parameter is a cell array, it should be a cell
%       array of strings containing reaction IDs. Only these reactions will
%       be included in a single SAMMI map.
%       *struct: If this model is a struct of length n, the model will be
%       parsed into n subgraphs. Each element of the struct should contain
%       two fields plus an additional optional one:
%           name: Name of the subgraph.
%           rxns: Reactions to be included in the subgraph.
%           flux: Optional field. Data to be mapped as reaction color.
%   data: Data to be mapped onto the model. Struct of length n. Defaults
%   to an empty array where no data is mapped. Element of the struct should
%   contain two fields:
%       type: A cell array of two strings. The firt string should be
%       either 'rxns', 'mets', or 'links' indicating which type of data
%       is to be mapped. The second string should be either 'color' or
%       'size', indicating how the data is to  be mapped. 'links' only
%       work with 'size', since link color is the same as the one of
%       the reaction it is assciated with.
%       data: a table object. VariableNames will be translated into
%       condition names, and RowNames should be reaction IDs for 'rxns'
%       and 'links' data, and metabolite IDs for 'mets' data. NaN values
%       will not be mapped.
%   secondaries: Cell array of strings of regular expressions. All
%   metabolites, in all subgraphs, matching any of the regular expressions
%   will be shelved. Default to empty array where no metabolites are
%   shelved.
%   options: Struct with the following fields:
%       htmlName: Name of the html file to be written and opened for the
%       visualization. Defaults to 'index_load'. Change this options to
%       write to a different html file that will not be overwritten by the
%       default option.
%       load: Load the html file in a new tab upon writing the file.
%       Default to true. If you would not like a new tab to open, set
%       this parameter to false and refresh a previously opened window. To
%       open a new window without re-running SAMMI use the openSammi
%       function.
%       jscode: String. Defaults to empty string. Additional JavaScript
%       code to run after loading the map. Can be any code to modify the
%       loaded map.
%
% OUTPUT:
%   No MATLAB output. Opens a browser window with the SAMMI visualization.
%
% EXAMPLES:
%   %1 Open model in single map
%   sammi(model)
%
%   %2 Open model as multiple subgraphs divided by subSystems
%   sammi(model,'subSystems')
%
%   %3 Open model as multiple subgraphs divided by subSystems, load two
%   %conditions with randomly generated data, and shelve hydrogen, water,
%   %and O2 upon loading.
%   rxntbl = array2table(randn(length(model.rxns),2),...
%       'VariableNames', {'condition1','condition2'},...
%       'RowNames', model.rxns);
%   data(1).type = {'rxns' 'color'};
%   data(1).data = rxntbl;
%   data(2).type = {'rxns' 'size'};
%   data(2).data = rxntbl;
%   secondaries = {'^h\[.\]$','^h20\[.\]$','^o2\[.\]$'};
%   sammi(model,'subSystems',data,secondaries)

if nargin < 2
    parser = [];
end
if nargin < 3
    data = [];
end
if nargin < 4
    secondaries = [];
end
if nargin < 5 || ~isfield(options,'htmlName')
    options.htmlName = 'index_load.html';
elseif isempty(regexp(options.htmlName,'\.html$'))
    options.htmlName = [options.htmlName '.html'];
end
if nargin < 5 || ~isfield(options,'load')
    options.load = true;
end
if nargin < 5 || ~isfield(options,'jscode')
    options.jscode = '';
end

%Read in index
sfolder = regexprep(which('sammi'),'sammi.m$','');
html = fileread([sfolder 'index.html']);

%Define options
if isstruct(parser)
    jsonstr = structParse(model,parser);
elseif ischar(parser) && exist(parser,'file') == 2 && ~isempty(regexp(parser,'\.json$','ONCE'))
    %Read map
    jsonstr = fileread(parser);
    jsonstr = strrep(jsonstr,'\','\\');
    %Add graph
    jsonstr = strcat('e = ',jsonstr,';\nreceivedTextSammi(JSON.stringify(e));');
elseif ischar(parser) && isfield(model,parser)
    ss = unique(model.(parser));
    if length(model.(parser)) == length(model.rxns)
        for i = 1:length(ss)
            dat(i).name = ss{i};
            dat(i).rxns = model.rxns(ismember(model.subSystems,ss{i}));
        end
    else
        for i = 1:length(ss)
            dat(i).name = ss{i};
            dat(i).rxns = model.rxns(sum(model.S(ismember(model.(parser),ss{i}),:)) ~= 0);
        end
    end
    jsonstr = structParse(model,dat);
elseif iscell(parser) || isempty(parser)
    if iscell(parser)
        %Keep only reactions we want
        model = removeRxns(model,model.rxns(~ismember(model.rxns,parser)));
    end
    %Convert model to sammi JSON string
    jsonstr = makeSAMMIJson(model);
    %Add graph
    jsonstr = strcat('e = ',jsonstr,';\nreceivedJSONwrapper(e)');
end

%Add data
for i = 1:length(data)
    if isequal(data(i).type{1},'rxns')
        if isequal(data(i).type{2},'color')
            datastring = makeSAMMIdataString(data(i).data);
            jsonstr = strcat(jsonstr,';\ndat = ',datastring,...
                ';\nreceivedTextFlux(dat)');
        elseif isequal(data(i).type{2},'size')
            datastring = makeSAMMIdataString(data(i).data);
            jsonstr = strcat(jsonstr,';\ndat = ',datastring,...
                ';\nreceivedTextSizeRxn(dat)');
        end
    end
    if isequal(data(i).type{1},'mets')
        if isequal(data(i).type{2},'color')
            datastring = makeSAMMIdataString(data(i).data);
            jsonstr = strcat(jsonstr,';\ndat = ',datastring,...
                ';\nreceivedTextConcentration(dat)');
        elseif isequal(data(i).type{2},'size')
            datastring = makeSAMMIdataString(data(i).data);
            jsonstr = strcat(jsonstr,';\ndat = ',datastring,...
                ';\nreceivedTextSizeMet(dat)');
        end
    end
    if isequal(data(i).type{1},'links')
        if isequal(data(i).type{2},'size')
            datastring = makeSAMMIdataString(data(i).data);
            jsonstr = strcat(jsonstr,';\ndat = ',datastring,...
                ';\nreceivedTextWidth(dat)');
        end
    end
end

%Shelve secondaries
if ~isempty(secondaries)
    secondaries = strrep(secondaries,'\','\\\\');
    jsonstr = strcat(jsonstr,';\nshelveList("(?:',strjoin(secondaries,')|(?:'),')");');
end

%Add last bit of code
jsonstr = strcat(jsonstr,';',options.jscode);

%Replace in html
html = strrep(html,'//MATLAB_CODE_HERE//',jsonstr);

%Account for speial characters
html = strrep(html,'%','%%');

%Write to file
fid = fopen([sfolder options.htmlName],'w');
fprintf(fid,html);
fclose(fid);

%Open window
if options.load
    web([sfolder options.htmlName],'-browser')
end
end

function jsonstr = structParse(model,parser)
    %Get only reactions we are using
    rx = {};
    for i = 1:length(parser); rx = unique(cat(1,rx,parser(i).rxns)); end
    model = removeRxns(model,model.rxns(~ismember(model.rxns,rx)));
    %Convert model to sammi JSON file
    jsonstr = makeSAMMIJson(model);
    %Add graph
    jsonstr = strcat('graph = ',jsonstr);
    %Make conversion vector
    convvec = makeSAMMIparseVector(parser);
    %Add parssing line
    jsonstr = strcat(jsonstr,';\ne = ',convvec,';\nfilterWrapper(e)');
end

author: opencobra

repo: cobratoolbox


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

    properties
        %Time
        t0(1,1) double = 0;
        tOffset(1,1) double = 0;

        %Constant
        const(1,1) double = 0;
        varConst(1,1) logical = false;
        constUb(1,1) double = 0;
        constLb(1,1) double = 0;

        %Neural Network
        net(1,1)
        varNetWB logical;
        netWbUb double = 0;
        netWbLb double = 0;
        hiddenSizes double = 5;
        gpuAvail(1,1) double = gpuDeviceCount("available");
    end

    methods
        function obj = FitNetModel(hiddenSizes)
            obj.hiddenSizes = hiddenSizes;
            obj.net = FitNetModel.createNetwork(hiddenSizes);

            netWb = abs(getwb(obj.net));
            obj.netWbUb = 100 * netWb;
            obj.netWbLb = 100 * (-netWb);
            obj.varNetWB = false(size(netWb));

            obj.gpuAvail = gpuDeviceCount("available");
        end

        function value = getValueAtTime(obj,ut)
            if(obj.gpuAvail > 0)
                gpuStr = 'yes';
            else
                gpuStr = 'no';
            end

            dt = (ut - obj.t0) + obj.tOffset;
            value = sim(obj.net,dt, 'UseGPU',gpuStr);
        end

        function wbParams = getWbParams(obj)
            wbParams = getwb(obj.net);
        end

        function value = getWbParamAtInd(obj, ind)
            wbParams = obj.getWbParams();
            value = wbParams(ind);
        end

        function setWbParamAtInd(obj, newVal, ind)
            wbParams = obj.getWbParams();
            wbParams(ind) = newVal;
            obj.net = setwb(obj.net, wbParams);
        end

        function listBoxStr = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.varNetWB))
                listBoxStr{i} = sprintf('Parameter %u', i); %#ok<AGROW>
            end
        end

        function newModel = deepCopy(obj)
            newModel = FitNetModel(obj.hiddenSizes);
            newModel.net = setwb(newModel.net, getwb(obj.net));

            newModel.t0 = obj.t0;
            newModel.tOffset = obj.tOffset;

            %Constant
            newModel.const = obj.const;
            newModel.varConst = obj.varConst;
            newModel.constUb = obj.constUb;
            newModel.constLb = obj.constLb;

            %Neural Network
            newModel.varNetWB = obj.varNetWB;
            newModel.netWbUb = obj.netWbUb;
            newModel.netWbLb = obj.netWbLb;
        end

        function t0 = getT0(obj)
            t0 = obj.t0;
        end

        function setT0(obj, newT0)
            obj.t0 = newT0;
        end

        function timeOffset = getTimeOffset(obj)
            timeOffset = obj.tOffset;
        end

        function setTimeOffset(obj, timeOffset)
            obj.tOffset = timeOffset;
        end

        function setConstValueForContinuity(obj, contValue)
            value = obj.getValueAtTime(obj.t0);

            obj.const = contValue - value;
        end

        function numVars = getNumVars(obj)
            numVars = 1 + numel(obj.varNetWB);
        end

        function x = getXsForVariable(obj)
            x = getwb(obj.net);
            x = x(obj.varNetWB);

            if(obj.varConst)
                x = [obj.const, x(:)'];
            end

            x = x(:)';
        end

        function updateObjWithVarValue(obj, x)
            if(obj.varConst)
                obj.const = x(1);
            end

            x = x(2:end);
            if(any(obj.varNetWB))
                wb = getwb(obj.net);
                wb(obj.varNetWB) = x(obj.varNetWB);
                obj.net = setwb(obj.net, wb);
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = [obj.constLb obj.netWbLb(:)'];
            ub = [obj.constUb obj.netWbUb(:)'];
        end

        function setBndsForVariable(obj, lb, ub)
            obj.constLb = lb(1);
            obj.constUb = ub(1);

            obj.netWbLb = lb(2:end);
            obj.netWbUb = ub(2:end);
        end

        function useTf = getUseTfForVariable(obj)
            useTf = [obj.varConst, obj.varNetWB(:)'];
        end

        function setUseTfForVariable(obj, useTf)
            obj.varConst = useTf(1);
            obj.varNetWB = useTf(2:end);
        end

        function nameStrs = getStrNamesOfVars(obj)
            x = getwb(obj.net);
            nameStrs = {'Constant Offset'};

            for(i=1:length(x)) %#ok<*NO4LP>
                nameStrs{end+1} = sprintf('Parameter %u', i); %#ok<AGROW>
            end
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = true;

            x = getwb(obj.net);
            for(i=1:length(x)) %#ok<*NO4LP>
                varsStoredInRad = horzcat(varsStoredInRad, false); %#ok<AGROW>
            end
        end

        function randomizeWB(obj)
            rands = obj.netWbLb + (obj.netWbUb - obj.netWbLb).*rand(size(obj.netWbLb));
            obj.net = setwb(obj.net, rands);
        end
    end

    methods(Static, Access=private)
        function net = createNetwork(hiddenSizes)
%             [x,t] = simplefit_dataset;

            net = fitnet(hiddenSizes);
            net = init(net);
%             net = configure(net,x,t);
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            obj.gpuAvail = gpuDeviceCount("available");
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function hand_translation ( )

%*****************************************************************************80
%
%% HAND_TRANSLATION applies a translation xy2 = xy + b to hand data.
%
%  Discussion:
%
%     This program assumes that the file 'HAND_NODES.TXT' is available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625.
%

%
%  Read the data.
%
  xy = load ( 'hand_nodes.txt' );
%
%  Make XY an array of column vectors.
%
  xy = xy';
%
%  Repeat the first column at the end so the polygon closes.
%
  xy = [ xy, [ xy(:,1) ] ];
%
%  Define the transformation as xy2 = xy + b.
%
  b = repmat ( [ 1.0; -0.5 ], 1, 60 );
%
%  Transform the data.
%
  xy2 = xy + b;
%
%  Clear the graphics frame.
%
  clf
%
%  Plot the original data.
%
  plot ( xy(1,:), xy(2,:), 'Color', 'r', 'LineWidth', 2 );
%
%  Plot the transformed data.
%
  hold on
  plot ( xy2(1,:), xy2(2,:), 'Color', 'b', 'LineWidth', 2 );
%
%  Annotate.
%
  axis equal
  grid on
  title ( 'Data (red), and data+[1.0;-0.5] (blue)' )

  hold off

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = seconds2human(secs, varargin)
%SECONDS2HUMAN( seconds )   Converts the given number of seconds into a
%                           human-readable string.
%
%   str = SECONDS2HUMAN(seconds) returns a human-readable string from a
%   given (usually large) amount of seconds. For example,
%
%       str = seconds2human(1463456.3)
%
%       str =
%       'About 2 weeks and 2 days.'
%
%   You may also call the function with a second input argument; either
%   'short' (the default) or 'full'. This determines the level of detail
%   returned in the string:
%
%       str = seconds2human(1463456.3, 'full')
%
%       str =
%       '2 weeks, 2 days, 22 hours, 30 minutes, 56 seconds.'
%
%   The 'short' format returns only the two largest units of time.
%
%   [secs] may be an NxM-matrix, in which case the output is an NxM cell
%   array of the corresponding strings.
%
%   NOTE: SECONDS2HUMAN() defines one month as an "average" month, which
%   means that the string 'month' indicates 30.471 days.
%
%   See also datestr, datenum, etime.

%   Author: Rody P.S. Oldenhuis
%   Delft University of Technology
%   E-mail: oldnhuis@dds.nl
%   Last edited 11/Feb/2010.

    % default error
    error(nargchk(1,2,nargin));%#ok

    % define some intuitive variables
    Seconds   = round(1                 );
    Minutes   = round(60     * Seconds  );
    Hours     = round(60     * Minutes  );
    Days      = round(24     * Hours    );
    Weeks     = round(7      * Days     );
    Months    = round(30.471 * Days     );
    Years     = round(365.26 * Days     );
    Centuries = round(100    * Years    );
    Millennia = round(10     * Centuries);

    % put these into an array, and define associated strings
    units   = [Millennia, Centuries, Years, Months, Weeks, ...
               Days, Hours, Minutes, Seconds];
    singles = {'millennium'; 'century'; 'year'; 'month'; ...
               'week'; 'day'; 'hour'; 'minute'; 'second'};
    plurals = {'millennia' ; 'centuries'; 'years'; 'months'; ...
               'weeks'; 'days'; 'hours'; 'minutes'; 'seconds'};

    % cut off all decimals from the given number of seconds
    assert(isnumeric(secs), 'seconds2human:seconds_mustbe_numeric', ...
        'The argument ''secs'' must be a scalar or matrix.');
    secs = round(secs);

    % parse second argument
    short = true;
    if (nargin > 1)
        % extract argument
        short = varargin{1};
        % check its type
        assert(ischar(short), 'seconds2human:argument_type_incorrect', ...
            'The second argument must be either ''short'' or ''full''.');
        % check its contents
        switch lower(short)
            case 'full' , short = false;
            case 'short', short = true;
            otherwise
                error('seconds2human:short_format_incorrect',...
                    'The second argument must be either ''short'' or ''full''.');
        end
    end

    % pre-allocate appropriate output-type
    numstrings = numel(secs);
    if (numstrings > 1), out = cell(size(secs)); end

    % build (all) output string(s)
    for j = 1:numstrings

        % initialize nested loop
        secsj   = secs(j);
        counter = 0;
        if short, string = 'About ';
        else      string = '';
        end

        % possibly quick exit
        if (secsj < 1), string = 'Less than one second.'; end

        % build string for j-th amount of seconds
        for i = 1:length(units)

            % amount of this unit
            amount = fix(secsj/units(i));

            % include this unit in the output string
            if amount > 0

                % increase counter
                counter = counter + 1;

                % append (single or plural) unit of time to string
                if (amount > 1)
                    string = [string, num2str(amount), ' ', plurals{i}];%#ok
                else
                    string = [string, num2str(amount), ' ', singles{i}];%#ok
                end

                % Finish the string after two units if short format is requested
                if (counter > 1 && short), string = [string, '.']; break, end%#ok

                % determine whether the ending should be a period (.) or a comma (,)
                if (rem(secsj, units(i)) > 0)
                    if short, ending = ' and ';
                    else ending = ', ';
                    end
                else ending = '.';
                end
                string = [string, ending];%#ok

            end

            % subtract this step from given amount of seconds
            secsj = secsj - amount*units(i);
        end

        % insert in output cell, or set output string
        if (numstrings > 1)
            out{j} = string;
        else
            out = string;
        end
    end % for

end % seconds2human

author: Sable

repo: mcbench-benchmarks


text:
function obj = matRad_exportDicomRTStruct(obj)
% matRad function to export dicom RT structure set.
% Class method of matRad_DicomExporter
%
% call
%   matRad_DicomExporter.matRad_exportDicomRTStruct()
%
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 matRad_cfg = MatRad_Config.instance();
matRad_cfg.dispInfo('Exporting DICOM RTStruct...\n');

env = matRad_getEnvironment();
isOctave = strcmp(env,'OCTAVE');

if isOctave
    matRad_cfg.dispWarning('RTStruct export currently not supported by matRad running in Octave due to crashing dicomwrite! Skipping...');
    return;
end

%% Metadata
%Class UID
ClassUID = '1.2.840.10008.5.1.4.1.1.481.3'; %RT Structure Set
meta.MediaStorageSOPClassUID = ClassUID;
meta.SOPClassUID             = ClassUID;
%TransferSyntaxUID = '1.2.840.10008.1.2.1'; %Explicit VR Little Endian - correct?
%meta.TransferSyntaxUID = TransferSyntaxUID;

%Identifiers
meta.SOPInstanceUID             = dicomuid;
meta.MediaStorageSOPInstanceUID = meta.SOPInstanceUID;
meta.SeriesInstanceUID          = dicomuid;
meta.SeriesNumber               = 1;
meta.InstanceNumber             = 1;

%Remaining Meta Data
meta.Modality = 'RTSTRUCT';
meta.Manufacturer = '';
meta.ReferringPhysicianName = obj.dicomName();
meta.OperatorsName = obj.OperatorsName;
meta.StationName = '';
meta = obj.assignDefaultMetaValue(meta,'ManufacturerModelName','matRad DicomExport');

meta.PatientName = obj.PatientName;
meta.PatientID = obj.PatientID;
meta.PatientBirthDate = obj.PatientBirthDate;
meta.PatientSex = obj.PatientSex;

%Name
meta = matRad_DicomExporter.assignDefaultMetaValue(meta,'StructureSetLabel','matRad_cst');
meta = matRad_DicomExporter.assignDefaultMetaValue(meta,'StructureSetName','matRad exported cst');

%ID of the Study
meta.StudyInstanceUID = obj.StudyInstanceUID;
meta.StudyID          = obj.StudyID;

%Dates & Times
currDate = now;
currDateStr = datestr(currDate,'yyyymmdd');
currTimeStr = datestr(currDate,'HHMMSS');
meta.InstanceCreationDate = currDateStr;
meta.InstanceCreationTime = currTimeStr;
meta.StudyDate = obj.StudyDate;
meta.StudyTime = obj.StudyTime;
meta = matRad_DicomExporter.assignDefaultMetaValue(meta,'StructureSetDate',currDateStr);
meta = matRad_DicomExporter.assignDefaultMetaValue(meta,'StructureSetTime',currTimeStr);

%Remaining stuff
meta.AccessionNumber = '';

%meta.PatientBirthDate = '';
%meta.PatientSex = 'O';
%meta.SoftwareVersion = '';

ct = obj.ct;

%Create X Y Z vectors if not present
if ~any(isfield(ct,{'x','y','z'}))
    %positionOffset = transpose(ct.cubeDim ./ 2);
    positionOffset = ct.cubeDim ./ 2;
    ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2);
    ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1);
    ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3);
end

%Since we are exporting HU directly --> no rescaling in any case
%meta.SliceThickness = ct.resolution.z;
%meta.PixelSpacing = [ct.resolution.y; ct.resolution.x];
%meta.ImageOrientationPatient = [1;0;0;0;1;0]; %lps

%meta.RescaleSlope = 1;
%meta.RescaleIntercept = 0;

obj.cst = matRad_computeVoiContoursWrapper(obj.cst,ct);

for i = 1:size(obj.cst,1)

    fprintf('Processinging VOI ''%s''...',obj.cst{i,2});
   %Select contours in axial slices
    contours        = obj.cst{i,7}(:,3);
    contourSliceIx  = find(~cellfun(@isempty,contours));
    contourSlicePos = ct.z(contourSliceIx);
    contours        = contours(contourSliceIx);

    % consider multiple contours of the same object
    for slice = 1:numel(contours)

       tmpContour = [];
       lower         = 1; % lower marks the beginning of a section
       cnt           = 0;

       while lower-1 ~= size(contours{slice,1},2)

          steps = contours{slice,1}(2,lower); % number of elements of current line section

          tmpContour(:,lower-cnt:lower-cnt+steps-1) = contours{slice,1}(:,lower+1:lower+steps);

          lower = lower+steps+1;
          cnt   = cnt + 1;
       end

       contours{slice,1} = tmpContour;

    end

    contours = cellfun(@(c) [ct.resolution.x; ct.resolution.y].* c + [ct.x(1) - ct.resolution.x;...
                             ct.y(1) - ct.resolution.y],contours,'UniformOutput',false);

    contours = cellfun(@(c,pos) addSlicePos(c,pos),contours,num2cell(contourSlicePos)','UniformOutput',false);
    %contours = cellfun(@transpose,contours,'UniformOutput',false);

    %Structure Definition
    ROISequenceItem.ROINumber = i;
    ROISequenceItem.ReferencedFrameOfReferenceUID = obj.FrameOfReferenceUID;
    ROISequenceItem.ROIName = obj.cst{i,2};
    ROISequenceItem.ROIGenerationAlgorithm = '';

    meta.StructureSetROISequence.(['Item_' num2str(i)]) = ROISequenceItem;

    %Contour Sequence
    if ~isOctave
        ROIContourSequenceItem.ROIDisplayColor = int32(round(255 * obj.cst{i,5}.visibleColor));
    end

    %Now create Contour subitems
    ContourSequence = struct;
    for c = 1:numel(contours)
        ContourItem = struct;

        %First store Image Sequence
        currCtSliceMeta = obj.ctSliceMetas(contourSliceIx(c));
        %currCtSliceMeta = matRad_DicomExporter.assignDefaultMetaValue(currCtSliceMeta,'SOPClassUID','1.2.840.10008.5.1.4.1.1.2');

        %Not sure about this
        ContourItem.ContourImageSequence.Item_1.ReferencedSOPClassUID = currCtSliceMeta.SOPClassUID; %We are referencing the CT
        ContourItem.ContourImageSequence.Item_1.ReferencedSOPInstanceUID = currCtSliceMeta.SOPInstanceUID; %TODO: ID of the respective ctSlice

        %Now remaing data
        ContourItem.ContourGeometricType = 'CLOSED_PLANAR';
        ContourItem.NumberOfContourPoints = size(contours{c},2);
        ContourItem.ContourData = contours{c}(:);

        %Now store Item
        ContourSequence.(['Item_' num2str(c)]) = ContourItem;
    end
    ROIContourSequenceItem.ContourSequence = ContourSequence;

    ROIContourSequenceItem.ReferencedROINumber = i;

    %Store to meta data
    meta.ROIContourSequence.(['Item_' num2str(i)]) = ROIContourSequenceItem;

    %RTROI Observation Sequence
    RTROIObservationsSequenceItem.ObservationNumber = i;
    RTROIObservationsSequenceItem.ReferencedROINumber = i;
    RTROIObservationsSequenceItem.ROIObservationLabel = obj.cst{i,2};

    if strcmp(obj.cst{i,3},'TARGET')
        if ~isempty(regexpi(obj.cst{i,2},['(' strjoin(obj.targetPtvDict) ')']))
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'PTV';
            fprintf('identified target type as PTV...');
        elseif ~isempty(regexpi(obj.cst{i,2},['(' strjoin(obj.targetGtvDict) ')']))
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'GTV';
            fprintf('identified target type as GTV...');
        elseif ~isempty(regexpi(obj.cst{i,2},['(' strjoin(obj.targetGtvDict) ')']))
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'CTV';
            fprintf('identified target type as CTV...');
        else
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'CTV';
            fprintf('Defaulting target type to CTV...');
        end
    else
        if ~isempty(regexpi(obj.cst{i,2},['(' strjoin(obj.externalContourDict) ')']))
            fprintf('automatically identified as External Contour...');
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'EXTERNAL';
        else
            RTROIObservationsSequenceItem.RTROIInterpretedType = 'AVOIDANCE';
        end
    end

    RTROIObservationsSequenceItem.ROIInterpreter = obj.dicomName();

    meta.RTROIObservationsSequence.(['Item_' num2str(i)]) = RTROIObservationsSequenceItem;

    fprintf('Done!\n');
    %matRad_progress(i,size(obj.cst,1));
end

for i = 1:numel(obj.ctSliceMetas)
    ImageSequenceItem.ReferencedSOPClassUID = obj.ctSliceMetas(i).SOPClassUID;
    ImageSequenceItem.ReferencedSOPInstanceUID = obj.ctSliceMetas(i).SOPInstanceUID;
    ContourImageSequence.(['Item_' num2str(i)]) = ImageSequenceItem;
end

RTReferencedSeriesSequenceItem.SeriesInstanceUID = obj.ctSliceMetas(1).SeriesInstanceUID;
RTReferencedSeriesSequenceItem.ContourImageSequence = ContourImageSequence;

RTReferencedStudySequenceItem.ReferencedSOPClassUID = '1.2.840.10008.3.1.2.3.2'; %Apparently this class UID is deprecated in DICOM standard - what to use instead?
RTReferencedStudySequenceItem.ReferencedSOPInstanceUID = obj.StudyInstanceUID;
RTReferencedStudySequenceItem.RTReferencedSeriesSequence.Item_1 = RTReferencedSeriesSequenceItem;

meta.ReferencedFrameOfReferenceSequence.Item_1.FrameOfReferenceUID = obj.FrameOfReferenceUID;
meta.ReferencedFrameOfReferenceSequence.Item_1.RTReferencedStudySequence.Item_1 = RTReferencedStudySequenceItem;

filename = 'RTstruct.dcm';
filepath = obj.dicomDir;
filename = fullfile(filepath,filename);

if isOctave
	dicomwrite(int16(zeros(2)),filename,meta);
else
    obj.rtssExportStatus = dicomwrite([],filename,meta,'CreateMode','copy');%,'TransferSyntax',TransferSyntaxUID);
end

obj.rtssMeta = meta;
end

function c = addSlicePos(c,slicePos)
    nPoints  = size(c,2);
    slicePos = slicePos * ones(1,nPoints);
    c(3,:)   = slicePos;
end

repo: matRad