text:
 function h = im(varargin)
%function h = im([options,] [xx,] [yy,] zz, [scale|clim,] [title])
%| show matrix zz as an image, possibly with (x,y) axes labeled by xx,yy
%|
%| options:
%|	subplot		an integer for subplot
%|	'notick'	no axis tick marks
%|	'black0'	make sure image value of 0 is black (max white)
%|	'blue0'		show zeros as blue
%|	'mid3'		show middle 3 planes of 3D volume
%|	'mip3'		show 3 MIP (max. intens. proj.) views of 3D volume
%|	'row', n	# of rows for 3D montages, for this call
%|	'col', n	# of cols for 3D montages, for this call
%|
%| after options:
%|	scale		scale image by this factor
%|	clim		limits for colorbar
%|	title		axis title, processed by titlef()
%|	'cbar'		add colorbar
%|
%| one argument commands:
%|	'off'		disable display,
%|	'off-quiet'	disable display, and do not print 'disabled' messages
%|	'on'		enable display, 'ison' report if enabled
%|	'clf'		clear display (if enabled)
%|	'colorneg'	show negatives as red
%|	'db40'		log grayscale over 40dB range
%|	'nan-warn'	warning if nan value(s), 'nan-fail' fail if nan
%|	'drawnow'	'drawnot' toggle immediate drawing
%|	'tickon'	'tickoff' toggle default tick drawing
%|	'tick1'		axis ticks [1 N] (default)
%|	'tick0'		axis ticks [0 N-1]
%|	'state'		report state
%|
%| other commands:
%|	'pl' sub_m sub_n specify subplot layout for subsequent plots
%|	'plc'		like 'pl' but first clf
%|	'pl-tight'	like 'pl' but tight plots with little (if any) spacing
%|	'subplot' plotindex	choose a subplot
%|	n2min	n		<= this min # for dim2 we plot instead (def: 1)
%|	zero blue	replace 1st colormap entry (zero?) with blue
%|	row n		# of rows for 3D montages - setting state
%|	col n		# of cols for 3D montages - setting state
%|
%| Copyright 1997, Jeff Fessler, University of Michigan

% handle states
persistent state
if ~isvar('state') || isempty(state)
	state = ir_im_reset;
end

% default is to give help
if ~nargin && ~nargout
	help(mfilename)
	if im, printm('im enabled'), else, printm('im disabled'), end
	fail(mfilename)
end

% for conditional of the form 'if im, ..., end'
if ~nargin && nargout
	h = state.display;
return
end

% im row 2 or im col 3
if nargin == 2 && (streq(varargin{1}, 'row') || streq(varargin{1}, 'col'))
	tmp = varargin{2};
	if ischar(tmp), tmp = str2num(tmp); end
	state.montage = {varargin{1}, tmp};
return
end

% process single string command arguments
if nargin == 1 && ischar(varargin{1})
	arg = varargin{1};

	switch arg
	case 'on'
		state.display = true;
		state.display_quiet = false;
		disp 'enabling images'
	case 'off'
		state.display = false;
		state.display_quiet = false;
		disp 'disabling images'
	case 'off-quiet'
		state.display = false;
		state.display_quiet = true;
	case 'clf'
		if state.display
			state.sub_m = [];
			clf
		end
	case 'drawnow'
		state.drawnow = true;
	case 'drawnot'
		state.drawnow = false;
	case 'subplot'
		if state.sub_m
			im_subplot(state, state.next_sub);
			state.next_sub = state.next_sub + 1;
		else
			fail('subplot not set')
		end
	case 'tickon'
		state.tick = true;
	case 'tickoff'
		state.tick = false;
	case 'tick1'
		state.ticks = 'tick1'; % [1 N]
	case 'tick0'
		state.ticks = 'tick0'; % [0 N-1]
	case 'tickc'
		state.ticks = 'tickc'; % [-(N-1)/2 0 (N-1)/2]
	case 'ison' % query
		if state.display
			h = true;
		else
			h = false;
		end
	case 'state' % query
		h = state;
	case 'nan-fail'
		state.nan_fail = true;
	case 'nan-warn'
		state.nan_fail = false;
	case 'blue0'
		state.blue0 = true;
	case 'colormap_keep'
		state.colormap_control = '';
	case 'colorneg'
		state.colorneg = true;
	case 'reset'
		state = ir_im_reset;
	case 'test'
		im_test, return
	otherwise
		if streq(arg, 'db', 2)
			state.db = sscanf(arg, 'db%g');
		else
			fail('unknown argument: "%s"', arg)
		end
	end
return
end

% 'pl' 'plc' 'pl-tight' option
if streq(varargin{1}, 'pl') || streq(varargin{1}, 'pl-tight') ...
	|| streq(varargin{1}, 'plc')
	if streq(varargin{1}, 'plc'), im clf, end
	if nargin == 3
		state.sub_m = ir_ensure_num(varargin{2});
		state.sub_n = ir_ensure_num(varargin{3});
	elseif nargin == 1
		state.sub_m = [];
		state.sub_n = [];
	else
		fail 'bad pl usage'
	end
	state.pl_tight = streq(varargin{1}, 'pl-tight');
	state.next_sub = 1;
return
end

% 'subplot' option
if streq(varargin{1}, 'subplot')
	state = im_subplot(state, ir_ensure_num(varargin{2}));
return
end

% 'n2min' option
if streq(varargin{1}, 'n2min')
	state.n2min = ir_ensure_num(varargin{2});
return
end

% 'zero' or 'zero blue'
if streq(varargin{1}, 'zero')
	cmap = get(gcf, 'colormap');
	cmap(1,:) = [0 0 1];
	set(gcf, 'colormap', cmap);
return
end

scale = 1;
titlearg = {};
opt.cbar = false;
clim = [];
xx = [];
yy = [];
opt.blue0 = false; % 1 to color 0 as blue
opt.colorneg = false; % 1 to put negatives in color and 0=blue
opt.db = 0; % nonzero to use a log scale
isxy = false; % 1 if user provides x and y coordinates
isplot = false;
opt.tick = state.tick;
opt.ticks = state.ticks;
opt.black0 = false;
opt.mid3 = false; % 1 if show mid-plane slices of 3D object
opt.mip3 = false; % 1 if show 3 MIP views of 3D object
opt.montage = state.montage;
opt.hsv = false;

if length(varargin) == 1 && isempty(varargin{1})
	fail 'empty argument?'
end

% optional arguments
zz_arg_index = 1;
while length(varargin)
	arg = varargin{1};
	if isempty(arg)
		0; % do nothing
	elseif max(size(arg)) == 1 % scalar
		if state.display
			arg1 = varargin{1};
			if arg1 > 99 % reset
				state.sub_m = [];
				state.sub_n = [];
			end
			if isempty(state.sub_m)
				if arg1 >= 111
					subplot(arg1)
				else
					printm('ignoring subplot %d', arg1)
				end
			else
				state = im_subplot(state, arg1);
			end
		end
	elseif streq(arg, 'notick')
		opt.tick = false;
	elseif streq(arg, 'tick1')
		opt.ticks = 'tick1'; % [1 N]
	elseif streq(arg, 'tick0')
		opt.ticks = 'tick0'; % [0 N-1]
	elseif streq(arg, 'tickc')
		opt.ticks = 'tickc'; % [-(N-1)/2 0 (N-1)/2]
	elseif streq(arg, 'colorneg')
		opt.colorneg = true;
	elseif streq(arg, 'db', 2)
		opt.db = sscanf(arg, 'db%g');
	elseif streq(arg, 'black0')
		opt.black0 = true;
	elseif streq(arg, 'blue0')
		opt.blue0 = true;
	elseif streq(arg, 'hsv')
		opt.hsv = true;
	elseif streq(arg, 'mid3')
		opt.mid3 = true;
	elseif streq(arg, 'mip3')
		opt.mip3 = true;
	elseif streq(arg, 'row') || streq(arg, 'col')
		opt.montage = {arg, varargin{2}};
		varargin = {varargin{2:end}};
		zz_arg_index = 1 + zz_arg_index;
	else
		break
	end

	varargin = {varargin{2:end}};
	zz_arg_index = 1 + zz_arg_index;
end

if length(varargin) < 1, ir_usage, end

if ~isempty(state.sub_m) % possibly redundant but ok
	nplot = state.sub_m * state.sub_n;
	im_subplot(state, 1 + mod(state.next_sub-1, nplot));
	state.next_sub = state.next_sub + 1;
end

% plotting vector(s)
if ndims(varargin{1}) <= 2 && min(size(varargin{1})) <= state.n2min ...
	&& length(varargin) < 3
		isplot = true;
		plot_data = varargin{1};

% xx, yy
elseif ndims(varargin{1}) <= 2 && min(size(varargin{1})) == 1
	xx = col(varargin{1});
	if length(varargin) < 2, help(mfilename), fail 'need both xx,yy', end
	if min(size(varargin{2})) ~= 1, fail 'both xx,yy need to be 1D', end
	yy = col(varargin{2});
	varargin = {varargin{3:end}};
	isxy = 1;
	zz_arg_index = 2 + zz_arg_index;
end

if ~state.display
	if ~state.display_quiet
		printm(['disabled: ' inputname(zz_arg_index)])
	end
return
end

% zz
if length(varargin)
	zz = varargin{1};
	if iscell(zz); zz = stackup(zz{:}); isplot = 0; end % trick
	zz = double(zz);
	if ndims(zz) > 2 && min(size(zz)) == 1
		zz = squeeze(zz); % handle [1 n2 n3] case as [n2 n3]
	end
	varargin = {varargin{2:end}};
else
	fail 'no image?'
end

if any(isnan(zz(:)))
	tmp = sprintf('image contains %d NaN values of %d!?', ...
		sum(isnan(zz(:))), numel(zz));
	if state.nan_fail, fail(tmp), else, warning(tmp), end
end

if any(isinf(zz(:)))
	tmp = sprintf('image contains %d Inf values of %d!?', ...
		sum(isinf(zz(:))), numel(zz));
	if state.nan_fail, fail(tmp), else, warning(tmp), end
end

if opt.mid3
	mip3_size = size(zz);
	pn = jf_protected_names;
	zz = pn.mid3(zz);
end

if opt.mip3
	mip3_size = size(zz);
	zz = jf_mip3(zz);
end

% title, scale, clim, cbar
while length(varargin)
	arg = varargin{1};
	if isempty(arg)
		% do nothing
	elseif ischar(arg)
		if streq(arg, 'cbar')
			opt.cbar = true;
		else
			titlearg = {arg};
		end
	elseif isnumeric(arg) % isa(arg, 'double')
		if max(size(arg)) == 1
			scale = arg;
		elseif all(size(arg) == [1 2])
			clim = arg;
		else
			fail 'nonscalar scale / nonpair clim?'
		end
%		pr scale
	else
		fail 'unknown arg'
	end
	varargin = {varargin{2:end}};
end

if isplot
	plot(plot_data, state.line1type)
	if ~isempty(titlearg), titlef(titlearg{:}), end
return
end

if issparse(zz), zz = full(zz); end
if ~isreal(zz)
	zz = abs(zz);
	printf('warn %s: magnitude of complex image', mfilename)
end

if opt.db || state.db
	if ~opt.db && state.db, opt.db = state.db, end
	if max(zz(:)) <= 0, fail 'db for negatives?', end
	zz = abs(zz); zz = max(zz, eps);
	zz = 10*log10(abs(zz) / max(abs(zz(:))));
	zz(zz < -opt.db) = -opt.db;
end

zmax = max(zz(:));
zmin = min(zz(:));

if opt.black0
	if ~isempty(clim)
		warning 'black0 overrules clim'
	end
	clim = [0 zmax];
end

if scale ~= 1
	if scale == 0
		zmin = 0;
		scale = 1;
	elseif scale < 0
		zmin = 0;
		scale = -scale;
	end
end

%% determine colormap and adjust data zz if necessary
function [zz, cmap] = im_cmap(zz, opt, state)
	cmap = [];
	if opt.blue0 || state.blue0
		cmap = [[0 0 1]; gray(256)];
%		colormap_gca(cmap)
%		colormap_cmap = cmap;
		zt = zz;
		zz(zt > 0) = 1+floor(255 * zz(zt > 0) / (abs(max(zt(:))) + eps));
		zz(zt == 0) = 1;

	elseif (opt.colorneg || state.colorneg) && ~ir_is_octave
		if 1 % original
			cmap = hot(512);	cmap = flipud(cmap(1:256,:));
			cmap = [cmap; [0 0 1]; gray(256)];
		else % +green -red
			tmp = [0:255]'/255;
			cmap = [flipud(tmp)*[1 0 0]; [0 0 1]; tmp*[0 1 0]];
		end
	%	colormap_gca(cmap)
	%	colormap_cmap = cmap;
		zt = zz;
		zz(zt > 0) = 257+1+floor(255 * zz(zt > 0) / (abs(max(zt(:))) + eps));
		zz(zt == 0) = 257;
		zz(zt < 0) = 1+floor(-255 * zz(zt < 0) / (abs(min(zt(:))) + eps));

	elseif opt.hsv
	%	colormap_gca(hsv(256))
		cmap = hsv(256);

	else
		switch state.colormap_control
		case 'gray'
		%	colormap_gca(gray(256))
			cmap = gray(256);
		case ''
			% do nothing
		otherwise
			fail('unknown colormap_control "%s"', state.colormap_control)
		end
	end
end

[zz, colormap_cmap] = im_cmap(zz, opt, state);

if scale ~= 1 % fix: use clim?
	n = size(colormap,1);
	if zmax ~= zmin
		zz = (n - 1) * (zz - zmin) / (zmax - zmin); % [0,n-1]
	else
		if zmin == 0
			zz(:) = 0;
			clim = [0 1];
		else
			zz(:) = n - 1;
		end
	end
	zz = 1 + round(scale * zz);
	zz = min(zz,n);
	zz = max(zz,1);
elseif zmin == zmax
	if zmin == 0
		clim = [0 1];
	else
		zz(:) = 1;
		clim = [0 1];
	end
end

if ndims(zz) < 3 % 2d
	zz = zz'; % trick: transpose for consistency with C programs
	if ir_is_octave && state.octave_yflip
		zz = flipdim(zz, 1); % trick: note zz' above
	end
	if isxy
		if length(xx) ~= size(zz,2), fail 'xx size', end
		if length(yy) ~= size(zz,1), fail 'yy size', end
		hh = image(xx, yy, zz, 'CDataMapping', 'scaled');
		if opt.tick
			switch opt.ticks
			case 'tick0'
				setgca(	'xtick', xx([1 end]), ...
					'ytick', yy([1 end]), ...
					'fontsize', ir_fontsize('im_axes'));
			case 'tickc'
				setgca('xtick', xx([1 ceil((end+1)/2) end]), ...
				'ytick', yy([1 ceil((end+1)/2) end]), ...
				'fontsize', ir_fontsize('im_axes'));
			end
		end
	else
		n1 = size(zz,2);
		n2 = size(zz,1);
		switch opt.ticks
		case 'tick1'
			xx = 1:n1;
			yy = 1:n2;
			xtick = [1 n1];
			ytick = [1 n2];
		case 'tick0'
			xx = 0:(n1-1);
			yy = 0:(n2-1);
			xtick = [0 n1-1];
			ytick = [0 n2-1];
		case 'tickc' % [- 0 +]
			xtick = [-1 0 1] * (n1-1)/2;
			ytick = [-1 0 1] * (n2-1)/2;
			xx = (-(n1-1)/2):((n1-1)/2);
			yy = (-(n2-1)/2):((n2-1)/2);
		end
		hh = image(xx, yy, zz, 'CDataMapping', 'scaled');
%		hh = image(zz, 'CDataMapping', 'scaled');
%		setgca('CLimMode','auto')

		% unclutter axes by only showing end limits
%		setgca('xtick', [1 n1], 'ytick', [1 n2])

		if opt.tick
			setgca('xtick', [], 'ytick', ytick)
			if is_pre_v7
				text(1, 1.04*(n2+0.5), '1', ...
					'horizontalalign', 'left', ...
					'verticalalign', 'top')
				text(n1, 1.04*(n2+0.5), num2str(n1), ...
					'horizontalalign', 'right', ...
					'verticalalign', 'top')
			elseif opt.mip3 || opt.mid3
				tmp = [1 mip3_size(1)+[0 mip3_size(3)]];
				setgca('xtick', tmp)
				tmp = [1 mip3_size(2)+[0 mip3_size(3)]];
				setgca('ytick', tmp)
			else
				setgca('xtick', xtick)
			end
			tmp = {'fontsize', ir_fontsize('im_axes')};
			setgca(tmp{:})
		end
		% problem with this is it fails to register
		% space used for 'ylabel'
		% so i stick with manual xtick since that is what overlaps
		% with the colorbar
		if 0 && opt.tick
			text(-0.01*n1, n2, '1', ...
				'verticalalign', 'bottom', ...
				'horizontalalign', 'right')
			text(-0.01*n1, 1, num2str(n2), ...
				'verticalalign', 'top', ...
				'horizontalalign', 'right')
		end
	end

else % 3d
	n1 = size(zz,1);
	n2 = size(zz,2);
	if 0 % white border
		zz(:,end+1,:) = max(zz(:));
		zz(end+1,:,:) = max(zz(:));
	end
	if ir_is_octave && state.octave_yflip
		zz = flipdim(zz, 2);
	end
	dimz = size(zz);
	zz = montager(zz, opt.montage{:});

	if isxy
		xx3 = [0:size(zz,1)/dimz(1)-1]*(xx(2)-xx(1))*dimz(1);
		xx3 = col(outer_sum(xx, xx3));
		yy3 = [0:size(zz,2)/dimz(2)-1]*(yy(2)-yy(1))*dimz(2);
		yy3 = col(outer_sum(yy, yy3));
		hh = image(xx3, yy3, zz', 'CDataMapping', 'scaled');
		if opt.tick && n1 > 1 && n2 > 1 % unclutter
			setgca('xtick', [xx(1) xx(end)], 'ytick', ...
				sort([yy(1) yy(end)]))
		end
		axis xy
	else
		hh = image(zz', 'CDataMapping', 'scaled');
		if opt.tick && n1 > 1 && n2 > 1 % unclutter
			setgca('xtick', [1 n1], 'ytick', [1 n2])
		end

		if state.line3plot % lines around each sub image
			m1 = size(zz,1) / n1;
			m2 = size(zz,2) / n2;
			plot_box = @(ox,oy,arg) plot(ox+[0 1 1 0 0]*n1+0.5, ...
					oy+[0 0 1 1 0]*n2+0.5, state.line3type, arg{:});
			n3 = dimz(3);
			for ii=0:n3-1
				i1 = mod(ii, m1);
				i2 = floor(ii / m1);
				hold on
				plot_box(i1*n1, i2*n2, ...
					{'linewidth', state.line3width})
				hold off
			end
		end
	end
end

if ~isempty(colormap_cmap)
	colormap_gca(colormap_cmap)
end

if (zmax == zmin)
%	fprintf('Uniform image %g [', zmin)
%	fprintf(' %g', size(zz))
%	fprintf(' ]\n')
	texts(0.5, 0.5, sprintf('Uniform: %g', zmin), 'center', 'color', 'blue')
end

if (opt.colorneg || state.colorneg) && ~ir_is_octave
	set(hh, 'CDataMapping', 'direct')
else
	if ~isempty(clim),
		setgca('CLim', clim)
	elseif ~ishold
		setgca('CLimMode', 'auto')
	end
end

setgca('tickdir', 'out')
%setgca('tickdir', 'in')

if nargout > 0
	h = hh;
end

% axis type depends on whether user provides x,y coordinates
if isxy
	if length(xx) > 1 && length(yy) > 1 && ...
		abs(xx(2)-xx(1)) == abs(yy(2)-yy(1)) % square pixels
		axis image
	end
	axis xy
else
	axis image
	if ir_is_octave && ~isvar('dimz') % 2d
		axis([0.5 n1+0.5 0.5 n2+0.5])
	end
end

if ~isempty(titlearg)
	titlef(titlearg{:})
else % default title
	tmp = inputname(zz_arg_index);
	tmp = sprintf('%s range: [%.3g %.3g]', tmp, zmin, zmax);
	title(tmp, 'interpreter', 'none')
end

if ~opt.tick
	setgca('xtick', [], 'ytick', [])
end

if opt.cbar, cbar, end

if state.drawnow, drawnow, end

end % im()

% ir_im_reset()
% reset state to its defaults
function state = ir_im_reset
	state.display = true;
	state.display_quiet = false;
	state.colorneg = false;
	state.db = 0;
	state.blue0 = false;
	state.nan_fail = false;
	state.drawnow = false;
	state.tick = true;
	state.ticks = 'tick1'; % [1 N] ticks, or [0 N-1] if 'tick0'
	state.octave_yflip = true; % 3.8.2 octave seems to ignore 'axis ij' and ydir
	state.montage = {}; % for montager
	state.next_sub = 1; % next subplot index % todo
	state.sub_m = []; % for subplot
	state.sub_n = [];
	state.n2min = 1;
	state.pl_tight = false;
	state.line3plot = true;
	state.line3type = 'y-';
	state.line3width = 1;
	state.line1type = '-'; % for 1D plots
	state.colormap_control = 'gray'; % default is to set it to gray(256) each time
end

function x = ir_ensure_num(x)
	if ischar(x), x = str2num(x); end
end

function colormap_gca(cmap)
	if ir_is_octave || isfreemat
		colormap(cmap)
	else
		colormap(gca, cmap)
	end
end

function state = im_subplot(state, num)
	if state.display
		if state.pl_tight
			num = num - 1;
			x = 1 / state.sub_n;
			y = 1 / state.sub_m;
			ny = floor(num / state.sub_n);
			nx = num - ny * state.sub_n;
			ny = state.sub_m - ny - 1;
			subplot('position', [nx*x ny*y x y])
		else
			if num > state.sub_m * state.sub_n
				warn('num=%d > nsubplot=%d; resetting', ...
					num, state.sub_m * state.sub_n)
				num = 1;
			end
			subplot(state.sub_m, state.sub_n, num)
		end
		state.next_sub = num;
	end
end

function setgca(varargin)
	set(gca, varargin{:})
end

function im_test
	im clf, im(rand(6))
	im clf, im pl-tight 2 3
	im(1, rand(3))
	im(2, rand(3))
	im(5, ones(3))
	prompt
	im clf, im(rand(2^7,2^7-2,4))
end


text:
function drawLine(p1, p2, varargin)
%DRAWLINE Draws a line from point p1 to point p2
%   DRAWLINE(p1, p2) Draws a line from point p1 to point p2 and holds the
%   current figure

plot([p1(1) p2(1)], [p1(2) p2(2)], varargin{:});

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [cout,hout] = tricontour(p,t,Hn,N)
% Contouring for functions defined on triangular meshes
%
%   TRICONTOUR(p,t,F,N)
%
% Draws contours of the surface F, where F is defined on the triangulation
% [p,t]. These inputs define the xy co-ordinates of the points and their
% connectivity:
%
%   P = [x1,y1; x2,y2; etc],            - xy co-ordinates of nodes in the
%                                         triangulation
%   T = [n11,n12,n13; n21,n23,n23; etc] - node numbers in each triangle
%
% The last input N defines the contouring levels. There are several
% options:
%
%   N scalar - N number of equally spaced contours will be drawn
%   N vector - Draws contours at the levels specified in N
%
% A special call with a two element N where both elements are equal draws a
% single contour at that level.
%
%   [C,H] = TRICONTOUR(...)
%
% This syntax can be used to pass the contour matrix C and the contour
% handels H to clabel by adding clabel(c,h) or clabel(c) after the call to
% TRICONTOUR.
%
% TRICONTOUR can also return 3D contours similar to CONTOUR3 by adding
% view(3) after the call to TRICONTOUR.
%
% Type "contourdemo" for some examples.
%
% See also, CONTOUR, CLABEL
% This function does NOT interpolate back onto a Cartesian grid, but
% instead uses the triangulation directly.
%
% If your going to use this inside a loop with the same [p,t] a good
% modification is to make the connectivity "mkcon" once outside the loop
% because usually about 50% of the time is spent in "mkcon".
%
% Darren Engwirda - 2005 (d_engwirda@hotmail.com)
% Updated 15/05/2006
% I/O checking
if nargin~=4
    error('Incorrect number of inputs')
end
if nargout>2
    error('Incorrect number of outputs')
end
% Error checking
if (size(p,2)~=2) || (size(t,2)~=3) || (size(Hn,2)~=1)
    error('Incorrect input dimensions')
end
if size(p,1)~=size(Hn,1)
    error('F and p must be the same length')
end
if (max(t(:))>size(p,1)) || (min(t(:))<=0)
    error('t is not a valid triangulation of p')
end
if (size(N,1)>1) && (size(N,2)>1)
    error('N cannot be a matrix')
end
% Make mesh connectivity data structures (edge based pointers)
[e,eINt,e2t] = mkcon(p,t);
numt = size(t,1);       % Num triangles
nume = size(e,1);       % Num edges
%==========================================================================
%                Quadratic interpolation to centroids
%==========================================================================
% Nodes
t1 = t(:,1); t2 = t(:,2); t3 = t(:,3);
% FORM FEM GRADIENTS
% Evaluate centroidal gradients (piecewise-linear interpolants)
x23 = p(t2,1)-p(t3,1);  y23 = p(t2,2)-p(t3,2);
x21 = p(t2,1)-p(t1,1);  y21 = p(t2,2)-p(t1,2);
% Centroidal values
Htx = (y23.*Hn(t1) + (y21-y23).*Hn(t2) - y21.*Hn(t3)) ./ (x23.*y21-x21.*y23);
Hty = (x23.*Hn(t1) + (x21-x23).*Hn(t2) - x21.*Hn(t3)) ./ (y23.*x21-y21.*x23);
% Form nodal gradients.
% Take the average of the neighbouring centroidal values
Hnx = 0*Hn; Hny = Hnx; count = Hnx;
for k = 1:numt
    % Nodes
    n1 = t1(k); n2 = t2(k); n3 = t3(k);
    % Current values
    Hx = Htx(k); Hy = Hty(k);
    % Average to n1
    Hnx(n1)   = Hnx(n1)+Hx;
    Hny(n1)   = Hny(n1)+Hy;
    count(n1) = count(n1)+1;
    % Average to n2
    Hnx(n2)   = Hnx(n2)+Hx;
    Hny(n2)   = Hny(n2)+Hy;
    count(n2) = count(n2)+1;
    % Average to n3
    Hnx(n3)   = Hnx(n3)+Hx;
    Hny(n3)   = Hny(n3)+Hy;
    count(n3) = count(n3)+1;
end
Hnx = Hnx./count;
Hny = Hny./count;
% Centroids [x,y]
pt = (p(t1,:)+p(t2,:)+p(t3,:))/3;
% Take unweighted average of the linear extrapolation from nodes to centroids
Ht = ( Hn(t1) + (pt(:,1)-p(t1,1)).*Hnx(t1) + (pt(:,2)-p(t1,2)).*Hny(t1) + ...
       Hn(t2) + (pt(:,1)-p(t2,1)).*Hnx(t2) + (pt(:,2)-p(t2,2)).*Hny(t2) + ...
       Hn(t3) + (pt(:,1)-p(t3,1)).*Hnx(t3) + (pt(:,2)-p(t3,2)).*Hny(t3) )/3;
% DEAL WITH CONTOURING LEVELS
if length(N)==1
    lev = linspace(max(Ht),min(Ht),N+1);
    num = N;
else
    if (length(N)==2) && (N(1)==N(2))
        lev = N(1);
        num = 1;
    else
        lev = sort(N);
        num = length(N);
        lev = lev(num:-1:1);
    end
end
% MAIN LOOP
c   = [];
h   = [];
in  = false(numt,1);
vec = 1:numt;
old = in;
for v = 1:num       % Loop over contouring levels

    % Find centroid values >= current level
    i     = vec(Ht>=lev(v));
    i     = i(~old(i));         % Don't need to check triangles from higher levels
    in(i) = true;

    % Locate boundary edges in group
    bnd  = [i; i; i];       % Just to alloc
    next = 1;
    for k = 1:length(i)
        ct    = i(k);
        count = 0;
        for q = 1:3     % Loop through edges in ct
            ce = eINt(ct,q);
            if ~in(e2t(ce,1)) || ((e2t(ce,2)>0)&&~in(e2t(ce,2)))
                bnd(next) = ce;     % Found bnd edge
                next      = next+1;
            else
                count = count+1;    % Count number of non-bnd edges in ct
            end
        end
        if count==3                 % If 3 non-bnd edges ct must be in middle of group
            old(ct) = true;         % & doesn't need to be checked for the next level
        end
    end
    numb = next-1; bnd(next:end) = [];

    % Skip to next lev if empty
    if numb==0
        continue
    end

    % Place nodes approximately on contours by interpolating across bnd
    % edges
    t1  = e2t(bnd,1);
    t2  = e2t(bnd,2);
    ok  = t2>0;

    % Get two points for interpolation. Always use t1 centroid and
    % use t2 centroid for internal edges and bnd midpoint for boundary
    % edges

    % 1st point is always t1 centroid
    H1 = Ht(t1);                                                % Centroid value
    p1 = ( p(t(t1,1),:)+p(t(t1,2),:)+p(t(t1,3),:) )/3;          % Centroid [x,y]

    % 2nd point is either t2 centroid or bnd edge midpoint
    i1        = t2(ok);                                         % Temp indexing
    i2        = bnd(~ok);
    H2        = H1;
    H2(ok)    = Ht(i1);                                         % Centroid values internally
    H2(~ok)   = ( Hn(e(i2,1))+Hn(e(i2,2)) )/2;                  % Edge values at boundary
    p2        = p1;
    p2(ok,:)  = ( p(t(i1,1),:)+p(t(i1,2),:)+p(t(i1,3),:) )/3;   % Centroid [x,y] internally
    p2(~ok,:) = ( p(e(i2,1),:)+p(e(i2,2),:) )/2;                % Edge [x,y] at boundary

    % Linear interpolation
    r     = (lev(v)-H1)./(H2-H1);
    penew = p1 + [r,r].*(p2-p1);

    % Do a temp connection between adjusted node & endpoint nodes in
    % ce so that the connectivity between neighbouring adjusted nodes
    % can be determined
    vecb    = (1:numb)';
    m       = 2*vecb-1;
    c1      = 0*m;
    c2      = 0*m;
    c1(m)   = e(bnd,1);
    c1(m+1) = e(bnd,2);
    c2(m)   = vecb;
    c2(m+1) = vecb;

    % Sort connectivity to place connected edges in sucessive rows
    [c1,i] = sort(c1); c2 = c2(i);

    % Connect adjacent adjusted nodes
    k    = 1;
    next = 1;
    while k<(2*numb)
        if c1(k)==c1(k+1)
            c1(next) = c2(k);
            c2(next) = c2(k+1);
            next     = next+1;
            k        = k+2;         % Skip over connected edge
        else
            k = k+1;                % Node has only 1 connection - will be picked up above
        end
    end
    ncc          = next-1;
    c1(next:end) = [];
    c2(next:end) = [];

    % Plot the contours
    % If an output is required, extra sorting of the
    % contours is necessary for CLABEL to work.
    if nargout>0

        % Form connectivity for the contour, connecting
        % its edges (rows in cc) with its vertices.
        ndx = repmat(1,nume,1);
        n2e = 0*penew;
        for k = 1:ncc
            % Vertices
            n1 = c1(k); n2 = c2(k);
            % Connectivity
            n2e(n1,ndx(n1)) = k; ndx(n1) = ndx(n1)+1;
            n2e(n2,ndx(n2)) = k; ndx(n2) = ndx(n2)+1;
        end
        bndn = n2e(:,2)==0;         % Boundary nodes
        bnde = bndn(c1)|bndn(c2);   % Boundary edges

        % Alloc some space
        tmpv = repmat(0,1,ncc);

        % Loop through the points at the current contour level (lev(v))
        % Try to assemble the CS data structure introduced in "contours.m"
        % so that clabel will work. Assemble CS by "walking" around each
        % subcontour segment contiguously.
        ce    = 1;
        start = ce;
        next  = 2;
        cn    = c2(1);
        flag  = false(ncc,1);
        x     = tmpv; x(1) = penew(c1(ce),1);
        y     = tmpv; y(1) = penew(c1(ce),2);
        for k = 1:ncc

            % Checked this edge
            flag(ce) = true;

            % Add vertices to patch data
            x(next) = penew(cn,1);
            y(next) = penew(cn,2);
            next    = next+1;

            % Find edge (that is not ce) joined to cn
            if ce==n2e(cn,1)
                ce = n2e(cn,2);
            else
                ce = n2e(cn,1);
            end

            % Check the new edge
            if (ce==0)||(ce==start)||(flag(ce))

                % Plot current subcontour as a patch and save handles
                x   = x(1:next-1);
                y   = y(1:next-1);
                z   = repmat(lev(v),1,next);
                h   = [h; patch('Xdata',[x,NaN],'Ydata',[y,NaN],'Zdata',z, ...
                                'Cdata',z,'facecolor','none','edgecolor','flat')]; hold on

                % Update the CS data structure as per "contours.m"
                % so that clabel works
                c = horzcat(c,[lev(v), x; next-1, y]);

                if all(flag)    % No more points at lev(v)
                    break
                else            % More points, but need to start a new subcontour

                    % Find the unflagged edges
                    edges = find(~flag);
                    ce    = edges(1);
                    % Try to select a boundary edge so that we are
                    % not repeatedly running into the boundary
                    for i = 1:length(edges)
                        if bnde(edges(i))
                            ce = edges(i); break
                        end
                    end
                    % Reset counters
                    start = ce;
                    next  = 2;
                    % Get the non bnd node in ce
                    if bndn(c2(ce))
                        cn = c1(ce);
                        % New patch vectors
                        x = tmpv; x(1) = penew(c2(ce),1);
                        y = tmpv; y(1) = penew(c2(ce),2);
                    else
                        cn = c2(ce);
                        % New patch vectors
                        x = tmpv; x(1) = penew(c1(ce),1);
                        y = tmpv; y(1) = penew(c1(ce),2);
                    end

                end

            else
                % Find node (that is not cn) in ce
                if cn==c1(ce)
                    cn = c2(ce);
                else
                    cn = c1(ce);
                end
            end

        end

    else        % Just plot the contours as is, this is faster...

        z = repmat(lev(v),2,ncc);

        patch('Xdata',[penew(c1,1),penew(c2,1)]', ...
              'Ydata',[penew(c1,2),penew(c2,2)]', ...
              'Zdata',z,'Cdata',z,'facecolor','none','edgecolor','flat');

        hold on

    end

end
% Assign outputs if needed
if nargout>0
    cout = c;
    hout = h;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [e,eINt,e2t] = mkcon(p,t)
numt = size(t,1);
vect = 1:numt;
% DETERMINE UNIQUE EDGES IN MESH

e       = [t(:,[1,2]); t(:,[2,3]); t(:,[3,1])];             % Edges - not unique
vec     = (1:size(e,1))';                                   % List of edge numbers
[e,j,j] = unique(sort(e,2),'rows');                         % Unique edges
vec     = vec(j);                                           % Unique edge numbers
eINt    = [vec(vect), vec(vect+numt), vec(vect+2*numt)];    % Unique edges in each triangle
% DETERMINE EDGE TO TRIANGLE CONNECTIVITY
% Each row has two entries corresponding to the triangle numbers
% associated with each edge. Boundary edges have one entry = 0.
nume = size(e,1);
e2t  = repmat(0,nume,2);
ndx  = repmat(1,nume,1);
for k = 1:numt
    % Edge in kth triangle
    e1 = eINt(k,1); e2 = eINt(k,2); e3 = eINt(k,3);
    % Edge 1
    e2t(e1,ndx(e1)) = k; ndx(e1) = ndx(e1)+1;
    % Edge 2
    e2t(e2,ndx(e2)) = k; ndx(e2) = ndx(e2)+1;
    % Edge 3
    e2t(e3,ndx(e3)) = k; ndx(e3) = ndx(e3)+1;
end
return

repo: ifem


text:
function y = end( x, k, n )

%Disciplined convex/geometric programming information for END:
%   The use of END as an array subscript (e.g., X(:,end)) is identical
%   with CVX variables as it is for numeric vectors and arrays.

sz = size( x );
nz = length( sz );
if k > nz,
    y = 1;
elseif k < n || nz <= n,
    y = sz( k );
else
    y = prod( sz( k : end ) );
end

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

author: yu-jiang

repo: radpbook

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


text:
classdef Skewness < Algorithm

    methods (Access = public)

        function obj = Skewness()
            obj.name = 'Skewness';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)
            %meanAlgorithm = dsp.Mean('RunningMean',false,'Dimension','All');
            %signalMean = step(meanAlgorithm,signal);
            signalMean = mean(signal);

            n = single(length(signal));

            upper = single(0);
            lower = single(0);
            for i = 1 : n
                temp = single(signal(i) - signalMean);
                temp2 = single(temp * temp);
                upper = upper + temp2 * temp;
                lower = lower + temp2;
            end
            lower = lower / n;
            result = upper / (n * lower^1.5);
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 6 * n;
            memory = 1;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end

end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function this = read_nifti_analyze(this, fileName, selectedVolumes)
% loads matrix into .data from nifti or analyze file using spm_read_vols
%
%   this = read_nifti_analyze(this, fileName, selectedVolumes)
%
% IN
%   fileName
%   selectedVolumes     [1,nVols] index of selected volumes to load,
%                       Inf for all; default: Inf
% OUT
%
% EXAMPLE
%   read_nifti_analyze
%
%   See also

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-04-16
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 2
    fileName = fullfile(this.parameters.path, ...
        this.parameters.unprocessedFile);
end

fileNameVolArray = tapas_uniqc_get_vol_filenames(fileName);

V = spm_vol(strvcat(fileNameVolArray));

hasSelectedVolumes = nargin > 2 && ~any(isinf(selectedVolumes));

if hasSelectedVolumes
    V = V(selectedVolumes);
end

% more than 4 dimensions not handles by spm_vol, read directly from
% file_array
if numel(V(1).private.dat.dim) > 4
    this.data = reshape(V.private.dat(:), size(V.private.dat));
else
    % use inbuilt SPM functionality
    try
        % this.data = transform_matrix_analyze2matlab(spm_read_vols(V));
        this.data = spm_read_vols(V);
        % maybe only header misalignment of volumes is the problem for nifti
        %...rename temporarily for loading
    catch err
        fnHdr = regexprep(fileName, '\.nii','\.mat');
        fnTmp = regexprep(fileName, '\.nii','\.tmp');
        if exist(fnHdr, 'file')
            movefile(fnHdr, fnTmp);
            warning('Headers of volumes not aligned, ignoring them...');
            V = spm_vol(fileName);
            this.data = spm_read_vols(V);
            movefile(fnTmp, fnHdr);
        else % nothing we can do, throw error
            throw(err);
        end
    end
end

author: translationalneuromodeling

repo: tapas


text:
function [params] = pt_calc(params)
    % Calculation of the probabilistic forecast test.
    %
    % [params] = pt_calc(params)
    %
    % Input parameters:
    %   params.mCatalog           Earthquake catalog
    %   params.mPolygon           Polygon (defined by ex_selectgrid)
    %   params.vX                 X-vector (defined by ex_selectgrid)
    %   params.vY                 Y-vector (defined by ex_selectgrid)
    %   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
    %   params.bRandom            Perform random simulation (true) or real calculation (false)
    %   params.nCalculation       Number of random simulations
    %   params.bMap               Calculate a map (true) or a cross-section (false)
    %   params.bNumber            Use constant number (true) or constant radius (false)
    %   params.nNumberEvents      Number of earthquakes if bNumber is true
    %   params.fRadius            Radius of gridnode if bNumber is false
    %   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
    %   params.fForecastPeriod    Forecasting period in years
    %   params.bLearningPeriod    Use different learning period than the rest of the catalog
    %   params.fLearningPeriod    Learning period in years
    %   params.bSignificance      Calculate significance during random simulation
    %                             using params.fRealProbability
    %   params.fRealProbability   Probability of calculation with real data
    %   params.nCalculateMC       Method to calculate magnitude of completeness (see also: help calc_Mc)
    %   params.nTestMethod        Method to calculate the Kagan & Jackson test (see also: help kj_poissonian)
    %   params.bMinMagMc          Use magnitude of completeness as lower limit of magnitude range for testing (true)
    %                             Use params.fMinMag as lower limit (false)
    %   params.fMinMag            Lower limit of magnitude range for testing
    %   params.fMaxMag            Upper limit of magnitude range for testing
    %
    % Output parameters:
    %   Same as input parameters including
    %   params.mValueGrid         Matrix of calculated Kagan & Jackson test values
    %   params.vRandomMeans       Vector of means of probability differences per simulation run
    %   params.vSignificanceLevel Vector of significance levels per simulation run
    %   params.fBValueOverall     Calculated overall b-value
    %   params.fStdDevOverall     Calculated standard deviation
    %   params.fMcOverall         Calculated magnitude of completeness
    %
    % Danijel Schorlemmer
    % July 10, 2002

    report_this_filefun();

    % Perform the calculation
    % -----------------------
    % Compute the b-values for both models
    [mValueGridH, vcsGridNamesH] = pt_calcinput(params.mLearningCatalog, params.mObservedCatalog, params.mPolygon, params.rOptions(1));
    [mValueGridN, vcsGridNamesN] = pt_calcinput(params.mLearningCatalog, params.mObservedCatalog, params.mPolygon, params.rOptions(2));
    % Apply test moethod settings
    if params.nTestMethod == 2    % Use overall b-value for test hypothesis if it's impossible to compute a b-value
        vSel = isnan(mValueGridH(:,1));
        mValueGridH(vSel,:) = mValueGridN(vSel,:);
    end
    % Init result matrix
    mValueGrid_ = [];
    % Loop over all grid nodes
    for nNode_ = 1:length(params.mPolygon(:,1))
        % Create node catalogs
        mLearningNodeCatalog_ = params.mLearningCatalog(params.rOptions(3).caLearningNodeIndices{nNode_}, :);
        mObservedNodeCatalog_ = params.mObservedCatalog(params.rOptions(3).caObservedNodeIndices{nNode_}, :);
        % Define magnitude range for testing
        fMinMag_ = max([mValueGridH(nNode_,3) mValueGridN(nNode_,3)]);
        if ~(params.bMinMagMc)
            fMinMag_ = max(params.fMinMag, fMinMag_);
        end
        % Calculate the probability-ratio
        [fDeltaProbability, fProbabilityN, fProbabilityH, vPredictionFMD, vObservedFMD, vMagnitudeBins] = pt_poissonian(mLearningNodeCatalog_, params.fLearningPeriodUsed, ...
            mObservedNodeCatalog_, params.fObservedPeriodUsed, params.rOptions(3).nMinimumNumber, mValueGridH(nNode_,1), mValueGridN(nNode_,1), ...
            mValueGridH(nNode_,3), mValueGridN(nNode_,3), fMinMag_, params.fMaxMag);
        if (((params.bRandomNode) | (params.bSaveRates)) & (~isnan(fProbabilityH)))
            nLen_ = length(vObservedFMD(:,1));
            vXMin_ = ones(nLen_, 1) * params.mPolygon(nNode_,1) - (params.rOptions(3).fSizeRectX/2);
            vXMax_ = ones(nLen_, 1) * params.mPolygon(nNode_,1) + (params.rOptions(3).fSizeRectX/2);
            vYMin_ = ones(nLen_, 1) * params.mPolygon(nNode_,2) - (params.rOptions(3).fSizeRectY/2);
            vYMax_ = ones(nLen_, 1) * params.mPolygon(nNode_,2) + (params.rOptions(3).fSizeRectY/2);
            vZMin_ = zeros(nLen_, 1);
            vZMax_ = zeros(nLen_, 1);
            vWeight_ = ones(nLen_, 1);
            vMagMin_ = vMagnitudeBins - 0.05;
            vMagMax_ = vMagnitudeBins + 0.05;
            vRatesH = [vXMin_ vXMax_ vYMin_ vYMax_ vZMin_ vZMax_ vMagMin_ vMagMax_ vPredictionFMD(:,1) vWeight_ vObservedFMD];
            vRatesN = [vXMin_ vXMax_ vYMin_ vYMax_ vZMin_ vZMax_ vMagMin_ vMagMax_ vPredictionFMD(:,2) vWeight_ vObservedFMD];
            if params.bSaveRates
                params.vRatesH = [params.vRatesH; vRatesH];
                params.vRatesN = [params.vRatesN; vRatesN];
            end
            if params.bRandomNode
                [rRelmTest] = relm_RTest4(vRatesH, vRatesN, params.nNumberCalculationNode, vMagnitudeBins(1), 1, 1, 0);
                fAlpha = rRelmTest.fAlpha;
                fBeta = rRelmTest.fBeta;
            else
                fAlpha = nan;
                fBeta = nan;
            end
        else
            fAlpha = nan;
            fBeta = nan;
        end
        mValueGrid_= [mValueGrid_; fDeltaProbability fProbabilityN fProbabilityH ...
            fAlpha fBeta];
    end % for nNode
    params.vcsGridNames = cellstr(char('Probability difference', ...
        'Null hypothesis', 'Test hypothesis', 'Alpha', 'Beta'));
    % Add the b-values of the hypothesis to the overall value grid
    params.mValueGrid = [mValueGrid_ mValueGridH mValueGridN];
    params.vcsGridNames = [params.vcsGridNames; vcsGridNamesH; vcsGridNamesN];
end

author: CelsoReyes


text:
function [Ain, Cin, b_in, f_in, center, res] = greedyROI(Y, K, params, ROI_list)
% component initialization using a greedy algorithm to identify neurons in 2d or 3d calcium imaging movies
%
% Usage:     [Ain, Cin, bin, fin, center, res] = greedyROI2d(data, K, params)
%
% Input:
% Y          d1 x d2 x (d3 x) T movie, raw data, each column is a vectorized movie
% K          number of neurons to extract (if K is a vector then the algorithm is run multiple times with different parameters)
% params     tuning parameter for fine-tuning the shape (optional)
%            params.nIter: number of iterations for shape tuning (default 5)
%            params.gSig: variance of Gaussian kernel to use (default 5) If params.gSig is a cell, then the algorithm is run with multiple times with different parameters
%            params.gSiz: size of kernel (default 2*gSiz+1)
%            params.nb: rank of background component (default 1)
%            params.save_memory: flag for processing data in chunks to save memory (default 0)
%            params.windowSiz: size of spatial window when computing the median (default 32 x 32)
%            params.chunkSiz: number of timesteps to be processed simultaneously if on save_memory mode (default: 100)
%            params.med_app: number of timesteps to be interleaved for fast (approximate) median calculation (default: 1, no approximation)
%            params.rolling_sum: flag for using rolling sum to detect new components (default: True)
%            params.rolling_length: length of rolling window (default: 100)
% ROI_list   Kn x 2 (or 3) matrix with user specified centroids. If this are present, then the algorithm only finds components around these centroids

%Output:
% Ain        (d) x K matrix, location of each neuron
% Cin        T x K matrix, calcium activity of each neuron
% center     K x 2 (or 3) matrix, inferred center of each neuron
% bin        (d) X nb matrix, initialization of spatial background
% fin        nb X T matrix, initalization of temporal background
% res        d1 x d2 x (d3 x) T movie, residual
%
% Author: Yuanjun Gao with modifications from Eftychios A. Pnevmatikakis and Weijian Yang

use_sum = false;
if nargin < 4 || isempty(ROI_list)
    user_ROIs = 0;
else
    user_ROIs = 1;
    K = size(ROI_list,1);
end

dimY = ndims(Y) - 1;  % dimensionality of imaged data (2d or 3d)
sizY = size(Y);
T = sizY(end);        % # of timesteps
dx = sizY(1:dimY);    % # of voxels in each axis
d = prod(dx);         % total # of voxels

if ~exist('K', 'var'),  %find K neurons
    K = 30;
    warning(['number of neurons are not specified, set to be the default value', num2str(K)]);
end

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

if ~isfield(params, 'gSig') || isempty(params.gSig);
    if dimY == 2; params.gSig = [5, 5]; else params.gSig = [5,5,5]; end
elseif length(params.gSig) == 1, params.gSig = params.gSig + zeros(1,dimY);
    if dimY == 3; params.gSig(3) = params.gSig(3)/2; end
end

if ~isfield(params, 'gSiz') || isempty(params.gSiz);
    if ~iscell(params.gSig)
        params.gSiz = ceil(2*params.gSig + 1);
    else
        for j = 1:length(params.gSig)
            params.gSiz{j,1} = 2*params.gSig{j}+1; %cellfun(@times,params.gSig{j},num2cell(ones(size(params.gSig{j}))*2));
        end
    end
elseif length(params.gSiz) == 1, params.gSiz = params.gSiz + zeros(1,dimY);
    if dimY == 3; params.gSiz(3) = ceil(params.gSiz(3)/2); end
end

if isfield(params,'ssub');
    if ~iscell(params.gSig); params.gSig(1:2) = params.gSig(1:2)/params.ssub; params.gSiz(1:2) = ceil(params.gSiz(1:2)/params.ssub);
    else
        for j = 1:length(params.gSig)
            params.gSig{j,1} = params.gSig{j}/params.ssub; %cellfun(@times,params.gSig{j},num2cell(ones(size(params.gSig{j}))/params.ssub));
            params.gSiz{j,1} = params.gSiz{j}/params.ssub; %cellfun(@times,params.gSiz{j},num2cell(ones(size(params.gSiz{j}))/params.ssub));
        end
    end
end

if ~isfield(params,'nb'), nb = 1; else nb = params.nb; end

if ~isfield(params, 'nIter'), nIter = 5;
else nIter = params.nIter; end

if ~isfield(params, 'save_memory'), save_memory = 0;
else save_memory = params.save_memory; end

if ~isfield(params, 'chunkSiz'), chunkSiz = 100; else chunkSiz = params.chunkSiz; end
if ~isfield(params, 'windowSiz'), windowSiz = 32; else windowSiz = params.windowSiz; end
if ~isfield(params, 'med_app'), med_app = 1; else med_app = params.med_app; end

if ~isfield(params,'rem_prct') || isempty(params.rem_prct); params.rem_prct = 20; end

Tint = 1:med_app:T;
% if save_memory
%     med = zeros(M,N);
%     for ii = 1:ceil(M/windowSiz)
%         intx = (ii-1)*windowSiz+1:min(ii*windowSiz,M);
%         for jj = 1:ceil(N/windowSiz)
%             inty = (jj-1)*windowSiz+1:min(jj*windowSiz,N);
%             med(intx,inty) = median(data(intx,inty,Tint),3);
%         end
%     end
% else
%if dimY == 2; med = median(Y(:,:,Tint), 3); else med = median(Y(:,:,:,Tint), 4); end
if dimY == 2; med = prctile(Y(:,:,Tint),params.rem_prct, 3); else med = prctile(Y(:,:,:,Tint),params.rem_prct,4); end

% end

Y = bsxfun(@minus, Y, med);
if iscell(params.gSig); params.gSig = cell2mat(params.gSig); params.gSiz = cell2mat(params.gSiz); end
if length(K) > 1  % order size of components to be found in descending order
    [~,ord] = sort(sum(params.gSig,2),'descend');
    K = K(ord);
    params.gSig = params.gSig(ord,:);
    params.gSiz = params.gSiz(ord,:);
end

if ~iscell(params.gSiz)
    Ain = spalloc(d,sum(K),K(:)'*ceil(params.gSiz(:,1).^2)); %zeros(M*N,sum(K));
else
    Ain = sparse(d,sum(K));
end
Cin = zeros(sum(K),T);
center = zeros(sum(K),dimY);

for r = 1:length(K)
    gSig = params.gSig(r,:);
    gSiz = params.gSiz(r,:);

    gHalf = floor(gSiz / 2); %half size of the kernel, used to calculate margin
    gSiz = 2 * gHalf + 1; %actual size

    Atemp = spalloc(d,K(r),K(r)*ceil(prod(gSiz))); %zeros(M, N, K(r));
    %basis = spalloc(M*N,K,K*prod(gSiz));
    trace = zeros(T, K(r));
    centers = zeros(K(r), dimY);

    %scan the whole image (only need to do this at the first iteration)
    rho = imblur(Y, gSig, gSiz, dimY, save_memory, chunkSiz); %covariance of data and basis

    if ~params.rolling_sum
        v = sum(rho.^2, dimY+1); %variance explained
    else
        % running maximum
        avg_fil = ones(1,params.rolling_length)/params.rolling_length;
        rho_s = filter(avg_fil,1,rho.^2,[],dimY+1);
        v = max(rho_s,[],dimY+1);
    end

    for k = 1:K(r),
        if user_ROIs
            iHat = ROI_list(k,:);
        else
            [~, ind] = max(v(:));
            %[iHat, jHat] = ind2sub([M, N], ind);
            iHat = zeros(1,dimY);
            if dimY == 2; [iHat(1),iHat(2)] = ind2sub(dx,ind); else [iHat(1),iHat(2),iHat(3)] = ind2sub(dx,ind); end
        end
        centers(k,:) = round(iHat);

        %iSig = [max(iHat - gHalf(1), 1), min(iHat + gHalf(1), M)]; iSigLen = iSig(2) - iSig(1) + 1;
        %jSig = [max(jHat - gHalf(2), 1), min(jHat + gHalf(2), N)]; jSigLen = jSig(2) - jSig(1) + 1;
        iSig = [max(iHat - gHalf,1)', min(iHat + gHalf,dx)']; iSigLen = iSig(:,2) - iSig(:,1) + 1;

        %fine tune the shape
        if dimY == 2
            dataTemp = Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), :);
            traceTemp = squeeze(rho(iHat(1), iHat(2), :));
        else
            dataTemp = Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), iSig(3,1):iSig(3,2), :);
            traceTemp = squeeze(rho(iHat(1), iHat(2), iHat(3), :));
        end
        [coef, score] = finetune(dataTemp, traceTemp, nIter);

        dataSig = bsxfun(@times, coef, reshape(score, [ones(1,dimY),T]));
        %[xSig,ySig] = meshgrid(iSig(1):iSig(2),jSig(1):jSig(2));
        basis = zeros(dx);
        if dimY == 2;  basis(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2)) = coef; else  basis(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), iSig(3,1):iSig(3,2)) = coef; end
        Atemp(:,k) = basis(:);
        %basis(sub2ind([M,N],xSig(:),ySig(:)),k) = coef(:);

        trace(:, k) = score';

        if dimY == 2;
            Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), :) = Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), :) - dataSig; %update residual
        else
            Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2), iSig(3,1):iSig(3,2), :) = Y(iSig(1,1):iSig(1,2), iSig(2,1):iSig(2,2),  iSig(3,1):iSig(3,2),:) - dataSig;
        end

        if mod(k,10) == 0; fprintf('found %i out of %i neurons..\n', k,K(r)); end

        %get next basis;
        if k < K(r)
            %iMod = [max(iHat - 2 * gHalf(1), 1), min(iHat + 2 * gHalf(1), M)]; iModLen = iMod(2) - iMod(1) + 1;%patches to modify
            %jMod = [max(jHat - 2 * gHalf(2), 1), min(jHat + 2 * gHalf(2), N)]; jModLen = jMod(2) - jMod(1) + 1;
            iMod = [max(iHat - 2*gHalf,1)', min(iHat + 2*gHalf,dx)']; iModLen = iMod(:,2) - iMod(:,1) + 1;
            %iLag = iSig - iMod(1) + 1; %relative location of iSig in the small patch
            %jLag = jSig - jMod(1) + 1;
            iLag = iSig - repmat(iMod(:,1),1,2) + 1;
            %dataTemp = zeros(iModLen, jModLen);
            dataTemp = zeros(iModLen(:)');
            %dataTemp(iLag(1):iLag(2), jLag(1):jLag(2)) = reshape(coef, [iSigLen, jSigLen]);
            if dimY == 2
                dataTemp(iLag(1,1):iLag(1,2), iLag(2,1):iLag(2,2)) = reshape(coef, iSigLen(:)');
            else
                dataTemp(iLag(1,1):iLag(1,2), iLag(2,1):iLag(2,2), iLag(3,1):iLag(3,2)) = reshape(coef, iSigLen(:)');
            end
            dataTemp = imblur(dataTemp, gSig, gSiz, dimY, 0, chunkSiz);
            rhoTemp = bsxfun(@times, dataTemp, reshape(score, [ones(1,dimY),T]));
            if dimY == 2
                rhoTemp = rho(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), :)  - rhoTemp;
                rho(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), :)  = rhoTemp;
                if ~params.rolling_sum
                    v(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2)) = sum(rhoTemp.^2, 3);
                else
                    rho_filt = filter(avg_fil,1,rhoTemp.^2,[],3);
                    v(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2)) = max(rho_filt,[],3);
                end
            else
                rhoTemp = rho(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), iMod(3,1):iMod(3,2), :)  - rhoTemp;
                rho(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), iMod(3,1):iMod(3,2), :)  = rhoTemp;
                if ~params.rolling_sum
                    v(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), iMod(3,1):iMod(3,2)) = sum(rhoTemp.^2, 4);
                else
                    rho_filt = filter(avg_fil,1,rhoTemp.^2,[],4);
                    v(iMod(1,1):iMod(1,2), iMod(2,1):iMod(2,2), iMod(3,1):iMod(3,2)) = max(rho_filt,[],4);
                end
            end
        end
    end
    Ain(:,sum(K(1:r-1))+1:sum(K(1:r))) = Atemp; %sparse(reshape(basis,d,K(r)));
    Cin(sum(K(1:r-1))+1:sum(K(1:r)),:) = trace';
    center(sum(K(1:r-1))+1:sum(K(1:r)),:) = centers;
end
res = reshape(Y,d,T) + repmat(med(:),1,T);

%% clear data matrix from local memory (avoid out-of-memory? see greedyROI_corr.m)
clear Y

%[b_in,f_in] = nnmf(max(res,0),nb);
%[b_in,f_in] = nnsvd(max(res,0),nb);
f_in = [mean(res);rand(nb-1,T)];

for nmfiter = 1:100
    b_in = max((res*f_in')/(f_in*f_in'),0);
    f_in = max((b_in'*b_in)\(b_in'*res),0);
end

function [ain,cin] = finetune(data,cin,nIter)

    %rank-1 semi-NMF to fine tune the inferred components
    %
    %Input:
    %data   d1 x d2 x (d3 x) T matrix, small patch containing one neuron
    %trace  initial value for trace
    %nIter  number of coordinate descent steps
    %
    %Output:
    %basis  d1 x d2 (x d3) matrix, result of the fine-tuned neuron shape
    %trace  1 x T matrix, result of the neuron

    if ~exist('nIter', 'var'), nIter = 1; end
    data = reshape(data,prod(iSigLen),T);
    for iter = 1:nIter
        %nc = norm(cin)^2;
        ain = max(data*cin,0)/norm(cin);
        an = norm(ain);
        if an > 0
            ain = ain/an;
        else
            fprintf('found degenerate component!\n')

            %break
        end
        cin = data'*ain;
    end
    ain = reshape(ain,iSigLen(:)');
    cin = cin(:)';
end

function data = imblur(data, sig, siz, nDimBlur, save_memory, chunkSiz)
%Gaussian blur for high dimensional data
%Input:
%data       original data
%sig        std of gaussian kernel
%siz        size of kernel
%nDimBlur   number of dims to blur (default: ndims(data) - 1)
%
%Output:
%data       result after the Gaussian blur

if ~exist('nDimBlur', 'var'), nDimBlur = ndims(data) - 1;
else nDimBlur = min(nDimBlur, ndims(data)); end

if length(sig) == 1, sig = sig * ones(1,nDimBlur); end
assert(nDimBlur == length(sig));

if length(siz) == 1, siz = siz * ones(1,nDimBlur); end
assert(nDimBlur == length(siz));

for i = 1:nDimBlur,
    if sig(i) > 0,
        x = -floor(siz(i) / 2):floor(siz(i) / 2);
        H = exp(-(x.^2/ (2 * sig(i)^2)));
        H = H' / norm(H(:));
        if nDimBlur > 1,
            indH = 1:nDimBlur; indH(i) = 1; indH(1) = i;
            H = permute(H, indH);
        end
        if save_memory
            L = size(data,ndims(data));
            for ci = 1:ceil(L/chunkSiz)
                int = (ci-1)*chunkSiz+1:min(ci*chunkSiz,L);
                data(:,:,int) = imfilter(data(:,:,int), H, 'same', 0);
            end
        else
            data = imfilter(data, H, 'same', 0);
        end
    end
end
end

end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function varargout = DSeg(varargin)
% DSEG MATLAB code for DSeg.fig
%      DSEG, by itself, creates a new DSEG or raises the existing
%      singleton*.
%
%      H = DSEG returns the handle to a new DSEG or the handle to
%      the existing singleton*.
%
%      DSEG('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DSEG.M with the given input arguments.
%
%      DSEG('Property','Value',...) creates a new DSEG or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before DSeg_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to DSeg_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help DSeg

% Last Modified by GUIDE v2.5 14-Jun-2011 16:56:58

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @DSeg_OpeningFcn, ...
                   'gui_OutputFcn',  @DSeg_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before DSeg is made visible.
function DSeg_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to DSeg (see VARARGIN)

% Choose default command line output for DSeg
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes DSeg wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = DSeg_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile({'*.*'},'Open Directory');
if isequal(filename,0) || isequal(pathname,0)
    return
end

set(handles.edit1,'String',strcat(pathname,filename));

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double

% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Im=get(handles.edit1,'String');
if length(Im)>3
flag=str2double(get(handles.edit2,'String'));
open=str2double(get(handles.edit3,'String'));
EdHist=str2double(get(handles.edit4,'String'));

[Segmented, Array]=DS(Im,flag,open,EdHist);
imwrite(Segmented,'Segmented.tif')
set(handles.listbox1,'String',num2str(Array));
end
% --- Executes on selection change in listbox1.
function listbox1_Callback(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns listbox1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from listbox1

% --- Executes during object creation, after setting all properties.
function listbox1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Sable

repo: mcbench-benchmarks


text:
function w = MakeSineWindow(boxlen,sigma)
% MakeSineWindow -- Make a nice window out of the IteratedSine function.
%  Inputs
%    boxlen   abscissa values for window evaluation -boxlen <= t < boxlen
%    sigma    scale of the window, controls the size of the
%             effective support
%  Outputs
%    w        values of the window
%  See Also
%    IteratedSine
%
% By Emmanuel Candes, 2003-2004

  ix = ((-boxlen:(boxlen -1)) + .5)./boxlen;
  w = IteratedSineWindow(ix./sigma);

author: thomaskuestner


text:
function g = translateKernGradient(kern, varargin)

% TRANSLATEKERNGRADIENT Gradient of TRANSLATE kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% input space translation
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO translateKernParamInit, kernGradient,
% translateKernDiagGradient, cmpndKernGradient, cmpndKernGradX, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

for i = 1:length(varargin)-1
  varargin{i} = varargin{i} - repmat(kern.centre, size(varargin{i}, 1), 1);
end
g = cmpndKernGradient(kern, varargin{:});

% Gradients with respect to X.
if length(varargin) == 2
  gKX = cmpndKernGradX(kern, varargin{1}, varargin{1});
  gKX = gKX*2;
  dgKX = cmpndKernDiagGradX(kern, varargin{1});
  for i = 1:size(gKX, 1)
    gKX(i, :, i) = dgKX(i, :);
  end
else
  gKX_12 = cmpndKernGradX(kern, varargin{1}, varargin{2});
  gKX_21 = cmpndKernGradX(kern, varargin{2}, varargin{1});
end

gcentre = zeros(1, kern.inputDimension);
if length(varargin) == 2
  for i = 1:size(varargin{1}, 1);
    for j = 1:kern.inputDimension
      gcentre(1, j) = gcentre(1, j) - gKX(:, j, i)'*varargin{end}(:, i);
    end
  end
else
  for i = 1:size(varargin{1}, 1)
    for j = 1:kern.inputDimension
      gcentre(1, j) = gcentre(1, j) - sum(gKX_12(:, j, i).*varargin{end}(i, :)');
    end
  end
  for i = 1:size(varargin{2}, 1)
    for j = 1:kern.inputDimension
      gcentre(1, j) = gcentre(1, j) - sum(gKX_21(:, j, i).*varargin{end}(:, i));
    end
  end
end
g = [g gcentre];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [mmHg] = kPa2mmHg(kPa)
% Convert pressure from kilopascals to millimeters of mercury.
% Chad Greene 2012
mmHg = kPa*7.50062;


text:
function ind = inpolygon(gb,xy)
% checks which gb data are within given polygon
%
% Syntax
%   ind = inpolygon(gb,[xmin,ymin,dx,dy]) % select indices by rectangle
%   ind = inpolygon(gb,[x1 y1; x2 y2; x3 y3; x4 y4]) % select indices by poylgon
%   gb = gb(ind) % select gb data by indices
%
% Input
%  gb    - @grainBoundary
%  xmin, xmax - lower left corner of a rectangle
%  dx, dy - extent of a rectangle
%  x, y  - vertices of a polygon
%
% Output
%  ind - logical
%
% See also
% inpolygon

% shortcut for simple rectangles
if numel(xy)==4

  xy(3:4) = xy(1:2) + xy(3:4);
  corners = [1 2; 3 2; 3 4; 1 4; 1 2];
  xy = xy(corners);

end

if ~isempty(gb.midPoint(:,1))
  %  check for inside
  ind = inpolygon(gb.midPoint(:,1),gb.midPoint(:,2),xy(:,1),xy(:,2));

end

author: mtex-toolbox

repo: mtex

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


text:
function [sameReactions,similarReactions] = mapKBaseToVMHReactions(translatedRxns)
% Part of the DEMETER pipeline. This functions aids in translating
% reactions from KBase to VMH nomenclature. Requires running the function
% propagateKBaseMetTranslationToRxns beforehand to translate metabolite
% IDs, which will then allow matching translated reactions to reactions
% that already exist in the VMH (Virtual Metabolic Human) database.
%
% USAGE:
%
%           [sameReactions,similarReactions] = mapKBaseToVMHReactions(translatedRxns)
%
% INPUTS
% translatedRxns:       Table with untranslated KBase reactions but
%                       translated metabolite IDs
%
% OUTPUT
% sameReactions:        Table with translated KBase reactions that already
%                       exist in the VMH database with corresponding IDs
% similarReactions:     Table with translated KBase reactions for which a
%                       reaction with the same formula but irreversible in
%                       VMH and reversible in KBase (or vice versa) exists
%
% .. Authors:
%       - Almut Heinken, 01/2021

sameReactions={'KBase_reaction','VMH_reaction','KBase_Formula','VMH_Formula'};
similarReactions={'KBase_reaction','VMH_reaction','KBase_Formula','VMH_Formula'};

% get VMH reaction database
database=loadVMHDatabase;

parsedVMHFormulas={};
cnt=1;

% parse the VMH reaction database
for i=1:size(database.reactions,1)
    % break down the formula
    [metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(database.reactions{i,3});
    % now put the formula back together, ensuring it is written the same
    % way as translated KBase formulas
    if revFlag==1
        revChar=' <=> ';
    elseif revFlag==0
        revChar=' -> ';
    end
    rxnForm='';
    % parse substrates
    substrates=find(stoichCoeffList<0);
    if ~isempty(substrates)
        for j=1:length(substrates)
            if j<length(substrates)
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(substrates(j)))) ' ' metaboliteList{substrates(j)} ' + ' ];
            else
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(substrates(j)))) ' ' metaboliteList{substrates(j)} revChar ];
            end
        end
    else
      % catch cases of exchange reactions
      rxnForm=[rxnForm revChar];
    end

    % parse products
    products=find(stoichCoeffList>0);
    if ~isempty(products)
        for j=1:length(products)
            if j<length(products)
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(products(j)))) ' ' metaboliteList{products(j)} ' + ' ];
            else
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(products(j)))) ' ' metaboliteList{products(j)}];
            end
        end
    end

    parsedVMHFormulas{cnt,1}=database.reactions{i,1};
    parsedVMHFormulas{cnt,2}=rxnForm;
    parsedVMHFormulas{cnt,3}=revFlag;
    cnt=cnt+1;
    if revFlag==1
        % also catch the reaction being written the other way around
        rxnForm=strsplit(rxnForm,' <=> ');
        parsedVMHFormulas{cnt,1}=database.reactions{i,1};
        parsedVMHFormulas{cnt,2}=[rxnForm{2} ' <=> ' rxnForm{1}];
        parsedVMHFormulas{cnt,3}=revFlag;
        cnt=cnt+1;
    end
end

% parse KBase reactions to translate in the same manner
parsedKBaseFormulas={};
cnt=1;

for i=1:size(translatedRxns,1)
    % break down the formula
    [metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(translatedRxns{i,3});
    % now put the formula back together, ensuring it is written the same
    % way as translated VMH formulas
    if revFlag==1
        revChar=' <=> ';
    elseif revFlag==0
        revChar=' -> ';
    end
    rxnForm='';
    % parse substrates
    substrates=find(stoichCoeffList<0);
    if ~isempty(substrates)
        for j=1:length(substrates)
            if j<length(substrates)
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(substrates(j)))) ' ' metaboliteList{substrates(j)} ' + ' ];
            else
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(substrates(j)))) ' ' metaboliteList{substrates(j)} revChar ];
            end
        end
    else
      % catch cases of exchange reactions
      rxnForm=[rxnForm revChar];
    end

    % parse products
    products=find(stoichCoeffList>0);
    if ~isempty(products)
        for j=1:length(products)
            if j<length(products)
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(products(j)))) ' ' metaboliteList{products(j)} ' + ' ];
            else
                rxnForm=[rxnForm num2str(abs(stoichCoeffList(products(j)))) ' ' metaboliteList{products(j)}];
            end
        end
    end

    parsedKBaseFormulas{cnt,1}=translatedRxns{i,1};
    parsedKBaseFormulas{cnt,2}=rxnForm;
    parsedKBaseFormulas{cnt,3}=revFlag;
    cnt=cnt+1;
    % also catch if irreversible versions of this reaction exist in the VMH
    % database -> may consider translation
    if revFlag==1
        % also catch the reaction being written the other way around
        rxnForm=strsplit(rxnForm,' <=> ');
        parsedKBaseFormulas{cnt,1}=translatedRxns{i,1};
        parsedKBaseFormulas{cnt,2}=[rxnForm{1} ' -> ' rxnForm{2}];
        parsedKBaseFormulas{cnt,3}=revFlag;
        cnt=cnt+1;
        parsedKBaseFormulas{cnt,1}=translatedRxns{i,1};
        parsedKBaseFormulas{cnt,2}=[rxnForm{2} ' -> ' rxnForm{1}];
        parsedKBaseFormulas{cnt,3}=revFlag;
        cnt=cnt+1;
    end
end

% now find the overlap between the VMH database formulas and translated
% KBase reactions
[C,IA,IB]=intersect(parsedVMHFormulas(:,2),parsedKBaseFormulas(:,2));

% now get the overlap. If the reversibility flags agree, they are the same
% reaction, otherwise, they are similar (but translation should be
% considered).

for i=1:length(C)
    VMHform=database.reactions{find(strcmp(database.reactions(:,1),parsedVMHFormulas{IA(i),1})),3};
    KBaseform=translatedRxns{find(strcmp(translatedRxns(:,1),parsedKBaseFormulas{IB(i),1})),3};

    toAdd={parsedKBaseFormulas{IB(i),1},parsedVMHFormulas{IA(i),1},KBaseform,VMHform};
    % if they are the same reaction
    if parsedVMHFormulas{IA(i),3} == parsedKBaseFormulas{IB(i),3}
        sameReactions(size(sameReactions,1)+1,:)=toAdd;
    else
        similarReactions(size(similarReactions,1)+1,:)=toAdd;
    end
end

% remove similar reactions that are already in same reactions
[C,IA]=intersect(similarReactions(:,1),sameReactions(:,1),'stable');
similarReactions(IA(2:end),:)=[];

end

author: opencobra

repo: cobratoolbox


text:
function [sDataCT] = correctPatient30_HN(pathDATA)
% -------------------------------------------------------------------------
% function [sDataCT] = correctPatient30_HN(pathDATA)
% -------------------------------------------------------------------------
% DESCRIPTION:
% Reason of existence of the function: the RTstruct of the CT scan of
% patient 30 is incomplete (most probably it was previously incorrectly
% processed by MIM). The problem is solved by copying and processing the
% 'contour' structure of the PET scan.
% -------------------------------------------------------------------------
% INPUTS:
% - pathDATA: Full path to the 'DATA' folder of the HN workspace.
% -------------------------------------------------------------------------
% OUTPUTS:
% - sData: Corrected sData.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathDATA)
load('Patient30_CT'), sDataCT = sData;
load('Patient30_PET'), sDataPET = sData;

% ADJUSTING CT CONTOURS ACCORDING TO PET CONTOURS.
sDataCT{2}.scan.contour = sDataPET{2}.scan.contour;
nContour = length(sDataPET{2}.scan.contour);
downF = sDataCT{2}.scan.pixelW/sDataPET{2}.scan.pixelW;
offsetStart = round(abs(sDataPET{3}(1).ImagePositionPatient(2)-sDataCT{3}(1).ImagePositionPatient(2))*sDataCT{2}.scan.pixelW);
offsetEnd = round((size(sDataCT{2}.scan.volume,1)*sDataCT{2}.scan.pixelW - size(sDataPET{2}.scan.volume,1)*sDataPET{2}.scan.pixelW)/sDataCT{2}.scan.pixelW)-offsetStart;
for i = 1:nContour
    sDataCT{2}.scan.contour(i).boxBound(1:2,1) = round((sDataPET{2}.scan.contour(i).boxBound(1:2,1)/downF + offsetStart)); sDataCT{2}.scan.contour(i).boxBound(1:2,2) = round((sDataPET{2}.scan.contour(i).boxBound(1:2,2)/downF + offsetEnd));
    boxBound = sDataCT{2}.scan.contour(i).boxBound;
    nSlices = size(sDataCT{2}.scan.contour(i).boxMask,3);
    sDataCT{2}.scan.contour(i).boxMask = zeros(boxBound(1,2)-boxBound(1,1)+1,boxBound(2,2)-boxBound(2,1)+1,nSlices);
    for j = 1:nSlices
        sDataCT{2}.scan.contour(i).boxMask(:,:,j) = imresize(sDataPET{2}.scan.contour(i).boxMask(:,:,j),[boxBound(1,2)-boxBound(1,1)+1,boxBound(2,2)-boxBound(2,1)+1],'nearest');
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [cIX,gIX,numK] = AllCentroidRegression(hfig)
M_0 = getappdata(hfig,'M_0');
thres_reg = getappdata(hfig,'thres_reg');

Reg = FindCentroid(hfig);

Corr = corr(Reg',M_0');

[corr_max,IX] = max(Corr,[],1);
% cIX = find(corr_max>thres_reg)';

if thres_reg>=0
    cIX = (find(corr_max>thres_reg))';
elseif thres_reg<0
    cIX = (find(corr_max<thres_reg))';
end

gIX = IX(cIX)';
numK = length(unique(gIX));
end

author: xiuyechen

repo: FishExplorer


text:
function out = ij2xy(in, maxy)
% out = xy2ij(in, maxy)
%
% Converts shapes from "image" coordinates to cartesian coordinates.
% In cartesian coordinates, the origin is on the lower-left corner of
% the image and y moves vertically while x moves horizontally.
% Image coordinates follow a matrix notation where the origin is in
% the upper left corner and the first coordinate 'i' moves vertically
% (in the opposite direction of y), while the second 'j' coordinate
% moves horizontally in the same direction as x.
%
%                         PARAMETERS
%
% in A Nx2xS matrix containing S shapes and N landmarks,
%    in images coordinates: [i j]
% maxy Number of rows of the image matrix, or height of the displayed image,
%      in pixels.
%
%                          RETURNS
%
% out A Nx2xS matrix containing S shapes and N landmarks,
%     in cartesian coordinates: [x y]
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

	out(:,1,:) = in(:,2,:);
	out(:,2,:) = repmat(maxy, [size(in, 1) 1 size(in,3)]) - in(:,1,:);


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [] = plotExperiment()

    % plots the results of an experiment
    %
    % Parameters:
    % parameter: name of parameter on x-axis (e.g. param1 or event_metrics)
    % lines: name of lines in plot (e.g. param1 or event_metrics)

    % get algorithm and experiment name
    result_settings = ReadYaml('input/result.yaml');
    algorithm = result_settings.algorithm;
    configuration = result_settings.configuration;
    experiment_name = result_settings.experiment;

    % load experiment results
    path_to_experiment = strcat('results/summary/', algorithm, '/', configuration, '/', experiment_name);
    load(strcat(path_to_experiment, '/experiment_result.mat'));

    % get parameter (x-axis on plot) and restriction (some parameters are
    % restricted to specific values, e.g. only household 2 is used)
    parameter = result_settings.parameter;
    restriction = result_settings.restriction;

    % get names and values of the experiment parameters
    yaml_file_of_experiment = strcat('input/autogen/experiments/yaml/', algorithm, '_', configuration, '/', experiment_name, '.yaml');
    experiment = ReadYaml(yaml_file_of_experiment);
    experiment_parameters = {};
    field_names = fieldnames(experiment);
    for field_name = field_names'
        values = experiment.(field_name{1});
        if length(values) > 1
            experiment_parameters{end+1} = field_name{1};
        end
    end
    parameterID = find(strcmpi(parameter,experiment_parameters));
    load(strcat(path_to_experiment, '/usedParameterValues.mat'));
    parameter_values_in_cell = usedParameterValues.(parameter);
    parameter_values = cellfun(@str2double, parameter_values_in_cell);

    % get indexes of the lines in the plot
    lines = usedParameterValues.(result_settings.lines{1});
    res = zeros(length(parameter_values), length(lines));
    result = experiment_result.(result_settings.domain);
    if ~isempty(result_settings.lines{2})
        lines_idx = find(ismember(usedParameterValues.(result_settings.lines{1}), result_settings.lines{2}));
        lines = lines(lines_idx);
    else
        lines_idx = 1:length(lines);
    end

    % for each line, get appropriate index of result (experiment_result)
    % idx(param1, param2, ..., param_n, event/appliance_name, combine_res,
    % event/consumption_metric)
    for i = lines_idx
        idx(1:length(experiment_parameters)) = {':'};

        % set index of specified event/appliance name, event/consumption metric
        % and 'combine_results' (method for combining results of different training periods)
        if strcmpi(result_settings.domain, 'events')
            idx(length(experiment_parameters)+1) = {find(ismember(usedParameterValues.event_names, result_settings.event_name))};
            idx(length(experiment_parameters)+3) = {find(ismember(usedParameterValues.event_metrics, result_settings.metric))};
        else
            idx(length(experiment_parameters)+1) = {find(ismember(usedParameterValues.appliance_names, result_settings.appliance_name))};
            idx(length(experiment_parameters)+3) = {find(ismember(usedParameterValues.appliance_metrics, result_settings.metric))};
        end
        idx(length(experiment_parameters)+2) = {find(ismember(usedParameterValues.combine_results, result_settings.combine_res))};

        % some parameters are restricted to specific values, adapt index
        % accordingly
        for k=1:length(restriction)
            param_restricted = restriction{k};
            param_restricted_name = param_restricted{1};
            param_restricted_values = param_restricted{2};
            param_restricted_idx = find(ismember(usedParameterValues.(param_restricted_name), param_restricted_values));
            index_of_restricted_parameter = getIndexOfParameterName(param_restricted_name, experiment_parameters);
            idx(index_of_restricted_parameter) = {param_restricted_idx};
        end

        % set index of 'lines parameter'
        index_of_lines_parameter = getIndexOfParameterName(result_settings.lines{1}, experiment_parameters);
        idx(index_of_lines_parameter) = {i};

        % get result and store it in the 'res' matrix
        for j = 1:size(result,parameterID)
            idx(parameterID) = {j};
            r = result(idx{:});
            idx2 = r ~= 0 & ~isnan(r);
            res(j,i) = sum(r(idx2)) / nnz(idx2 ~=0);
        end
    end

    % plot the result (res matrix)
    set(0, 'DefaultAxesFontSize',12)
    set(0, 'DefaultTextFontName','Times New Roman');
    figure();
    [parameter_values_sorted, sortIdx] = sort(parameter_values);
     res_sorted = res(sortIdx, :);
    plot(parameter_values_sorted, res_sorted(:,lines_idx));

    % set legend of plot
    if iscellstr(lines)
        legend(lines);
    else
        legend(strcat(result_settings.lines{1}, {' '}, cellstr(num2str(lines'))'));
    end

    % set labels of axises in plot
    xlabel(parameter)
    if strcmpi(result_settings.lines{1}, 'event_metrics') || strcmpi(result_settings.lines{1}, 'appliance_metrics')
        ylabel('performance')
    else
        ylabel(result_settings.metric)
    end

    % set title of plot
    %title(['Performance of ''', algorithm ,''' algorithm with ''', configuration, ''' configuration:'],'FontWeight', 'Bold');

    grid on;

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [status,exitflag] = scipRetCode(code)
%Returns Status associated with SCIP Return Code

%Note v3.2.0 some numbers changed! Below as per scip.m
%   Return Status:
%       0 - Unknown
%       1 - User Interrupted
%       2 - Node Limit Reached
%       3 - Total Node Limit Reached
%       4 - Stall Node Limit Reached
%       5 - Time Limit Reached
%       6 - Memory Limit Reached
%       7 - Gap Limit Reached
%       8 - Solution Limit Reached
%       9 - Solution Improvement Limit Reached
%      10 - Restart Limit Reached
%      11 - Problem Solved to Optimality
%      12 - Problem is Infeasible
%      13 - Problem is Unbounded
%      14 - Problem is Either Infeasible or Unbounded

switch(code)
    case 1, status = 'User Exited'; exitflag = -5;
    case 2, status = 'Node Limit Reached'; exitflag = 0;
    case 3, status = 'Total Node Limit Reached'; exitflag = 0;
    case 4, status = 'Stall Node Limit Reached'; exitflag = 0;
    case 5, status = 'Time Limit Reached'; exitflag = 0;
    case 6, status = 'Memory Limit Reached'; exitflag = -4;
    case 7, status = 'Gap Limit Reached'; exitflag = 0;
    case 8, status = 'Solution Limit Reached'; exitflag = 0;
    case 9, status = 'Solution Improvement Limit Reached'; exitflag = 0;
    case 10, status = 'Restart Limit Reached'; exitflag = 0;
    case 11, status = 'Globally Optimal'; exitflag = 1;
    case 12, status = 'Infeasible'; exitflag = -1;
    case 13, status = 'Unbounded'; exitflag = -2;
    case 14, status = 'Unbounded or Infeasible'; exitflag = -2;
    otherwise, status = 'Unknown'; exitflag = -4;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function AddToolbarShortcut(varargin)
%ADDTOOLBARSHORTCUT is deprecated.  Use AddShortcut instead.

% $Author: rcotton $	$Date: 2010/08/20 10:21:23 $	$Revision: 1.7 $
% Copyright: Health and Safety Laboratory 2010

AddShortcut(varargin{:});

end


text:
%GETNEXTCLUSTERS Takes in a cluster graph and returns the indices
%   of the nodes between which the next message should be passed.
%
%   [i j] = GetNextClusters(P,Messages,oldMessages,m,useSmart)
%
%   INPUT
%     P - our cluster graph
%     Messages - the current values of all messages in P
%     oldMessages - the previous values of all messages in P. Thus,
%         oldMessages(i,j) contains the value that Messages(i,j) contained
%         immediately before it was updated to its current value
%     m - the index of the message we are passing (ie, m=0 indicates we have
%         passed 0 messages prior to this one. m=5 means we've passed 5 messages
%     useSmart - indicates whether we should use the Naive or Smart message
%         passing order
%
%
%   Output [i j]
%     i = the origin of the m+1th message
%     j = the destination of the m+1th message
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function [i j] = GetNextClusters(P,Messages,oldMessages,m,useSmart)

if(~exist('useSmart','var')||~useSmart)
  [i j] = NaiveGetNextClusters(P,m);
else
  [i j] = SmartGetNextClusters(P,Messages,oldMessages,m);
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [xdata,ydata] = jdfigure2data(figurename)
%%%derive the xdata and ydata from a figure(2x dimension)

%%Input:
%%%%----figurename       :      name of the figure      [required]

%%Output:
%%%%----xdata            :      data on x axis
%%%%----ydata            :      data on y axis
%%=====================================================================
    open(figurename);
    h_line = get(gca,'Children');
    xdata = get(h_line,'Xdata');
    ydata = get(h_line,'Ydata');
end

author: jindongwang

repo: activityrecognition


text:
 function y = jinc(x)
%function y = jinc(x)
% jinc(x) = J_1(pi x) / (2 x),
% where J_1 is Bessel function of the first kind of order 1.
% This is the 2D Fourier transform of a disk of diameter 1,
% so its DC value is the area of that disk which is pi/4.
% Equivalently it is the Hankel transform of rect = @(r) abs(r) < 1/2;
% Jeff Fessler, University of Michigan

x = abs(x); % kludge for bessel with negative arguments, perhaps not needed
y = pi/4 + 0 * x; % jinc(0) = pi/4
ig = x ~= 0;
y(ig) = besselj(1, pi * x(ig)) ./ (2 * x(ig));

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%%% This script alets the user enter the name of the images (base name, numbering scheme,...

% Checks that there are some images in the directory:

l_ras = dir('*ras');
s_ras = size(l_ras,1);
l_bmp = dir('*bmp');
s_bmp = size(l_bmp,1);
l_tif = dir('*tif');
s_tif = size(l_tif,1);
l_pgm = dir('*pgm');
s_pgm = size(l_pgm,1);
l_ppm = dir('*ppm');
s_ppm = size(l_ppm,1);
l_jpg = dir('*jpg');
s_jpg = size(l_jpg,1);

s_tot = s_ras + s_bmp + s_tif + s_pgm + s_jpg + s_ppm;

if s_tot < 1,
   fprintf(1,'No image in this directory in either ras, bmp, tif, pgm, ppm or jpg format. Change directory and try again.\n');
   break;
end;

% IF yes, display the directory content:

dir;

Nima_valid = 0;

while (Nima_valid==0),

   fprintf(1,'\n');
   calib_name = input('Basename camera calibration images (without number nor suffix): ','s');

   format_image = '0';

	while format_image == '0',

   	format_image =  input('Image format: ([]=''r''=''ras'', ''b''=''bmp'', ''t''=''tif'', ''p''=''pgm'', ''j''=''jpg'', ''m''=''ppm'') ','s');

		if isempty(format_image),
   		format_image = 'ras';
		end;

      if lower(format_image(1)) == 'm',
         format_image = 'ppm';
      else
         if lower(format_image(1)) == 'b',
            format_image = 'bmp';
         else
            if lower(format_image(1)) == 't',
               format_image = 'tif';
            else
               if lower(format_image(1)) == 'p',
                  format_image = 'pgm';
               else
                  if lower(format_image(1)) == 'j',
                     format_image = 'jpg';
                  else
                     if lower(format_image(1)) == 'r',
                        format_image = 'ras';
                     else
                        disp('Invalid image format');
                        format_image = '0'; % Ask for format once again
                     end;
                  end;
               end;
            end;
         end;
      end;
   end;

   check_directory;

end;

if (Nima_valid~=0),
   % Reading images:

   ima_read_calib_no_read; % may be launched from the toolbox itself

   fprintf(1,'\n');

   fprintf(1,'To display the thumbnail images of all the calibration images, you may run mosaic_no_read (may be slow)\n');

end;

repo: ekfmonoslam


text:
function pos2 = aberat (pos1, ve, tlight)

% this function corrects position vector for aberration of light.
% algorithm includes relativistic terms.  adapted from murray (1981)
% mon. notices royal ast. society 195, 639-648.

% input

%  pos1   = position vector of observed object, with reespect to
%           origin at observer (or the geocenter), components in au

%  ve     = velocity vector of observer (or the geocenter),
%           with respect to origin at solar system barycenter,
%           components in au/day (in)

%  tlight = light time from body to observer (or the geocenter) in days

% output

%  pos2 = position vector of observed object, with respect to
%         origin at observer (or the geocenter), corrected
%         for aberration, components in au

% ported from NOVAS 3.1

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

% light-time for one astronomical unit in seconds, from de-405

ausec = 499.0047838061d0;

% speed of light in au/day

c = 86400.0d0 / ausec;

tl = tlight;

p1mag = tl * c;

if (tl == 0.0d0)

    p1mag = sqrt(pos1(1)^2 + pos1(2)^2 + pos1(3)^2);

    tl = p1mag / c;
end

vemag = sqrt(ve(1)^2 + ve(2)^2 + ve(3)^2);

beta = vemag / c;

rdotv = pos1(1) * ve(1) + pos1(2) * ve(2) + pos1(3) * ve(3);

cosd = rdotv / (p1mag * vemag);

gammai = sqrt(1.0d0 - beta^2);

p = beta * cosd;

q = (1.0d0 + p / (1.0d0 + gammai)) * tl;

r = 1.0d0 + p;

for j = 1:1:3

    pos2(j) = (gammai * pos1(j) + q * ve(j)) / r;

end

author: goGPS-Project


text:
function [f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,callfun)
%ASSERT_SIGRESHAPE_PRE  Preprocess and handle dimension input.
%
%   Input parameters:
%      f            : signal, possibly ND-array
%      L            : L parameter
%      dim          : dim parameter
%      callfun      : Name of calling function
%   Output parameters:
%      f            : Input signal as matrix
%      L            : Verified L
%      Ls           : Length of signal along dimension to be processed
%      W            : Number of transforms to do.
%      dim          : Verified dim
%      permutedsize : pass to assert_sigreshape_post
%      order        : pass to assert_sigreshape_post

if ~isnumeric(f)
  error('%s: The input must be numeric.',callfun);
end;

D=ndims(f);

% Dummy assignment.
order=1;

if isempty(dim)
  dim=1;

  if sum(size(f)>1)==1
    % We have a vector, find the dimension where it lives.
    dim=find(size(f)>1);
  end;

else
  if (numel(dim)~=1 || ~isnumeric(dim))
    error('%s: dim must be a scalar.',callfun);
  end;
  if rem(dim,1)~=0
    error('%s: dim must be an integer.',callfun);
  end;
  if (dim<1) || (dim>D)
    error('%s: dim must be in the range from 1 to %d.',callfun,D);
  end;

end;

if (numel(L)>1 || ~isnumeric(L))
  error('%s: L must be a scalar.',callfun);
end;
if (~isempty(L) && rem(L,1)~=0)
  error('%s: L must be an integer.',callfun);
end;

if dim>1
  order=[dim, 1:dim-1,dim+1:D];

  % Put the desired dimension first.
  f=permute(f,order);

end;

Ls=size(f,1);

% If L is empty it is set to be the length of the transform.
if isempty(L)
  L=Ls;
end;

% Remember the exact size for later and modify it for the new length
permutedsize=size(f);
permutedsize(1)=L;

% Reshape f to a matrix.
if ~isempty(f)
  f=reshape(f,size(f,1),numel(f)/size(f,1));
end;
W=size(f,2);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat