text:
function D=Delta(Time_to_Maturity,Stock_Value,Stock_Volatility,Strike,Riskfree_Rate)

Time_to_Maturity=max(Time_to_Maturity,10^(-9));  % avoid division by zero message

d1 = 1/(Stock_Volatility*sqrt(Time_to_Maturity))*( log(Stock_Value/Strike) + ...
                (Riskfree_Rate* Stock_Volatility^2/2)*Time_to_Maturity  );

D =normcdf(d1,0,1);


text:
function mapVals = grayAmplitudeVals(V, ampType, varargin);
% Return event-related amplitudes for a gray view
%
%  mapVals = grayAmplitudeVals(V, ampType, [options]);
%
% ras, 09/30/2009
if notDefined('V'),				V = getSelectedGray;		end
if notDefined('ampType'),       ampType = 'z-score';		end

%% checks
% check that a GLM has been run
if ~isequal(viewGet(V, 'DTName'), 'GLMs')
	error('Need to be in the GLMs data type.')
end

% get scan num, # of non-null conditions
scan = viewGet(V, 'CurScan');
stim = er_concatParfiles(V);
N = length(stim.condNames) - 1;  % omit baseline condition

%% params
% default params
preserveCoords = 0;
saveAmps = 0;

% parse options (these will override the dialog values)
for ii = 1:2:length(varargin)
	val = varargin{ii+1};
	eval( sprintf('%s = val;', varargin{ii}) );
end

%% if the beta maps haven't been xformed from inplanes, do it now
testFile = fullfile( dataDir(V), sprintf('betas-predictor%i.mat', N) );
if ~exist(testFile, 'file')
	hI = initHiddenInplane('GLMs', 1);

	for i = 1:N     % loop across conditions
		mapPath = fullfile(dataDir(hI), 'RawMaps', ...
			sprintf('betas_predictor%i.mat', i));
		hI = loadParameterMap(hI, mapPath);
		V = ip2volParMap(hI, V, 0, 1, 'linear'); % trilinear interpolation

		mapPath = fullfile(dataDir(hI), 'RawMaps', ...
			sprintf('stdDev_predictor%i.mat', i));
		hI = loadParameterMap(hI, mapPath);
		V = ip2volParMap(hI, V, 0, 1, 'linear');
	end

	hI = loadParameterMap(hI, 'Residual Variance.mat');
	V = ip2volParMap(hI, V, 0, 1, 'linear');
end

if notDefined('whichConds'),
	whichConds=1:1:N;
else
	N=length(whichConds);
end

%% get the amplitude estimate for each voxel
switch lower(ampType)
	case 'subtracted-betas'
		for i = 1:N
			n=whichConds(i);
			mapName = sprintf('betas-predictor%i', n);
			V = loadParameterMap(V, mapName);
			mapVals(i,:) = V.map{scan};
		end

		% remove the "cocktail blank" or mean of all conditions for each voxel
		mapVals = mapVals - repmat( nanmean(mapVals), [N 1] );

	case {'raw-betas' 'rawbetas' 'raw betas' 'betas'}
		for i = 1:N
			n=whichConds(i);
			mapName = sprintf('betas-predictor%i', n);
			V = loadParameterMap(V, mapName);

			mapVals(i,:) = V.map{scan};
		end

	case {'z-score' 'z score' 'zscore'}
		% load the betas
		for i = 1:N
			n=whichConds(i);
			mapName = sprintf('betas-predictor%i', n);
			V = loadParameterMap(V, mapName);
			mapVals(i,:) = V.map{scan};
		end

		% load the residual variance map
		V = loadParameterMap(V, 'Residual Variance');
		sigmaVals = V.map{scan};

		% convert from res. var -> std. dev
		% (we need some info from the GLM model: degrees of freedom)
		modelFile = sprintf('Inplane/GLMs/Scan%i/glmSlice1.mat', scan);
		load(modelFile, 'dof');
		sigmaVals = sqrt( sigmaVals .^ 2 ./ dof );

		% remove the "cocktail blank" or DC component for each voxel
		mapVals = mapVals - repmat( nanmean(mapVals), [N 1] );

		% now, normalize by dividing by the estimated residual variance at
		% each voxel
		mapVals = mapVals ./ repmat( sigmaVals, [N 1] );

	otherwise
		error('Invalid amplitude type.')
end

%% return only values in the ROI? (Maybe we don't want to do this)
if isempty(V.ROIs)
    I = 1:size(mapVals, 2);
else
    I = roiIndices(V, V.ROIs(end).coords, preserveCoords);
end
mapVals = mapVals(:,I);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function mysavefig(fname, varargin)

% Usage: savefig(filename, fighdl, format, options)
%
% Saves a pdf, eps, png, jpeg, and/or tiff of the contents of the fighandle's (or current) figure.
% It saves an eps of the figure and the uses Ghostscript to convert to the other formats.
% The result is a cropped, clean picture. There are options for using rgb or cmyk colours,
% or grayscale. You can also choose the resolution.
%
% The advantage of savefig is that there is very little empty space around the figure in the
% resulting files, you can export to more than one format at once, and Ghostscript generates
% trouble-free files.
%
% If you find any errors, please let me know! (peder at axensten dot se)
%
% filename: File name without suffix.
%
% fighdl:  (default: gcf) Integer handle to figure.
%
% options: (default: '-r300', '-lossless', '-rgb') You can define your own
%          defaults in a global variable savefig_defaults, if you want to, i.e.
%          savefig_defaults= {'-r200','-gray'};.
% 'eps':   Output in Encapsulated Post Script (no preview yet).
% 'pdf':   Output in (Adobe) Portable Document Format.
% 'png':   Output in Portable Network Graphics.
% 'jpeg':  Output in Joint Photographic Experts Group format.
% 'tiff':  Output in Tagged Image File Format (no compression: huge files!).
% '-rgb':  Output in rgb colours.
% '-cmyk': Output in cmyk colours (not yet 'png' or 'jpeg' -- '-rgb' is used).
% '-gray': Output in grayscale (not yet 'eps' -- '-rgb' is used).
% '-lossless':  Use lossless compression, works on most formats. same as '-c0', below.
% '-c<float>':  Set compression for non-indexed bitmaps in PDFs -
%               0: lossless; 0.1: high quality; 0.5: medium; 1: high compression.
% '-r<integer>':  Set resolution.
% '-crop': Removes points and line segments outside the viewing area --
% permanently.
%          Only use this on figures where many points and/or line segments are outside
%          the area zoomed in to. This option will result in smaller vector files (has no
%          effect on pixel files).
% '-dbg':  Displays gs command line(s).
%
% EXAMPLE:
% savefig('nicefig', 'pdf', 'jpeg', '-cmyk', '-c0.1', '-r250');
% Saves the current figure to nicefig.pdf and nicefig.png, both in cmyk and at 250 dpi,
%          with high quality lossy compression.
%
% REQUIREMENT: Ghostscript. Version 8.57 works, probably older versions too, but '-dEPSCrop'
%          must be supported. I think version 7.32 or newer is ok.
%
% HISTORY:
% Version 1.0, 2006-04-20.
% Version 1.1, 2006-04-27:
% - No 'epstopdf' stuff anymore! Using '-dEPSCrop' option in gs instead!
% Version 1.2, 2006-05-02:
% - Added a '-dbg' option (see options, above).
% - Now looks for a global variable 'savefig_defaults' (see options, above).
% - More detailed Ghostscript options (user will not really notice).
% - Warns when there is no device for a file-type/color-model combination.
% Version 1.3, 2006-06-06:
% - Added a check to see if there actually is a figure handle.
% - Now works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too).
% - Now compatible with Ghostscript 8.54, released 2006-06-01.
% Version 1.4, 2006-07-20:
% - Added an option '-soft' that enables anti-aliasing on pixel graphics (on by default).
% - Added an option '-hard' that don't do anti-aliasing on pixel graphics.
% Version 1.5, 2006-07-27:
% - Fixed a bug when calling with a figure handle argument.
% Version 1.6, 2006-07-28:
% - Added a crop option, see above.
% Version 1.7, 2007-03-31:
% - Fixed bug: calling print with invalid renderer value '-none'.
% - Removed GhostScript argument '-dUseCIEColor' as it sometimes discoloured things.
% Version 1.8, 2008-01-03:
% - Added MacIntel: 'MACI'.
% - Added 64bit PC (I think, can't test it myself).
% - Added option '-nointerpolate' (use it to prevent blurring of pixelated).
% - Removed '-hard' and '-soft'. Use '-nointerpolate' for '-hard', default for '-soft'.
% - Fixed the gs 8.57 warning on UseCIEColor (it's now set).
% - Added '-gray' for pdf, but gs 8.56 or newer is needed.
% - Added '-gray' and '-cmyk' for eps, but you a fairly recent gs might be needed.
% Version 1.9, 2008-07-27:
% - Added lossless compression, see option '-lossless', above. Works on most formats.
% - Added lossy compression, see options '-c<float>...', above. Works on 'pdf'.
%   Thanks to Olly Woodford for idea and implementation!
% - Removed option '-nointerpolate' -- now savefig never interpolates.
% - Fixed a few small bugs and removed some mlint comments.
%
% TO DO: (Need Ghostscript support for these, so don't expect anything soon...)
% - svg output.
% - '-cmyk' for 'jpeg' and 'png'.
% - Preview in 'eps'.
% - Embedded vector fonts, not bitmap, in 'eps'.
%
% Copyright (C) Peder Axensten (peder at axensten dot se), 2006.

% KEYWORDS:     eps, pdf, jpg, jpeg, png, tiff, eps2pdf, epstopdf, ghostscript
%
% INSPIRATION:  eps2pdf (5782), eps2xxx (6858)
%
% REQUIREMENTS: Works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too).

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

  % parse my parameters
  doremove = false(1,length(varargin));
  fontname = [];
  fontsize = [];
  fontmode='scaled';
  width = -1;
  height = -1;
  i = 0;
  while true,
    i = i + 1;
    if i >= length(varargin), break; end
    if ~ischar(varargin{i}), continue; end
    switch lower(varargin{i}),
      case 'fontmode'
    fontmode = lower(varargin{i+1});
    if (~strcmp(fontmode,{'scaled','fixed'}))
      error('FontMode must be ''scaled'' or ''fixed''.');
    end
      case 'fontname'
        fontname = varargin{i+1};
        if ~any(strcmpi(fontname,listfonts)),
          error('FontName must be a member of listfonts');
        end
        doremove([i,i+1]) = true;
        i = i+1;
      case 'fontsize'
        fontsize = LocalToNum(varargin{i+1});
        if(~LocalIsPositiveScalar(fontsize))
          error('FontSize must be a numeric scalar > 0');
        end
        doremove([i,i+1]) = true;
        i = i+1;
      case 'width',
        width = LocalToNum(varargin{i+1});
        if(~LocalIsPositiveScalar(width))
          error('Width must be a numeric scalar > 0');
        end
        doremove([i,i+1]) = true;
        i = i+1;
      case 'height',
          height = LocalToNum(varargin{i+1});
          if(~LocalIsPositiveScalar(height))
            error('Height must be a numeric scalar > 0');
          end
        doremove([i,i+1]) = true;
        i = i+1;
    end
  end
  varargin(doremove) = [];
  realnargin = nargin - nnz(doremove);

  % Get handle.
	fighdl=		get(0, 'CurrentFigure'); % See gcf.						% Get figure handle.
	if((realnargin >= 2) && (numel(varargin{1}) == 1) && isnumeric(varargin{1}))
		fighdl=		varargin{1};
		varargin=	{varargin{2:end}};
	end
	if(isempty(fighdl)), error('There is no figure to save!?');	end
  H = fighdl;

  allText   = findall(H, 'type', 'text');
  allAxes   = findall(H, 'type', 'axes');
  allFont   = [allText; allAxes];

  old.objs = {};
  old.prop = {};
  old.values = {};

  % Process size parameters
  paperPos = get(H, 'PaperPosition');
  old = LocalPushOldData(old, H, 'PaperPosition', paperPos);
  figureUnits = get(H, 'Units');
  set(H, 'Units', get(H,'PaperUnits'));
  figurePos = get(H, 'Position');
  aspectRatio = figurePos(3)/figurePos(4);
  set(H, 'Units', figureUnits);
  if (width == -1) & (height == -1)
    width = paperPos(3);
    height = paperPos(4);
  elseif (width == -1)
    width = height * aspectRatio;
  elseif (height == -1)
    height = width / aspectRatio;
  end
  set(H, 'PaperPosition', [0 0 width height]);
  paperPosMode = get(H, 'PaperPositionMode');
  old = LocalPushOldData(old, H, 'PaperPositionMode', paperPosMode);
  set(H, 'PaperPositionMode', 'manual');
  oldPosition = get(H,'Position');
  newPosition = oldPosition;
  newPosition(3) = width;
  newPosition(4) = height;
  set(H,'Position',newPosition);
  old = LocalPushOldData(old,H,'Position',oldPosition);

  % Process font parameters
  if ~isempty(fontname)
    set(allFont,'fontname',fontname);
  end
  if (~isempty(fontmode))
    oldfonts = LocalGetAsCell(allFont,'FontSize');
    switch (fontmode)
      case 'fixed'
        oldfontunits = LocalGetAsCell(allFont,'FontUnits');
        old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits);
        set(allFont,'FontUnits','points');
        if (fontsize == -1)
          set(allFont,'FontSize',7);
        else
          set(allFont,'FontSize',fontsize);
        end
      case 'scaled'
        if (fontsize == -1)
          wscale = width/figurePos(3);
          hscale = height/figurePos(4);
          scale = min(wscale, hscale);
        else
          scale = 1;
        end
        newfonts = LocalScale(oldfonts,scale,5);
        set(allFont,{'FontSize'},newfonts);
      otherwise
        error('Invalid FontMode parameter');
    end
    % make sure we push the size after the units
    old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts);
  end

	op_dbg=		false;												% Default value.

	% Compression
	compr=			[' -dUseFlateCompression=true -dLZWEncodePages=true -dCompatibilityLevel=1.6' ...
					 ' -dAutoFilterColorImages=false -dAutoFilterGrayImages=false ' ...
					 ' -dColorImageFilter=%s -dGrayImageFilter=%s'];	% Compression.
	lossless=		sprintf (compr, '/FlateEncode', '/FlateEncode');
	lossy=			sprintf (compr, '/DCTEncode',   '/DCTEncode'  );
	lossy=			[lossy ' -c ".setpdfwrite << /ColorImageDict << /QFactor %g ' ...
					 '/Blend 1 /HSample [%s] /VSample [%s] >> >> setdistillerparams"'];

	% Create gs command.
	cmdEnd=			' -sDEVICE=%s -sOutputFile="%s"';					% Essential.
	epsCmd=			'';
	epsCmd=	[epsCmd ' -dSubsetFonts=true -dEmbedAllFonts=false -dNOPLATFONTS'];% Future support?
	epsCmd=	[epsCmd ' -dUseCIEColor=true -dColorConversionStrategy=/UseDeviceIndependentColor' ...
					' -dProcessColorModel=/%s'];						% Color conversion.
	pdfCmd=	[epsCmd ' -dAntiAliasColorImages=false' cmdEnd];
	epsCmd=	[epsCmd cmdEnd];

	% Get file name.
	if((realnargin < 1) || isempty(fname) || ~ischar(fname))				% Check file name.
		error('No file name specified.');
	end
	[pathstr, namestr] = fileparts(fname);
	if(isempty(pathstr)), fname= fullfile(cd, namestr);	end

	set(fighdl, 'Units', 'centimeters')									% Set paper stuff.
	sz=			get(fighdl, 'Position');
	sz(1:2)=	0;
	set(fighdl, 'PaperUnits', 'centimeters', 'PaperSize', sz(3:4), 'PaperPosition', sz);

	% Set up the various devices.
	% Those commented out are not yet supported by gs (nor by savefig).
	% pdf-cmyk works due to the Matlab '-cmyk' export being carried over from eps to pdf.
	device.eps.rgb=		sprintf(epsCmd, 'DeviceRGB',	'epswrite', [fname '.eps']);
	device.jpeg.rgb=	sprintf(cmdEnd,	'jpeg', 					[fname '.jpeg']);
%	device.jpeg.cmyk=	sprintf(cmdEnd,	'jpegcmyk', 				[fname '.jpeg']);
	device.jpeg.gray=	sprintf(cmdEnd,	'jpeggray',					[fname '.jpeg']);
	device.pdf.rgb=		sprintf(pdfCmd, 'DeviceRGB',	'pdfwrite', [fname '.pdf']);
	device.pdf.cmyk=	sprintf(pdfCmd, 'DeviceCMYK',	'pdfwrite', [fname '.pdf']);
	device.pdf.gray=	sprintf(pdfCmd, 'DeviceGray',	'pdfwrite', [fname '.pdf']);
	device.png.rgb=		sprintf(cmdEnd,	'png16m', 					[fname '.png']);
%	device.png.cmyk=	sprintf(cmdEnd,	'png???', 					[fname '.png']);
	device.png.gray=	sprintf(cmdEnd,	'pnggray', 					[fname '.png']);
	device.tiff.rgb=	sprintf(cmdEnd,	'tiff24nc',					[fname '.tiff']);
	device.tiff.cmyk=	sprintf(cmdEnd,	'tiff32nc', 				[fname '.tiff']);
	device.tiff.gray=	sprintf(cmdEnd,	'tiffgray', 				[fname '.tiff']);

	% Get options.
	global savefig_defaults;											% Add global defaults.
 	if( iscellstr(savefig_defaults)), varargin=	{savefig_defaults{:}, varargin{:}};
	elseif(ischar(savefig_defaults)), varargin=	{savefig_defaults, varargin{:}};
	end
	varargin=	{'-r300', '-lossless', '-rgb', varargin{:}};			% Add defaults.
	res=		'';
	types=		{};
	crop=		false;
	for n= 1:length(varargin)											% Read options.
		if(ischar(varargin{n}))
			switch(lower(varargin{n}))
			case {'eps','jpeg','pdf','png','tiff'},		types{end+1}=	lower(varargin{n});
			case '-rgb',				color=	'rgb';	deps= {'-depsc2'};
			case '-cmyk',				color=	'cmyk';	deps= {'-depsc2', '-cmyk'};
			case '-gray',				color=	'gray';	deps= {'-deps2'};
			case '-lossless',			comp=			0;
			case '-crop',				crop=			true;
			case '-dbg',				op_dbg=			true;
			otherwise
				if(regexp(varargin{n}, '^\-r[0-9]+$')), 	 res=  varargin{n};
				elseif(regexp(varargin{n}, '^\-c[0-9.]+$')), comp= str2double(varargin{n}(3:end));
				else	warning('pax:savefig:inputError', 'Unknown option in argument: ''%s''.', varargin{n});
				end
			end
		else
			warning('pax:savefig:inputError', 'Wrong type of argument: ''%s''.', class(varargin{n}));
		end
	end
	types=		unique(types);
	if(isempty(types)), error('No output format given.');	end

	if (comp == 0)														% Lossless compression
		gsCompr=		lossless;
	elseif (comp <= 0.1)												% High quality lossy
		gsCompr=		sprintf(lossy, comp, '1 1 1 1', '1 1 1 1');
	else																% Normal lossy
		gsCompr=		sprintf(lossy, comp, '2 1 1 2', '2 1 1 2');
	end

	% Generate the gs command.
	switch(computer)													% Get gs command.
		case {'MAC','MACI'},			gs= '/usr/local/bin/gs';
		case {'PCWIN','PCWIN64'},		gs= 'gswin32c.exe';
		otherwise,						gs= 'gs';
	end
	gs=		[gs		' -q -dNOPAUSE -dBATCH -dEPSCrop'];					% Essential.
	gs=		[gs		' -dUseFlateCompression=true'];						% Useful stuff.
	gs=		[gs		' -dAutoRotatePages=/None'];						% Probably good.
	gs=		[gs		' -dHaveTrueTypes'];								% Probably good.
	gs=		[gs		' ' res];											% Add resolution to cmd.

	if(crop && ismember(types, {'eps', 'pdf'}))							% Crop the figure.
		fighdl= do_crop(fighdl);
	end

	% Output eps from Matlab.
	renderer=	['-' lower(get(fighdl, 'Renderer'))];					% Use same as in figure.
	if(strcmpi(renderer, '-none')), renderer=	'-painters';	end		% We need a valid renderer.
	print(fighdl, deps{:}, '-noui', renderer, res, [fname '-temp']);	% Output the eps.

	% Convert to other formats.
	for n= 1:length(types)												% Output them.
		if(isfield(device.(types{n}), color))
			cmd=		device.(types{n}).(color);						% Colour model exists.
		else
			cmd=		device.(types{n}).rgb;							% Use alternative.
			if(~strcmp(types{n}, 'eps'))	% It works anyways for eps (VERY SHAKY!).
				warning('pax:savefig:deviceError', ...
						'No device for %s using %s. Using rgb instead.', types{n}, color);
			end
		end
		cmp=	lossless;
		if (strcmp(types{n}, 'pdf')),	cmp= gsCompr;		end			% Lossy compr only for pdf.
		if (strcmp(types{n}, 'eps')),	cmp= '';			end			% eps can't use lossless.
		cmd=	sprintf('%s %s %s -f "%s-temp.eps"', gs, cmd, cmp, fname);% Add up.
			status= system(cmd);										% Run Ghostscript.
			if (op_dbg || status), 		display (cmd),		end
	end
	delete([fname '-temp.eps']);										% Clean up.

  % Restore figure settings
  for n = 1:length(old.objs)
    set(old.objs{n}, old.prop{n}, old.values{n});
  end

end % ends savefig function

function fig= do_crop(fig)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	Remove line segments that are outside the view.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	haxes=	findobj(fig, 'Type', 'axes', '-and', 'Tag', '');
	for n=1:length(haxes)
		xl=		get(haxes(n), 'XLim');
		yl=		get(haxes(n), 'YLim');
		lines=	findobj(haxes(n), 'Type', 'line');
		for m=1:length(lines)
			x=				get(lines(m), 'XData');
			y=				get(lines(m), 'YData');

			inx=			(xl(1) <= x) & (x <= xl(2));	% Within the x borders.
			iny=			(yl(1) <= y) & (y <= yl(2));	% Within the y borders.
			keep=			inx & iny;						% Within the box.

			if(~strcmp(get(lines(m), 'LineStyle'), 'none'))
				crossx=		((x(1:end-1) < xl(1)) & (xl(1) < x(2:end))) ...	% Crossing border x1.
						|	((x(1:end-1) < xl(2)) & (xl(2) < x(2:end))) ...	% Crossing border x2.
						|	((x(1:end-1) > xl(1)) & (xl(1) > x(2:end))) ...	% Crossing border x1.
						|	((x(1:end-1) > xl(2)) & (xl(2) > x(2:end)));	% Crossing border x2.
				crossy=		((y(1:end-1) < yl(1)) & (yl(1) < y(2:end))) ...	% Crossing border y1.
						|	((y(1:end-1) < yl(2)) & (yl(2) < y(2:end))) ...	% Crossing border y2.
						|	((y(1:end-1) > yl(1)) & (yl(1) > y(2:end))) ...	% Crossing border y1.
						|	((y(1:end-1) > yl(2)) & (yl(2) > y(2:end)));	% Crossing border y2.
				crossp=	[(	(crossx & iny(1:end-1) & iny(2:end)) ...	% Crossing a x border within y limits.
						|	(crossy & inx(1:end-1) & inx(2:end)) ...	% Crossing a y border within x limits.
						|	crossx & crossy ...							% Crossing a x and a y border (corner).
						 ),	false ...
						];
				crossp(2:end)=	crossp(2:end) | crossp(1:end-1);		% Add line segment's secont end point.

				keep=			keep | crossp;
			end
			set(lines(m), 'XData', x(keep))
			set(lines(m), 'YData', y(keep))
		end
	end
end

function outData = LocalPushOldData(inData, objs, prop, values)
outData.objs = {inData.objs{:}, objs};
outData.prop = {inData.prop{:}, prop};
outData.values = {inData.values{:}, values};
end

function cellArray = LocalGetAsCell(fig,prop)
cellArray = get(fig,prop);
if (~isempty(cellArray)) & (~iscell(cellArray))
  cellArray = {cellArray};
end
end

function newArray = LocalScale(inArray, scale, minValue)
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  newArray{k} = max(minValue,scale*inArray{k}(1));
end
end

function bool = LocalIsPositiveScalar(value)
bool = isnumeric(value) & ...
       prod(size(value)) == 1 & ...
       value > 0;
end

function value = LocalToNum(value)
if ischar(value)
  value = str2num(value);
end
end

author: kristinbranson

repo: JAABA


text:
function val = source(x_g, param)

% Evaluate source terms (if any) at points x_g

% characteristic functions for source terms
    chi1_g = (x_g(:,1) < 0.2*param.L).*(x_g(:,2) < 0.5*param.W);
    chi3_g = (x_g(:,1) > 0.8*param.L).*(x_g(:,2) > 0.5*param.W);
%   chi2_g = ~chi1_g & ~chi3_g;

    val =   0*chi1_g +   0*chi3_g;
%

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function MS = createMesh3D(varargin)
% MeshStructure = createMesh3D(Nx, Ny, Nz, Width, Height, Depth)
% creates a uniform 3D mesh:
% Nx is the number of cells in x (horizontal) direction
% Ny is the number of cells in y (vertical) direction
% Nz is the number of cells in z (perpendicular) direction
% Lx is the domain length in x direction
% Ly is the domain length in y direction
% Lz is the domain length in z direction
%
% SYNOPSIS:
%   MeshStructure = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
%
% PARAMETERS:
%   Nx: number of cells in the x direction
%   Lx: domain length in x direction
%   Ny: number of cells in the y direction
%   Ly: domain length in y direction
%   Nz: number of cells in the z direction
%   Lz: domain length in z direction
%
% RETURNS:
%   MeshStructure.
%                 dimensions=3 (3D problem)
%                 numbering: shows the indexes of cellsn from left to right
%                 and top to bottom and back to front
%                 cellsize: x, y, and z elements of the cell size =[Lx/Nx,
%                 Ly/Ny, Lz/Nz]
%                 cellcenters.x: location of each cell in the x direction
%                 cellcenters.y: location of each cell in the y direction
%                 cellcenters.z: location of each cell in the z direction
%                 facecenters.x: location of interface between cells in the
%                 x direction
%                 facecenters.y: location of interface between cells in the
%                 y direction
%                 facecenters.z: location of interface between cells in the
%                 z direction
%                 numberofcells: [Nx, Ny, Nz]
%
%
% EXAMPLE:
%   Nx = 2;
%   Lx = 1.0;
%   Ny = 3;
%   Ly = 2.0;
%   Nz = 4;
%   Lz = 3.0;
%   m = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz);
%   [X, Y, Z] = ndgrid(m.cellcenters.x, m.cellcenters.y, m.cellcenters.z);
%   [Xf, Yf, Zf] = ndgrid(m.facecenters.x, m.facecenters.y, m.facecenters.z);
%   plot3(X(:), Y(:), Z(:), 'or')
%   hold on;
%   plot3(Xf(:), Yf(:), Zf(:), '+b')
%   legend('cell centers', 'cell corners');
%
% SEE ALSO:
%     createMesh1D, createMesh2D, createMeshCylindrical1D, ...
%     createMeshCylindrical2D, createCellVariable, createFaceVariable

% Written by Ali A. Eftekhari
% See the license file

if nargin==6
  % uniform 1D mesh
  Nx=varargin{1};
  Ny=varargin{2};
  Nz=varargin{3};
  Width=varargin{4};
  Height=varargin{5};
  Depth=varargin{6};
  % cell size is dx
  dx = Width/Nx;
  dy = Height/Ny;
  dz = Depth/Nz;
  G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
  CellSize.x= dx*ones(Nx+2,1);
  CellSize.y= dy*ones(Ny+2,1);
  CellSize.z= dz*ones(Nz+2,1);
  CellLocation.x= [1:Nx]'*dx-dx/2;
  CellLocation.y= [1:Ny]'*dy-dy/2;
  CellLocation.z= [1:Nz]'*dz-dz/2;
  FaceLocation.x= [0:Nx]'*dx;
  FaceLocation.y= [0:Ny]'*dy;
  FaceLocation.z= [0:Nz]'*dz;
elseif nargin==3
  % nonuniform 1D mesh
  facelocationX=varargin{1};
  facelocationY=varargin{2};
  facelocationZ=varargin{3};
  facelocationX=facelocationX(:);
  facelocationY=facelocationY(:);
  facelocationZ=facelocationZ(:);
  Nx = length(facelocationX)-1;
  Ny = length(facelocationY)-1;
  Nz = length(facelocationZ)-1;
  G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
  CellSize.x= [facelocationX(2)-facelocationX(1); ...
    facelocationX(2:end)-facelocationX(1:end-1); ...
    facelocationX(end)-facelocationX(end-1)];
  CellSize.y= [facelocationY(2)-facelocationY(1); ...
    facelocationY(2:end)-facelocationY(1:end-1); ...
    facelocationY(end)-facelocationY(end-1)];
  CellSize.z= [facelocationZ(2)-facelocationZ(1); ...
    facelocationZ(2:end)-facelocationZ(1:end-1); ...
    facelocationZ(end)-facelocationZ(end-1)];
  CellLocation.x= 0.5*(facelocationX(2:end)+facelocationX(1:end-1));
  CellLocation.y= 0.5*(facelocationY(2:end)+facelocationY(1:end-1));
  CellLocation.z= 0.5*(facelocationZ(2:end)+facelocationZ(1:end-1));
  FaceLocation.x= facelocationX;
  FaceLocation.y= facelocationY;
  FaceLocation.z= facelocationZ;
end
c=G([1,end], [1,end], [1, end]);
e1=G([1, end], [1, end], 2:Nz+1);
e2=G([1, end], 2:Ny+1, [1, end]);
e3=G(2:Nx+1, [1, end], [1, end]);
MS=MeshStructure(3, ...
  [Nx,Ny, Nz], ...
  CellSize, ...
  CellLocation, ...
  FaceLocation, ...
  c(:), ...
  [e1(:); e2(:); e3(:)]);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [ r, s, area ] = node_reference_q8 ( )

%*****************************************************************************80
%
%% NODE_REFERENCE_Q8 returns the basis nodes for an 8 node quadrilateral.
%
%  Discussion:
%
%    This element is known as the quadratic "serendipity" element.
%
%  Reference Element Q8:
%
%    |
%    1  4--7--3
%    |  |     |
%    |  |     |
%    S  8     6
%    |  |     |
%    |  |     |
%    0  1--5--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real R(8), S(8), the coordinates of the basis nodes.
%
%    Output, real AREA, the area of the element.
%
  r(1:8) = [ 0.0, 1.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0 ];
  s(1:8) = [ 0.0, 0.0, 1.0, 1.0, 0.0, 0.5, 1.0, 0.5 ];

  area = 1.0;

  return
end


text:
function volCoords = meshCursor2Volume(volView, msh)
%
%    volCoords = meshCursor2Volume(volView, [msh=viewGet(volView,'currentmesh')]);
%
%  Gets the current mrMesh cursor location and transforms it to mrVista
%  volume coordiantes.
%
% There are currently 3 ways to map the cursor to the VOLUME:
%
% 1. map the 3d cursor location directly to volume coords. This is the
% least accurate, since even smoothing the mesh will throw it off. If the
% mesh is not a real surface (eg. it's a cut plane or a DTI fiber), then
% this method will be used, since we have no other choice.
%
% 2. get the vertex index of the cursor and then find the 3d location of
% that vertex in the volume space. This avoids the smoothing problem in
% method 1, but is not consistent with all other mapping code that uses the
% vertexToGray map. This method will only be used if the vertex-to-gray
% transform field of the msh struct doesn't exist or is empty.
%
% 3. get the vertex index of the cursor and use the vertex-to-gray
% transform associated with the mesh to map these to layer-1 gray nodes.
% Then, we can just take the 3d coordinates of those layer 1 nodes. This is
% the most accurate method and is consistent with all our other mapping
% code. This method will be used by default, if it can.
%
% HISTORY:
%  2006.06.01 RFD: wrote it.
%  2006.10.24 RAS: doesn't crash, but warns, if cursor is outside volume
%  range.
%  2008.12.16 RFD & DY:
if(~exist('msh','var')||isempty(msh))
    msh = viewGet(volView,'currentmesh');
end

if(isfield(msh,'vertexGrayMap') && ~isempty(msh.vertexGrayMap))
    vertInd = mrmGet(msh,'cursorVertex');
    if vertInd < 1      % cursor not pointing to a volume vertex
        volCoords = [];
        warning('[%s]: Mesh Cursor is outside volume range.', mfilename); %#ok<*WNTAG>
        return
    end
    layerOneVolumeIndex = msh.vertexGrayMap(1,vertInd);

    % If the cursor position falls within the zone where functional data
    % was collected, translate this to a volume coordinate by grabbing the
    % position of the layer 1 node
    if(layerOneVolumeIndex>0 && layerOneVolumeIndex<=size(volView.coords,2) && ~isnan(layerOneVolumeIndex))
        volCoords = volView.coords(:,layerOneVolumeIndex)';
    % If not (and we have not computed layer 1 nodes for that position),
    % just get the cursor position at the gray/white boundary (which is
    % about 1mm off from where the layer 1 node position would be)
    else
        volCoords = round(msh.initVertices([2 1 3],vertInd)');
    end
else
    volCoords = round(mrmGet(msh,'cursor'));
end

return;


text:
function params = eventParamsEdit(params);
% GUI to edit/set event-related parameters.
%
% params = eventParamsEdit([params]);
%
% If omitted, params uses the default params in eventParamDefault.
% So, this code can be used to modify existing event-related parameters,
% or get a new set.
%
% ras 08/05.
if notDefined('params'), params = eventParamsDefault; end

% This is going to be complicated, and may be better done in GUIDE.

return


text:
function acc=svmTest(W,H,text_feature,data,train_ratio,C)
    load([data,'/group.txt']);
    numOfNode = size(group,1);
    if strcmp(data,'citeseer')
        numOfGroup = 6;
    elseif strcmp(data,'cora')
        numOfGroup = 7;
    else
        numOfGroup = 19;
    end

    group(:,1) = group(:,1) + 1;
    if strcmp(data,'wiki')==0
        group(:,2) = group(:,2) + 1;
    end
    group = sparse(group(:,1),group(:,2),ones(size(group(:,1))),numOfNode,numOfGroup);

    grouptmp=group;
    acc=0;
    features = [W' text_feature*H'];    % representation learned by TADW
    for i=1:size(features,2)
        if (norm(features(:,i))>0)
            features(:,i) = features(:,i)/norm(features(:,i));
        end
    end

    for i=1:10  % do the procedure for 10 times and take the average
        rp = randperm(numOfNode);
        testId = rp(1:floor(numOfNode*(1-train_ratio)));

        groupTest = group(testId,:);
        group(testId,:)=[];

        trainId = [1:numOfNode]';
        trainId(testId,:)=[];

        result=SocioDim(features, group, trainId, testId, C);
        [res b] = evaluate(result,groupTest);
        acc=acc+res.micro_F1;
        group=grouptmp;
    end
    acc=acc/10;

repo: TADW


text:
function mu = triangmf(z,a,b,c)
%TRIANGMF Triangular membership function.
%   MU = TRIANGMF(Z,A,B,C) computes a fuzzy membership function with a
%   triangular shape. Z is the input variable and can be a vector of any
%   length. A, B, and C are scalar parameters, such that B >= A and C >=
%   B, that define the triangular shape.
%
%       MU = 0                          Z < A
%       MU = (Z - A) ./ (B - A)         A <= Z < B
%       MU = 1 - (Z - B) ./ (C - B)     B <= Z < C
%       MU = 0                          C <= Z
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

mu = zeros(size(z));

low_side  = (a <= z) & (z < b);
high_side = (b <= z) & (z < c);

mu(low_side)  = (z(low_side) - a)./(b - a);
mu(high_side) = 1 - (z(high_side) - b)./(c - b);

author: dipum

repo: dipum-toolbox

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


text:
function [v,x,t,m,ze]=quadpeak(z)
%PEAK2DQUAD find quadratically-interpolated peak in a N-D array
%
%  Inputs:  Z(m,n,...)   is the input array (ignoring trailing singleton dimensions)
%                        Note: a row vector will have 2 dimensions
%
% Outputs:  V        is the peak value
%           X(:,1)  is the position of the peak (in fractional subscript values)
%           T        is -1, 0, +1 for maximum, saddle point or minimum
%           M        defines the fitted quadratic: z = [x y ... 1]*M*[x y ... 1]'
%           ZE       the estimated version of Z

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

persistent wz a
% first calculate the fixed matrix, a (can be stored if sz is constant)
sz=size(z);         % size of input array
psz=prod(sz);       % number of elements in input array
dz=numel(sz);       % number of input dimensions
mz=find(sz>1);      % non-singleton dimension indices
nm=numel(mz);       % number of non-singleton dimensions
vz=sz(mz);          % size of squeezed input array
dx=max(mz);         % number of output dimensions
if ~nm              % if the input array is a scalar
    error('Cannot find peak of a scalar');
end
nc=(nm+1)*(nm+2)/2;  % number of columns in A matrix
if min(vz)<3
    error('Need at least 3 points in each non-singleton dimension');
end
if isempty(wz) || numel(wz)~=numel(vz) || ~all(wz==vz)
    wz=vz;
    a=ones(psz,nc);
    ix=(0:psz-1)';
    for i=1:nm
        jx=floor(ix/sz(mz(i)));
        a(:,i+nc-nm-1)=1+ix-jx*sz(mz(i));
        ix=jx;
        a(:,(i^2-i+2)/2:i*(i+1)/2)=a(:,nc-nm:i+nc-nm-1).*repmat(a(:,i+nc-nm-1),1,i);
    end
    a=(a'*a)\a';        % converts to polynomial coeficients {x^2 xy y^2 x y 1]
end

% now find the peak

c=a*z(:);   % polynomial coefficients for this data
w=zeros(nm+1,nm+1);
i=1:(nm+1)*(nm+2)/2;
j=floor((sqrt(8*i-7)-1)/2);
w(i+j.*(2*nm+1-j)/2)=c;
w=(w+w.')/2; % make it symmetrical
mr=w(1:nm,1:nm);
we=w(1:nm,nm+1);
y=-(mr\we);
v=y'*we+w(nm+1,nm+1);  % value at peak

% insert singleton dimensions into outputs

x=zeros(dx,1);
x(mz)=y;
m=zeros(dx+1,dx+1);
mz(nm+1)=dx+1;
m(mz,mz)=w;
if nargout>2
    ev=eig(mr);
    t=all(ev>0)-all(ev<0);
end
if nargout>4
    ze=zeros(sz);
    scp=cumprod([1 sz(1:end-1)]);
    ivec=fix(repmat((0:psz-1)',1,dz)./repmat(scp,psz,1));
    xe=[1+ivec-repmat(sz,psz,1).*fix(ivec./repmat(sz,psz,1)) ones(psz,1)];
    ze=reshape(sum((xe*m).*xe,2),sz);
end

if ~nargout && nm<=2
    % plot the data
    desc={'Maximum','Saddle Point','Minimum'};
    if nargout<=2
        ev=eig(mr);
        t=all(ev>0)-all(ev<0);
    end
    if nm==1
        xax=linspace(1,psz,100);
        plot(xax,c(1)*xax.^2+c(2)*xax+c(3),'-r',1:psz,z(:),'ob',x,v,'^k');
        set(gca,'xlim',[0.9 psz+0.1]);
        ylabel('z');
        xlabel(sprintf('x%d',mz(1)));
        title(sprintf('\\Delta = %s: z(%.2g) = %.2g',desc{t+2},y(1),v));
    else
        ngr=17;
        xax=repmat(linspace(1,vz(1),ngr)',1,ngr);
        yax=repmat(linspace(1,vz(2),ngr),ngr,1);
        zq=(c(1)*xax+c(2)*yax+c(4)).*xax+(c(3)*yax+c(5)).*yax+c(6);
        hold off
        mesh(xax,yax,zq,'EdgeColor','r');
        hold on
        plot3(repmat((1:vz(1))',1,vz(2)),repmat(1:vz(2),vz(1),1),reshape(z,vz),'ob',y(1),y(2),v,'^k');
        hold off
        set(gca,'xlim',[0.9 vz(1)+0.1],'ylim',[0.9 vz(2)+0.1]);
        xlabel(sprintf('x%d',mz(1)));
        ylabel(sprintf('x%d',mz(2)));
        zlabel('z');
        title(sprintf('\\Delta = %s: z(%.2g,%.2g) = %.2g',desc{t+2},y(1),y(2),v));
    end
end

author: jtkim-kaist

repo: VAD

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


text:
function varargout = exit_modaldlg(varargin)
% EXIT_MODALDLG M-file for exit_modaldlg.fig
%      EXIT_MODALDLG by itself, creates a new EXIT_MODALDLG or raises the
%      existing singleton*.
%
%      H = EXIT_MODALDLG returns the handle to a new EXIT_MODALDLG or the handle to
%      the existing singleton*.
%
%      EXIT_MODALDLG('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in EXIT_MODALDLG.M with the given input arguments.
%
%      EXIT_MODALDLG('Property','Value',...) creates a new EXIT_MODALDLG or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before exit_modaldlg_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to exit_modaldlg_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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 14-Nov-2005 01:25:26

%* Author: Ranjith Unnikrishnan                                          *
%* Carnegie Mellon University, Vision and Mobile Robotics Laboratory     *
%* THE MATERIAL EMBODIED IN THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"     *
%* AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,      *
%* INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR      *
%* FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL CARNEGIE MELLON  *
%* UNIVERSITY BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,            *
%* SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY         *
%* KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,        *
%* LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF     *
%* THIRD PARTIES, WHETHER OR NOT CARNEGIE MELLON UNIVERSITY HAS BEEN     *
%* ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON        *
%* ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE     *
%* POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.                      *
%

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @exit_modaldlg_OpeningFcn, ...
                   'gui_OutputFcn',  @exit_modaldlg_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 exit_modaldlg is made visible.
function exit_modaldlg_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 exit_modaldlg (see VARARGIN)

% Choose default command line output for exit_modaldlg
handles.output = 'Yes';

% Update handles structure
guidata(hObject, handles);

% Insert custom Title and Text if specified by the user
% Hint: when choosing keywords, be sure they are not easily confused
% with existing figure properties.  See the output of set(figure) for
% a list of figure properties.
if(nargin > 3)
    for index = 1:2:(nargin-3),
        if nargin-3==index, break, end
        switch lower(varargin{index})
         case 'title'
          set(hObject, 'Name', varargin{index+1});
         case 'string'
          set(handles.text1, 'String', varargin{index+1});
        end
    end
end

% Determine the position of the dialog - centered on the callback figure
% if available, else, centered on the screen
FigPos=get(0,'DefaultFigurePosition');
OldUnits = get(hObject, 'Units');
set(hObject, 'Units', 'pixels');
OldPos = get(hObject,'Position');
FigWidth = OldPos(3);
FigHeight = OldPos(4);
if isempty(gcbf)
    ScreenUnits=get(0,'Units');
    set(0,'Units','pixels');
    ScreenSize=get(0,'ScreenSize');
    set(0,'Units',ScreenUnits);

    FigPos(1)=1/2*(ScreenSize(3)-FigWidth);
    FigPos(2)=2/3*(ScreenSize(4)-FigHeight);
else
    GCBFOldUnits = get(gcbf,'Units');
    set(gcbf,'Units','pixels');
    GCBFPos = get(gcbf,'Position');
    set(gcbf,'Units',GCBFOldUnits);
    FigPos(1:2) = [(GCBFPos(1) + GCBFPos(3) / 2) - FigWidth / 2, ...
                   (GCBFPos(2) + GCBFPos(4) / 2) - FigHeight / 2];
end
FigPos(3:4)=[FigWidth FigHeight];
set(hObject, 'Position', FigPos);
set(hObject, 'Units', OldUnits);

% Show a question icon from dialogicons.mat - variables questIconData
% and questIconMap
load dialogicons.mat

IconData=questIconData;
questIconMap(256,:) = get(handles.figure1, 'Color');
IconCMap=questIconMap;

Img=image(IconData, 'Parent', handles.axes1);
set(handles.figure1, 'Colormap', IconCMap);

set(handles.axes1, ...
    'Visible', 'off', ...
    'YDir'   , 'reverse'       , ...
    'XLim'   , get(Img,'XData'), ...
    'YLim'   , get(Img,'YData')  ...
    );

% Make the GUI modal
set(handles.figure1,'WindowStyle','modal')

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

% --- Outputs from this function are returned to the command line.
function varargout = exit_modaldlg_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;

% The figure can be deleted now
delete(handles.figure1);

% --- 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)

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

% --- 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)

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isequal(get(handles.figure1, 'waitstatus'), 'waiting')
    % The GUI is still in UIWAIT, us UIRESUME
    uiresume(handles.figure1);
else
    % The GUI is no longer waiting, just close it
    delete(handles.figure1);
end

% --- Executes on key press over figure1 with no controls selected.
function figure1_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Check for "enter" or "escape"
if isequal(get(hObject,'CurrentKey'),'escape')
    % User said no by hitting escape
    handles.output = 'No';

    % Update handles structure
    guidata(hObject, handles);

    uiresume(handles.figure1);
end

if isequal(get(hObject,'CurrentKey'),'return')
    uiresume(handles.figure1);
end

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over pushbutton1.
function pushbutton1_ButtonDownFcn(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)

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function [sizeArray, planC] = getUniformizedSize(planC)
%"getUniformizedSize"
%   Return the size of the uniformized dataset ([x y z]) and store it in
%   planC{indexS.scan}.uniformScanInfo.size for future reference.
%
%   If the uniformized data does not exist, returns and stores [0 0 0].
%
% JRA 11/14/03
%
% Usage: [sizeArray, planC] = getUniformizedSize(planC)

indexS = planC{end};

try
	uniformScanInfo = planC{indexS.scan}(1).uniformScanInfo;
	scanInfo = planC{indexS.scan}(1).scanInfo(1);

	%Find number of slices in whole uniformized set.
	nCTSlices = abs(uniformScanInfo.sliceNumSup - uniformScanInfo.sliceNumInf) + 1;
	nSupSlices = size(planC{indexS.scan}(1).scanArraySuperior, 3);
	if isempty(planC{indexS.scan}(1).scanArraySuperior), nSupSlices = 0;, end
	nInfSlices = size(planC{indexS.scan}(1).scanArrayInferior, 3);
	if isempty(planC{indexS.scan}(1).scanArrayInferior), nInfSlices = 0;, end
	zSize = nCTSlices + nSupSlices + nInfSlices;
	xSize = scanInfo(1).sizeOfDimension2;
	ySize = scanInfo(1).sizeOfDimension1;

    uniformScanInfo.size = [xSize ySize zSize];
    planC{indexS.scan}(1).uniformScanInfo = uniformScanInfo;
catch
    planC{indexS.scan}(1).uniformScanInfo.size = [0 0 0];
end

sizeArray = planC{indexS.scan}(1).uniformScanInfo.size;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [y, iter, lambda, status] = minimize_cubic_newton(H, g, sigma, options)
% Minimize a cubicly regularized quadratic via Newton root finding.
%
% [y, iter, lambda, status] = minimize_cubic_newton(H, g, sigma, options)
%
% Inputs: a symmetric matrix H of size n, a nonzero vector g of length n,
% a positive real sigma and an options structure. The code expects H to
% be tridiagonal, stored as a sparse matrix.
%
% The main output is a vector y of length n, which should minimize
%
%   f(y) = g'*y + (1/2)*y'*H*y + (1/3)*sigma*norm(y)^3.
%
% This is achieved by reducing the problem to a univariate root finding
% problem, where the unknown is a scalar lambda. This root is computed
% using a Newton method.
%
% Other outputs are iter (the number of Newton iterations completed),
% lambda (a real scalar, see below) and status. The latter is 0 if the
% target tolerance was reached, 1 if subsequent iterations induce no
% significant change, and -1 if the algorithm return because it reached
% the maximum number of iterations (see the options structure.)
% Non-negative status values are considered successes.
%
% The options structure must contain the following fields (between
% parentheses are some recommended values):
%   options.verbosity (3): to control how much information this function
%   prints to the command window. Anything below 6 silences the function.
%   options.maxiter_newton (100): maximum number of Newton iterations.
%   options.tol_newton (1e-16): tolerance on the root finding accuracy. See
%   in code for details.
%
% The code is based on Section 6 in
% Cartis, Gould and Toint, "Adaptive cubic regularisation methods for
% unconstrained optimization. Part I: motivation, convergence and numerical
% results", Mathematical Programming, 2011.
% https://link.springer.com/article/10.1007/s10107-009-0286-5
%
% Theorem 3.1 in the referenced paper states y is optimal if and only
% if it there exists a real lambda such that
%
% (H + lambda*I)y = -g,  lambda = sigma*||y||  and  H + lambda*I is psd,
%
% where psd means positive semidefinite. The other way around, if we
% find the corresponding scalar lambda, than we can recover y by
% solving a linear system (though this system might not have a unique
% solution in full generality.) Thus, the general strategy is to search
% for lambda rather than for y.
%
% See also: arc arc_lanczos

% This file is part of Manopt: www.manopt.org.
% Original authors: May 1, 2018,
%    Naman Agarwal, Brian Bullins, Nicolas Boumal and Coralia Cartis.
% Contributors:
% Change log:

    n = size(H, 1);

    % Pick an initial lambda that is cheap to compute and that surely makes
    % the shifted H positive definite.
    lambda = norm(H, 1) + 2;
    H_shifted = H + lambda*speye(n);

    % Compute the smallest eigenvalue of H, as we know the target lambda
    % must be at least as large as the negative of that, so that the
    % shifted H will be positive semidefinite.
    %
    % Since H ought to be sparse and tridiagonal, and since we only need
    % its smallest eigenvalue, this computation could be sped up
    % significantly. It does not appear to be a bottleneck, and eig is
    % simple and reliable, so we keep this for now.
    lambda_min = min(eig(H));
    left_barrier = max(0, -lambda_min);

    % Counter 'iter' holds the number of fully executed Newton iterations.
    iter = 0;
    while true

        if iter >= options.maxiter_newton
            % Iterations exceeded maximum number allowed.
            status = -1;
            return;
        end

        % If lambda has the correct value and the shifted H is positive
        % definite, then this y is a minimizer.
        y = -(H_shifted\g);
        ynorm = norm(y);

        % If the following quantity is zero, we have found a solution.
        phi = 1/ynorm - sigma/lambda;

        % Check if it is close enough to zero to stop.
        if abs(phi) <= options.tol_newton*ynorm
            status = 0;
            return;
        end
        psi = ynorm^2;

        % TODO: clarify this part of the code (see referenced paper).
        % The following is a Newton type of step on the equation
        % sigma/lambda = 1/sqrt(psi(lambda_prev)) ...
        %          - (lambda - lambda_prev)((psi'(lambda_prev))/2(psi)^1.5)
        delta_y = -(H_shifted\y);
        psi_prime = 2*(y'*delta_y);
        p0 = 2*sigma*(psi^(1.5));
        p1 = -2*psi - lambda*psi_prime;
        p2 = psi_prime;
        r = roots([p2 p1 p0]);
        del_lambda = max(r) - lambda;
        iter = iter + 1;

        % If the Newton step would bring us left of the left barrier, jump
        % instead to the midpoint between the left barrier and the current
        % lambda.
        if lambda + del_lambda <= left_barrier
            del_lambda = -.5*(lambda - left_barrier);
        end

        % If the step is so small that it numerically does not make a
        % difference when added to the current lambda, we stop.
        if abs(del_lambda) <= eps(lambda)
            status = 1;
            return;
        end

        % Update lambda
        H_shifted = H_shifted + del_lambda*speye(n);
        lambda = lambda + del_lambda;

        if options.verbosity >= 6
            fprintf(['lambda %.12e, ||y|| %.12e, lambda/sigma %.12e, ' ...
                     'phi %.12e\n\n'], lambda, ynorm, lambda / sigma, phi);
        end

    end

end

author: NicolasBoumal

repo: manopt


text:
function [hRect,api] = imrectLite( hParent, pos, lims, ar, varargin )
% A 'lite' version of imrect [OBSOLETE: use imrectRot].
%
% There are significant problems using imrect in a larger GUI, mostly
% due to the fact that imrect is such a heavy object.  This creates a
% simple, fast, and in some ways nicer version of imrect.
%
% The implementation uses two graphics object: a rectangle and a patch.
% The patch is just there to catch mouse clicks in case the rectangle does
% not have a face color.  Any valid property of the rectangle may be
% altered (such as 'FaceColor', 'Curvature', ... ), EXCEPT for its
% 'ButtonDownFcn', 'DeleteFcn' and 'Position' properties.  To change the
% 'Position' of the rectangle, use the output structure 'api' (described
% below).
%
% USAGE
%  [hRect,api] = imrectLite( [hParent],[pos],[lims],[ar],[varargin] )
%
% INPUTS
%  hParent    - [gca] object parent, typically an axes object (gca)
%  pos        - [] initial pos vector [x y w h] or [] or [x0 y0]
%  lims       - [] 2x2 matrix of [xMin xMax; yMin yMax]
%  ar         - [0] aspect ratio (if positive w/h is constrained to ar)
%  varargin   - [] parameters to RECTANGLE graphics object
%
% OUTPUTS
%  hRect      - handle to rectangle, use to alter color etc.
%  api        - interface allowing access to created object
%  .getPos()       - get position - returns 4 elt pos
%  .setPos(pos)    - set position (while respecting constraints)
%  .setAr(ar)      - set aspect ratio to fixed value
%  .setPosLock(b)  - if lock set (b==true), rectangle cannot change
%  .setSizLock(b)  - if lock set (b==true), rectangle cannot change size
%  .setPosChnCb(f) - whenever pos changes (even slightly), calls f(pos)
%  .setPosSetCb(f) - whenever pos finished changing, calls f(pos)
%  .uistack(...)   - calls 'uistack( [hRect hPatch], ... )', see uistack
%
% EXAMPLE
%  figure(1), imshow cameraman.tif
%  rectProp = {'EdgeColor','g','LineWidth',4,'Curvature',[.1 .1]};
%  lims = [get(gca,'xLim'); get(gca,'yLim')];
%  [h,api] = imrectLite( gca, [], lims, [],  rectProp{:} );
%  api.setPosChnCb( @(pos) disp(['        ' num2str(pos)]) );
%  api.setPosSetCb( @(pos) disp(['FINAL = ' num2str(pos)]) );
%
% See also IMRECT, RECTANGLE, PATCH
%
% Piotr's Image&Video Toolbox      Version 2.35
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% default arguments / globals
if( nargin<1 || isempty(hParent) ); hParent=gca; end
if( nargin<2 ); pos=[]; end
if( nargin<3 ); lims=[]; end
if( nargin<4 || isempty(ar) ); ar=0; end
[posChnCb,posSetCb]=deal([]);
posLock=false; sizLock=false; minSiz=1;

% get figure and axes handles
hAx = ancestor(hParent,'axes'); assert(~isempty(hAx));
hFig = ancestor(hAx,'figure'); assert(~isempty(hFig));
set( hFig, 'CurrentAxes', hAx );

% create rectangle as well as patch objects
if(isempty(pos)); vis='off'; else vis='on'; end
hPatch = patch('FaceColor','none','EdgeColor','none');
hRect = rectangle(varargin{:},'Visible',vis);
set([hRect hPatch],'ButtonDownFcn',{@btnDwn,0},'DeleteFcn',@deleteFcn);

% set / get position
if(length(pos)==4), setPos(pos); else
  btnDwn([],[],-1); waitfor(hFig,'WindowButtonUpFcn','');
end

% create api
api = struct('getPos',@getPos, 'setPos',@setPos, ...
  'setPosChnCb',@setPosChnCb, 'setPosSetCb',@setPosSetCb, ...
  'setPosLock',@setPosLock, 'setSizLock',@setSizLock, ...
  'setAr',@setAr, 'uistack',@uistack1 );

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

  function setPos( posNew, varargin )
    if(isempty(hRect) || isempty(hPatch)); return; end
    % constrain position to fall within lims
    pos = constrainPos( posNew, varargin{:} );
    % compute rectangle vertices
    xs = pos([1 1 1 1]); xs(2:3)=xs(2:3)+pos(3);
    ys = pos([2 2 2 2]); ys(3:4)=ys(3:4)+pos(4);
    vert = [xs' ys' [1;1;1;1]]; face=1:4;
    % draw objects
    set(hPatch,'Faces',face,'Vertices',vert);
    set(hRect,'Position',pos);
  end

  function pos = constrainPos( pos, anchor, pos0 )
    % constrain position to fall within lims
    if( ~isempty(lims) )
      posStr=pos(1:2);  posEnd=pos(1:2)+pos(3:4);
      posStr = min( max(posStr,lims(1:2)), lims(3:4)-minSiz );
      posEnd = max( min(posEnd,lims(3:4)), posStr+minSiz );
      pos = [posStr posEnd-posStr];
    end

    % now constrain for aspect ratio
    if( ar<=0 ); return; end;
    w=pos(3); h=pos(4);
    if( nargin<3 )
      w = min( w, h*ar );
      h = min( h, w/ar );
      pos(3:4) = [w h];
    else
      if( isempty(lims) )
        if( anchor(1)==0 ); w=h*ar; else h=w/ar; end;
      else
        mid = pos(1:2) + pos(3:4)/2;
        lr='MRL'; lr=lr(sign(anchor(1)-pos(1))+2);
        tb='MBT'; tb=tb(sign(anchor(2)-pos(2))+2);
        if( strcmp(tb,'M') )
          maxH=min(lims(4)-mid(2), mid(2)-lims(2))*2;
          switch lr
            case 'L'; maxW=min(lims(3)-pos(1),maxH*ar);
            case 'R'; maxW=min(pos(1)+pos(3)-lims(1),maxH*ar);
          end
          w=min(w,maxW); h=w/ar;
        else
          switch lr
            case 'L'; maxW=pos(1)+pos(3)-lims(1);
            case 'M'; maxW=min(lims(3)-mid(1), mid(1)-lims(1))*2;
            case 'R'; maxW=lims(3)-pos(1);
          end
          maxH=min(pos(2)+pos(4)-lims(2),maxW/ar);
        end
        h=min(h,maxH); w=h*ar;
      end
      if( strcmp(lr,'M') )
        pos0(1)=pos0(1)-(w-pos0(3))/2; pos0(3)=w;
      elseif( strcmp(tb,'M') )
        pos0(2)=pos0(2)-(h-pos0(4))/2; pos0(4)=h;
      end
      sgnSiz = [w h] .* (2*(pos(1:2)==anchor)-1);
      pos = shiftPos( anchor, pos0, sgnSiz );
    end
  end

  function pos = shiftPos( anchor, pos0, del )
    % shift pos by del unless anchor(i)==0
    pos = pos0;
    for i=1:2, if(anchor(i)==0), continue; end
      pos(i)=anchor(i); pos(i+2)=max(minSiz,abs(del(i)));
      if(del(i)<0); pos(i)=pos(i)-pos(i+2); end;
    end
  end

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

  function pnt = getCurPnt()
    pnt = get(hAx,'CurrentPoint');
    pnt = pnt([1,3]);
  end

  function d = axisUnitsPerCentimeter()
    % approximage absolute axes size
    units=get(hAx,'units');
    set(hAx,'units','centimeters');
    cms = get(hAx,'Position');
    cms = max(cms(3:4));
    set(hAx,'units',units);
    % axes size in axes units
    xLim=get(gca,'xLim'); yLim=get(gca,'yLim');
    sizInner = max([xLim(2)-xLim(1),yLim(2)-yLim(1)]);
    % units per centemeter
    d = sizInner/cms;
  end

  function [anchor,cursor] = getSide( pnt )
    % get pnt in coordintates of center of rectangle
    rs=pos(3:4)/2; pnt=pnt-pos(1:2)-rs;
    % side(i) -1=near lf/tp bnd; 0=in center; +1:near rt/bt bnd
    t = axisUnitsPerCentimeter()*.15; side=[0 0];
    for i=1:2
      ti = min(t,rs(i)/2);
      if( pnt(i)<-rs(i)+ti ), side(i)=-1; % near lf/tp boundary
      elseif( pnt(i)>rs(i)-ti), side(i)=1; % near rt/bt boundary
      end
    end
    % compute anchor in orig coord system(opposite of selected side)
    anchor=-side.*rs+pos(1:2)+rs; anchor(side==0)=0;
    % select cursor based on position
    cs={'bottom','botr','right','topr','top','topl','left','botl'};
    a=mod(round((atan2(side(1),side(2))/pi)*4),8)+1; cursor=cs{a};
  end

  function btnDwn( h, evnt, flag )
    if(isempty(hRect) || isempty(hPatch)); return; end;
    if(posLock); return; end; if(sizLock); flag=1; end;
    if( flag==-1 ) % create new rectangle
      if(isempty(pos)), anchor=ginput(1); else anchor=pos(1:2); end
      setPos( [anchor 1 1] );
      set(hRect,'Visible','on'); cursor='botr'; flag=0;
    elseif( flag==0 ) % resize (or possibly drag) rectangle
      pnt=getCurPnt(); [anchor,cursor] = getSide( pnt );
      if(all(anchor==0)), btnDwn(h,evnt,1); return; end
    elseif( flag==1 ) % move rectangle
      anchor=getCurPnt(); cursor='fleur';
    else
      assert(false);
    end
    set( hFig, 'Pointer', cursor );
    set( hFig, 'WindowButtonMotionFcn',{@drag,flag,anchor,pos} );
    set( hFig, 'WindowButtonUpFcn', @stopDrag );
  end

  function drag( h, evnt, flag, anchor, pos0 ) %#ok<INUSL>
    if(isempty(hRect) || isempty(hPatch)); return; end;
    pnt = getCurPnt(); del = pnt-anchor;
    if( flag==1 ) % shift rectangle by del
      setPos( [pos0(1:2)+del pos0(3:4)] );
    else % resize rectangle
      setPos(shiftPos(anchor,pos0,del),anchor,pos0);
    end
    drawnow
    if(~isempty(posChnCb)); posChnCb(pos); end;
  end

  function stopDrag( h, evnt ) %#ok<INUSD>
    set( hFig, 'WindowButtonMotionFcn','');
    set( hFig, 'WindowButtonUpFcn','');
    set( hFig, 'Pointer', 'arrow' );
    if(~isempty(posSetCb)); posSetCb(pos); end;
  end

  function deleteFcn( h, evnt ) %#ok<INUSD>
    [posChnCb,posSetCb]=deal([]); hs=[hPatch hRect];
    for i=1:length(hs), if(ishandle(hs(i))); delete(hs(i)); end; end
    [hRect,hPatch]=deal([]);
  end

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

  function pos0 = getPos(), pos0=pos; end

  function setAr( ar1 ), ar=ar1; setPos(pos); end

  function setPosChnCb( f ), posChnCb=f; end

  function setPosSetCb( f ), posSetCb=f; end

  function setPosLock( b ), posLock=b; end

  function setSizLock( b ), sizLock=b; end

  function uistack1( varargin )
    if(isempty(hRect) || isempty(hPatch)); return; end;
    uistack( [hRect hPatch], varargin{:} );
  end

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function ell_demo ( )

%*****************************************************************************80
%
%% ELL_DEMO demonstrates MESH2D on the L-shaped region.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELL_DEMO:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use MESH2D to mesh the L-shaped region.\n' );
  fprintf ( 1, '  Use MESH2D_TO_MEDIT to write the mesh to a file.\n' );

  warning off

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Set maximum element size HDATA.HMAX = 0.1\n' );

  v = [ 0.0, 0.0; 2.0, 0.0; 2.0, 1.0; 1.0, 1.0; 1.0, 2.0; 0.0, 2.0 ];

  hdata = [];
  hdata.hmax = 0.1;

  [ p, t ] = mesh2d ( v, [], hdata );

  [ nv, ~ ] = size ( v );
  [ np, ~ ] = size ( p );
  [ nt, ~ ] = size ( t );
  fprintf ( 1, '  %d boundary vertices input, %d nodes and %d triangles created\n', nv, np, nt );
  pause
%
%  Now call MESH2D_TO_MEDIT to write the data to a file.
%
  filename = 'ell.mesh';
  mesh2d_to_medit ( p, t, filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote mesh data to MEDIT mesh file "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELL_DEMO:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function h = hfun1 ( x, y )

%*****************************************************************************80
%
%% HFUN1 is a size-function for the L-shaped region.
%
%  Discussion:
%
%    The smallest size is at (1.0,1.0), and sizes increase as their distance
%    from that point increases.
%
  h = 0.01 + 0.1 * sqrt ( ( x - 1.0 ).^2  + ( y - 1.0 ).^2 );

  return
end


text:
function mtrFigureSmoothness(kSmooth,saveDir)

%Definitions
%numSamples = 2000;
bSaveImages = 1;

% Create temp directory currently if none provided
% if ieNotDefined(tempDir)
%     tempDir = 'temp';
% end
% mkdir(tempDir);
% genPlotsSphDist(kSmooth,numSamples,tempDir);
genPlotsSphFunc(kSmooth);
set(gcf,'Position',[504   751   354   195]);
if bSaveImages
    figFilename = fullfile(saveDir,['smoothPdf.png']);
    set(gcf,'PaperPositionMode','auto');
    print('-dpng', figFilename);
end
return;

function genPlotsSphFunc(kSmooth)
 % Going to color a sphere according to the desired distribution
figure;
%cmap = [autumn(255); [.25 .25 .25]];
colormap(gray);
[px py pz] = sphere(80);
testVec = [px(:) py(:) pz(:)]';
watValues = reshape(watsonPdf([0,0,1]',testVec,kSmooth),size(pz));
% Get hemisphere
pz(pz<0) = 0;
surf(px,py,pz,reshape(watValues,size(px)),'EdgeAlpha',0,'FaceAlpha',1);
colorbar;
caxis([0,max(watValues(:))]);
axis equal
axis vis3d
grid on;

return;

function genPlotsSphDist(kSmooth,numSamples,tempDir)

% Create script for sampling vectors according to desired distribution
tempScriptFile = fullfile(tempDir,'tempRun.sh');
fid = fopen(tempScriptFile,'wt');
fprintf(fid,'#!/bin/bash\n');
fprintf(fid,'sphsample -pdf watson -k %g -ns %g -outputfile tempDist.Bdouble\n',kSmooth,numSamples);
fclose(fid);
input(['Run ' tempScriptFile ' then press [return] when finished to continue.']);

% Cleanup and load the sample data
delete(tempScriptFile);
figure;
cmap = [autumn(255); [.25 .25 .25]];
sampleFile = fullfile(tempDir,'tempDist.Bdouble');
displaySphDist(sampleFile,numSamples);
%set(gca,'Postion',[0.13   0.11   0.775   0.815]);
title(['\DeltaS = ' num2str(round(180 / (pi*sqrt(kSmooth)))) ' \circ']);
delete(sampleFile);
return;

function displaySphDist(sampleFile,numSamples)

fid = fopen(sampleFile,'rb','b');
d = fread(fid,'double'); fclose(fid);
vecs = reshape(d,3,numSamples);
sub_vecs = vecs(:,vecs(3,:)>0);
plot3(sub_vecs(1,:), sub_vecs(2,:), sub_vecs(3,:), '.')
hold on
[px py pz] = sphere;
surf(px,py,pz,repmat(256,size(pz)),'EdgeAlpha',1,'FaceAlpha',0);
%alpha(0.3)
axis equal
axis vis3d
grid off;
% xlabel('x');
% ylabel('y');
% zlabel('z');

return;


text:
%% The aliquot parts of a number
%
% This demo file teaches about the aliquot parts of a number,
% and how to use the functions I've provided.
%
% Author: John D'Errico
%
% e-mail: woodchips@rochester.rr.com
%
% Release: 1.0
%
% Release date: 9/21/08

%% The aliquot parts of a number are its proper divisors
% For example, the number 12 has the list of prime factors
factor(12)

%%
% but its proper divisors are given by
aliquotparts(12)

%%
% See that N will not be listed as a proper divisor of itself, not
% for any value of N. In fact, the number 1 has no proper divisors.
aliquotparts(1)

%%
% Of course, prime numbers can have only 1 as a proper divisor.
aliquotparts(17)

%% The sum of the aliquot parts, or sum of proper divisors
aliquotsum([2 4 5 6 8 12 27 135 40320])

%% aliquotsum is efficient and vectorized
% In this test, the sum of the divisors for EVERY number from 1 to 100000
% is computed in short order.
N = (1:100000)';
tic,pdsum = aliquotsum(N);toc
%%
% Show the numbers up to 100000 with the largest relative sums of their proper divisors.
[pdsum,tags] = sort(pdsum./N,'descend');
disp('         N,      sigma(N)./N')
disp([N(tags(1:10)),pdsum(1:10)])

%% You can also count the number of divisors
% This is just the sum of the zero'th powers of the divisors. Logically,
% we expect all the prime numbers to have only one divisor.
N = (1:20)';
[N,aliquotsum(N,0)]

%%
% What number no larger than 100000 has the most proper divisors?
N = (1:100000)';
d = aliquotsum(N,0);
[maxdivisors,Nmax] = max(d)
%%
% As you might expect, that number has multiple small factors
factor(Nmax)

%% Or you can form the sum of the p'th powers
% Here we will compute the sum of the squares of the aliquot parts.
N = (1:20)';
[N,aliquotsum(N,2)]

%% Perfect numbers
% Perfect numbers have their aliquot sums equal to the
% number itself.
aliquotsum([6 28 496 8128])

%% Abundant numbers
% Just under 25% of all numbers are abundant
N = 100000;
100*sum(aliquotsum(1:N)>(1:N))/N

%% Hyper-abundant numbers
% Few (roughly 2%) of numbers are hyper-abundant.
% I've defined hyper-abundancy as an aliquot sum of the number N
% that is at least twice as large as N.
N = (1:100000)';
%%
% As a percentage...
100*sum(aliquotsum(N)>(2*N))./100000

%%
% Show the numbers up to 100000 with the largest relative sums of their proper divisors.
pdsum = aliquotsum(N);
[pdsum,tags] = sort(pdsum./N,'descend');
disp('         N,      sigma(N)./N')
disp([N(tags(1:10)),pdsum(1:10)])

%% Collosally-abundant numbers?
% How hyper-abundant can they get? Is there a maximum?
N = 1000000;
pdsum = aliquotsum(1:N);
[maxsum,whichN] = max(pdsum)
%%
% It does not appear that there is any maximum value
% to this ratio. Certainly it exceeds 3. Mathworld suggests
% <http://mathworld.wolfram.com/ColossallyAbundantNumber.html |(Collossally abundant)|>
% that the ratio can get fairly large. Add 1 to get the ratio they list.
maxsum/whichN

%% Odd abundant numbers are less common than the abundant numbers in general
% The smallest odd abundant number is 945.
N = (1:2:10000)';
pdsum = aliquotsum(N);
ind = find(pdsum>N);
N(ind)'

%% The odd abundant numbers are also "less" abundant
% The third column here is the extent of the overabundancy.
% Since 1.0 there would indicate a perfect number, these odd abundants
% are all clearly just barely so. I wonder, is there a limit to the
% extent of abundancy for the odd numbers? Are any hyper-abundant,
% as I've defined it above?
N = (1:2:1000000)';
pdsum = aliquotsum(N);
maxOddAbundancy = max(pdsum./N)

%% Why are the odd abundants so rare?
% We can get some idea as to the reason why odd abundant numers
% are so rare, by looking at the factors of a highly abundant even
% number.
%
% Try 840 for example. Clearly it has multiple powers of 2 in its
% factorization.
factor(840)
%%
% 840 is also quite abundant.
aliquotsum(840)
%%
% If we look at the factors of the most highly abundant numbers,
% you will generally find many spare powers of 2. But how about
% the odd abunbdants? There are no powers of 2 allowed in an odd number,
% so pick some of the odd abundant numbers and look at their factors.
% As expected, we find some extra powers of 3, along with some other
% moderately small odd primes.
factor(945)
factor(9765)
%%
% Pick a rather large odd number that will have very many factors.
% Unfortunately, this is close to as far as aliquotsum will go, due to
% the limits of MATLAB's operating precision. I'll admit the ratio
% is starting to approach 2.
N = 3*3*3*5*5*7*7*11*13*17*19
aliquotsum(N)/N

%% Deficient numbers
% All primes are deficient, as are all pure powers of prime numbers.
% As well, all the divisors of any perfect number are known to be deficient.
aliquotsum(primes(50))

%% Maximally deficient numbers
% Naturally, the prime numbers represent the most deficient numbers
% possible, since each prime has only one divisor less than itself,
% and that is the number 1. If we ignore the prime numbers, which
% numbers are the next most deficient? Do you see a pattern?
N = setdiff(1:50,primes(50))';
[N,aliquotsum(N)]

%% Amicable numbers
% A pair of numbers is amicable if their aliquot sums are
% equal to the other member of the pair.
aliquotsum([220 284])

%% Sociable numbers form cliques, or amicable cycles
% Perfect numbers are self-amicable, with a cycle length of 1.
% Amicable pairs have a cycle length of 2.
% This cycle has length 5 before it returns to the start.
aliquotsum([12496 14288 15472 14536 14264])

%% Amicable cycles
% Perfect numbers are self-amicable. The amicablecycles function:
%
%  cycles = amicablecycles(N,L)
%
% will find all cycles that start with a number
% as large as N. The length of the cycle is L.
amicablecycles(20000,1)
%%
% Amicable pairs have a cycle length of 2.
amicablecycles(20000,2)

%%
% There are no known sociable cliques of length 3 or 4
amicablecycles(20000,3)
amicablecycles(20000,4)

%%
% The smallest cycle of length 5 starts at 12496
amicablecycles(20000,5)

%% Odd amicables
% You can force amicablecycles to look only at certain starting points.
% For example, to find the odd amicable cycles only
amicablecycles(1:2:200000,2)

%% Some interesting things to try
% What fun things can you find to do with these numbers?
% Can you find any odd perfect numbers? Perhaps a quasi-perfect
% number?
%
% Some simple problems for the student, and a few that may not be so
% simple.
%
% 1. Which numbers have many distinct divisors?
%
% 2. What is the smallest number to have exactly 31 divisors?
%
% 3. Does the product of the first k primes necessarily have the largest number of distinct divisors?
%
% 4. Can you construct a number less than 2*3*5*7*11 that has more divisors than this number?
%
% 5. For which values of non-prime numbers N is the lower bound sqrt(N)+1 realized for the aliquot sum?
%
% 6. Can you prove that sqrt(N)+1 forms a lower bound for the aliquot sum for the non-prime numbers N?
%
% 7. Is there a simple relationship for the upper bound of the aliquot sum
% as a function of N? Consider N*log(N) as a start.
%
% 8. There are many abundant numbers. As was shown above, roughly 25% of
% all numbers are abundant. Can you generate the list of all abundant
% numbers up to some value? (Start with the aliquotsum function.) What is the
% smallest abundant number?
%
% 9. The odd abundant numbers were shown to be somewhat rare, and always
% seemed to have multiple powers of 3 in their representations. Are there
% any abundants that do not have either 2 or 3 as a factor? Can this
% happen? Note that such an investigation might involve numbers that are
% too large for the tools I've provided to work with properly.
N = vpi(5)*5*5*5*5*5*7*7*7*11*13
aliquotsum(N)
%%
% 10. Can you represent numbers as the sum of two abundant numbers? Which
% numbers are not so representable? Is there a smallest number such that
% all those above it are the sum of two abundant numbers?
%
% 11. An interesting generalization of a perfect number is what I'll call the p-perfect
% number. Thus, if a perfect number is an integer such that the sum of its
% proper divisors adds up to the number itself, a p-perfect number N might
% have the sum of the p'th powers of the divisors adding up to N^p. Do any
% such numbers exist? Can you show this cannot happen? Or are there any
% numbers with the squares of the divisors that adds up to the original
% number?
%
% 12. How about p-amicability? Can you think of a way to generalize this
% concept?

author: Sable

repo: mcbench-benchmarks


text:
function [BG, FG, L_hat, S_hat, T_hat, t_hat, ...
    P_track_full, T_calc]= ...
    MEDRoP(M1, P_init, mu, ev_thresh, alpha, K)
%%%This is the MEDRoP algorithm. This is the main function for the real video Background-
%%%Foreground separation problem.

%This folder contains the code accompanying pre-print.

%%%                          Inputs                         %%%
%%%     M - measurement matrix                              %%%
%%%     ev_thres - threshold for subspace change detection  %%%
%%%     P_init - an initial estimate of the subspace        %%%
%%%     t_train - the dimension of the training data        %%%

%%%                       Algorithm parameters              %%%
%%%     alpha - frame length                                %%%
%%%     mu - column-averages of data                        %%%
%%%     K - number of projection PCA steps                  %%%
%%%     omega - threshold for non-zero value in S           %%%

%%%                          Outputs                        %%%
%%%     BG - Estimated Background                           %%%
%%%     FG - Estimated Foreground                           %%%
%%%     L_hat - estimate of the low rank matrix             %%%
%%%     P_hat - estimate of the subspace in which L lies    %%%
%%%     S_hat - estimate of the sparse signal               %%%
%%%     t_hat - estimate of subspace change times           %%%

%% Initializations

P_hat_old = P_init;
[~, r] = size(P_init);
P_hat = P_hat_old;

[n, t_max] = size(M1);
T_hat = zeros(n, t_max);
S_hat = zeros(n, t_max);
L_hat = zeros(n, t_max);

FG = zeros(n, t_max);
BG = zeros(n, t_max);

t_hat = [];
M = M1 - repmat(mu, 1, t_max);
k = 0;
cnt = 1;
ph = 0;     %ph - 0 => detect, 1 => ppca
opts.delta = 0.4;

phi_t = speye(n) - P_hat * P_hat';
%% Main Algorithm Loop
for ii = 2 : t_max
        %% Estimate support
        Atf.times = @(x) x - (P_hat * (P_hat' * x));
        Atf.trans = @(y) y - (P_hat * (P_hat' * y));
        phi.times = @(x) x - (P_hat_old * (P_hat_old' * x));
        y_t = Atf.times(M(:, ii));
        opts.tol   = 1e-3;
        opts.print = 0;

        opts.delta = norm(Atf.times(L_hat(:, ii - 1)));

        x_t_hat_cs = yall1(Atf, y_t, opts);
        omega = sqrt(M(:, ii)' * M(:, ii) / n);

        t_hat_temp = find(abs(x_t_hat_cs) > omega);
        T_hat(t_hat_temp, ii) = 255;

        LS.times = @(x) phi(:, t_hat_temp) * x;
        LS.trans = @(y) phi(:, t_hat_temp)' * x;

        %% Estimate signal components

        S_hat(t_hat_temp, ii) = cgls(phi_t(:, t_hat_temp), y_t, 0, 1e-3);
        L_hat(:, ii) = M(:, ii) - S_hat(:, ii);
        FG(t_hat_temp, ii) = M1(t_hat_temp, ii);
        BG(:, ii) = L_hat(:, ii) + mu;
        %% Subspace update
        if(~mod(ii - 1 , alpha))
            MM = (1 / sqrt(alpha)) * phi.times(L_hat(:, ii - alpha + 1 : ii));

            if(~ph)     %%detect phase
                aa = svds(MM, 1);
                if(aa >= sqrt(ev_thresh))
                    ph = 1;
                    t_hat = [t_hat, ii];
                    k = 0;
                end
            else        %%ppca phase
                P_hat = simpleEVD((L_hat(:, ii - alpha + 1 : ii)), r);
                phi_t = speye(n) - P_hat * P_hat';
                k = k + 1;

                if(k==K + 1)
                    P_hat_old = P_hat;
                    k = 1;
                    ph = 0;
                    phi_t = speye(n) - P_hat * P_hat';
                end
            end
        end

    %% Return subspace estimates
    if((ii == 0) || (~(mod(ii - 1, alpha))))
        P_track_full{cnt} = P_hat;
        %P_track_new{cnt} = P_hat_new;
        T_calc(cnt) = ii;
        cnt = cnt + 1;
    end
end
%L_hat = L_hat + repmat(mu, n, 1);
end

author: andrewssobral

repo: lrslibrary


text:
function [seq_ids, frame_ids] = load_dataset_indexes(filename)

fid = fopen(filename, 'r');
C = textscan(fid, '%s');
indexes = C{1};
fclose(fid);

num = numel(indexes);
seq_ids = zeros(num, 1);
frame_ids = zeros(num, 1);

for i = 1:num
    str = indexes{i};
    pos = strfind(str, '/');
    seq_ids(i) = str2double(str(1:pos-1));
    frame_ids(i) = str2double(str(pos+1:end));
end

author: yuxng


text:
function [txt_han,txt_lft] = coefdisp(txt_lft,txt_btm,coeffs,paren,data)
% COEFDISP Display of coefficients.

% Author: Craig Borghesani
% Date: 11/15/94
% Revised:
% Copyright (c) 1999, Prentice-Hall

len_coeff = length(coeffs);
txtct = 1;
ers = 'norm';
fontname = 'times';

for loc_coeff = 1:len_coeff,

 cur_coeff = coeffs(loc_coeff);
 power = len_coeff - loc_coeff;

 if cur_coeff~=0,
  coeff_str = num2str(abs(cur_coeff));
  if power > 0 & abs(cur_coeff) == 1,
   coeff_str = [];
  end

  if cur_coeff > 0 & loc_coeff > 1,
   coeff_str = ['+',coeff_str];
  elseif cur_coeff < 0,
   coeff_str = ['-',coeff_str];
  end

  if power > 0,
   cur_str = [coeff_str,'s'];
  else
   cur_str = coeff_str;
  end

  if paren,
   if loc_coeff == 1,
    cur_str = ['(',cur_str];
   elseif loc_coeff == len_coeff,
    cur_str = [cur_str,')'];
   end
  end

  txt_han(txtct) = text(txt_lft,txt_btm,cur_str,'erase',ers,...
              'horizontalalignment','left','verticalalignment','bottom',...
              'fontname',fontname);
  set(txt_han(txtct),'userdata',data);
  txt_ext = get(txt_han(txtct),'extent');
  txt_lft = txt_lft + txt_ext(3);
  txtct=txtct+1;

  if power > 1,
   sup_btm = txt_btm + txt_ext(4)/3;
   pow_str = num2str(power);
   txt_han(txtct) = text(txt_lft,sup_btm,pow_str,'erase',ers,...
      'horizontalalignment','left','verticalalignment','bottom',...
      'fontsize',8,'fontname',fontname);
   set(txt_han(txtct),'userdata',data);
   txt_ext = get(txt_han(txtct),'extent');
   txt_lft = txt_lft + txt_ext(3);
   txtct=txtct+1;
  end
 end
end


text:
function M = symfixedrankYYfactory(n, k)
% Manifold of n-by-n symmetric positive semidefinite matrices of rank k.
%
% function M = symfixedrankYYfactory(n, k)
%
% The geometry is based on the paper,
% M. Journee, P.-A. Absil, F. Bach and R. Sepulchre,
% "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices",
% SIAM Journal on Optimization, 2010.
%
% Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
%
% A point X on the manifold is parameterized as YY^T where Y is a matrix of
% size nxk. The matrix Y (nxk) is a full column-rank matrix. Hence, we deal
% directly with Y.
%
% Notice that this manifold is not complete: if optimization leads Y to be
% rank-deficient, the geometry will break down. Hence, this geometry should
% only be used if it is expected that the points of interest will have rank
% exactly k. Reduce k if that is not the case.
%
% An alternative, complete, geometry for positive semidefinite matrices of
% rank k is described in Bonnabel and Sepulchre 2009, "Riemannian Metric
% and Geometric Mean for Positive Semidefinite Matrices of Fixed Rank",
% SIAM Journal on Matrix Analysis and Applications.

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:
%  July 10, 2013 (NB)
%       Added vec, mat, tangent, tangent2ambient ;
%       Correction for the dimension of the manifold.

M.name = @() sprintf('YY'' quotient manifold of %dx%d PSD matrices of rank %d', n, k);

M.dim = @() k*n - k*(k-1)/2;

% Euclidean metric on the total space
M.inner = @(Y, eta, zeta) trace(eta'*zeta);

M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));

M.dist = @(Y, Z) error('symfixedrankYYfactory.dist not implemented yet.');

M.typicaldist = @() 10*k;

M.proj = @projection;
    function etaproj = projection(Y, eta)
        % Projection onto the horizontal space
        YtY = Y'*Y;
        SS = YtY;
        AS = Y'*eta - eta'*Y;
        Omega = lyap(SS, -AS);
        etaproj = eta - Y*Omega;
    end

M.tangent = M.proj;
M.tangent2ambient = @(Y, eta) eta;

M.retr = @retraction;
    function Ynew = retraction(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Ynew = Y + t*eta;
    end

M.egrad2rgrad = @(Y, eta) eta;
M.ehess2rhess = @(Y, egrad, ehess, U) M.proj(Y, ehess);

M.exp = @exponential;
    function Ynew = exponential(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end

        Ynew = retraction(Y, eta, t);
        warning('manopt:symfixedrankYYfactory:exp', ...
            ['Exponential for symmetric, fixed-rank ' ...
            'manifold not implemented yet. Used retraction instead.']);
    end

% Notice that the hash of two equivalent points will be different...
M.hash = @(Y) ['z' hashmd5(Y(:))];

M.rand = @random;

    function Y = random()
        Y = randn(n, k);
    end

M.randvec = @randomvec;
    function eta = randomvec(Y)
        eta = randn(n, k);
        eta = projection(Y, eta);
        nrm = M.norm(Y, eta);
        eta = eta / nrm;
    end

M.lincomb = @lincomb;

M.zerovec = @(Y) zeros(n, k);

M.transp = @(Y1, Y2, d) projection(Y2, d);

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

end

% Linear conbination of tangent vectors
function d = lincomb(Y, a1, d1, a2, d2) %#ok<INUSL>

if nargin == 3
    d  = a1*d1;
elseif nargin == 5
    d = a1*d1 + a2*d2;
else
    error('Bad use of symfixedrankYYfactory.lincomb.');
end

end


text:
function [as,bs,vars] = sumARMA(varargin)

%SUMARMA Sum of ARMA-processes
%   [ar_s,ma_s,var_s] = sumARMA(ar1,ma1,var1,ar2,ma2,var2) yields
%   the sum of 2 ARMA-processes. The list of arguments can be
%   extended with more ARMA-processes.

%S. de Waele, 2001.

%Technical note:
%If this method does not work, the following procedure can be followed:
%- Calculate the covariance functions and add them;
%- Transform this covariance function into an AR process.

if rem(nargin,3),
    error('Incorrect number of inputs.')
end
np = nargin/3;

%Initialisation
a  = cell(np,1);
b  = cell(np,1);
var= cell(np,1);
gain = cell(np,1);
nomc = cell(np,1);
for i = 1:np,
    a{i}    = varargin{1+3*(i-1)};
    b{i}    = varargin{2+3*(i-1)};
    var{i}  = varargin{3+3*(i-1)};
    gain{i} = pgain(b{i},a{i});
    nomc{i} = 1;
end
den = 1;
vars = 0;
nom = 0;

%Calculation of the contributions of the processes to the nominator
for i = 1:np,
    for j = 1:np,
        if j~=i,
            nomc{j} = conv(nomc{j},a{i});
        else
            nomc{j} = conv(nomc{j},b{j});
        end %if j~i,
    end %for j = 1:np,
    vars = vars+var{i};
    den = conv(den,a{i});
end %for i = 1:np,
m = length(nomc{1});
for i = 2:np, m = max( [m length(nomc{i})] ); end
nom = zeros(1,m);

for i = 1:np,
    l = length(nomc{i})-1;
    nomc{i} = conv(nomc{i},fliplr(nomc{i}));
    nomc{i} = nomc{i}(l+1:2*l+1);
    nom(1:l+1) = nom(1:l+1)+var{i}*nomc{i}/gain{i};
end

as = den;
bs = cor2ma(nom,200*length(nom));


text:
function [errorCode] = pco_sdk_example_swtrig(exposure_time)
% set variables grab and display a single image
%
%   [errorCode] = pco_sdk_example_swtrig(exposure_time)
%
% * Input parameters :
%    exposure_time          camera exposure time (default=10ms)
%
% * Output parameters :
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab single image from a recording pco.camera
%using functions PCO_AddBufferEx and PCO_WaitforBuffer
%display the grabbed image
%
%function workflow
%open camera
%set variables
%start camera
%grab image
%show image
%stop camera
%close camera
%

glvar=struct('do_libunload',1,'do_close',0,'camera_open',0,'out_ptr',[]);

if(~exist('exposure_time','var'))
 exposure_time = 10;
end

%reduce_display_size=1: display only top-left corner 800x600Pixel
reduce_display_size=0;

pco_camera_load_defines();
subfunc=pco_camera_subfunction();
subfunc.fh_lasterr(0);

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

subfunc.fh_stop_camera(out_ptr);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,out_ptr,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
if(errorCode)
 pco_errdisp('PCO_GetCameraDescription',errorCode);
end

bitpix=uint16(cam_desc.wDynResDESC);

cam_type=libstruct('PCO_CameraType');
set(cam_type,'wSize',cam_type.structsize);
[errorCode,out_ptr,cam_type] = calllib('PCO_CAM_SDK', 'PCO_GetCameraType', out_ptr,cam_type);
if(errorCode)
 pco_errdisp('PCO_GetCameraType',errorCode);
end

interface=uint16(cam_type.wInterfaceType);
subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2)
subfunc.fh_set_triggermode(out_ptr,TRIGGER_MODE_SOFTWARETRIGGER);

if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_BINARYANDASCII);
end

%set bitalignment LSB
bitalign=uint16(BIT_ALIGNMENT_LSB);
[errorCode,out_ptr] = calllib('PCO_CAM_SDK', 'PCO_SetBitAlignment', out_ptr,bitalign);
 pco_errdisp('PCO_SetBitAlignment',errorCode);

%save actual RecoderSubmode
prev_rec_submode=uint16(10);
[errorCode,out_ptr,prev_rec_submode] = calllib('PCO_CAM_SDK', 'PCO_GetRecorderSubmode', out_ptr,prev_rec_submode);
pco_errdisp('PCO_GetRecorderSubmode',errorCode);

%set RECORDER_SUBMODE_RING_BUFFER
[errorCode,out_ptr] = calllib('PCO_CAM_SDK', 'PCO_SetRecorderSubmode', out_ptr,RECORDER_SUBMODE_RINGBUFFER);
pco_errdisp('PCO_SetRecorderSubmode',errorCode);

%set highest Pixelrate
index=1;
for n=2:4
 if(cam_desc.dwPixelRateDESC(n)>cam_desc.dwPixelRateDESC(index))
  index=n;
 end
end

[errorCode,out_ptr] = calllib('PCO_CAM_SDK', 'PCO_SetPixelRate', out_ptr,cam_desc.dwPixelRateDESC(index));
if(errorCode)
 pco_errdisp('PCO_SetPixelRate',errorCode);
else
 disp(['PixelRate set to ',int2str(cam_desc.dwPixelRateDESC(index))]);
end

%if PCO_ArmCamera does fail no images can be grabbed
errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_ArmCamera',errorCode);
 ME = MException('PCO_ERROR:ArmCamera','Cannot continue script with not armed camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

%adjust transfer parameter if necessary
subfunc.fh_set_transferparameter(out_ptr);

subfunc.fh_get_triggermode(out_ptr);
subfunc.fh_show_frametime(out_ptr);

metadatasize=0;
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
 wMetaDataMode=uint16(0);
 wMetaDataSize=uint16(0);
 wMetaDataVersion=uint16(0);

 [errorCode,~,wMetaDataMode,wMetaDataSize,wMetaDataVersion] = calllib('PCO_CAM_SDK', 'PCO_GetMetaDataMode',out_ptr,wMetaDataMode,wMetaDataSize,wMetaDataVersion);
 pco_errdisp('PCO_GetMetaDataMode',errorCode);
 if(errorCode)
  pco_errdisp('PCO_GetMetaDataMode',errorCode);
  return;
 end

 disp(['Camera supports Metadata Version: ',num2str(wMetaDataVersion)]);

 if(wMetaDataMode~=0)
  metadatasize=wMetaDataSize;
  disp(['Metadata is enabled. Number of Bytes: ',num2str(metadatasize)]);
 end
end

act_xsize=uint16(0);
act_ysize=uint16(0);
ccd_xsize=uint16(0);
ccd_ysize=uint16(0);

%use PCO_GetSizes because this always returns accurate image size for next recording
[errorCode,out_ptr,act_xsize,act_ysize]  = calllib('PCO_CAM_SDK', 'PCO_GetSizes', out_ptr,act_xsize,act_ysize,ccd_xsize,ccd_ysize);
if(errorCode)
 pco_errdisp('PCO_GetSizes',errorCode);
end

disp(['sizes: horizontal ',int2str(act_xsize),' vertical ',int2str(act_ysize)]);

lineadd=0;
if(metadatasize>0)
 xs=uint32(fix((double(bitpix)+7)/8));
 xs=xs*uint32(act_xsize);
 i=uint16(floor(metadatasize/double(xs)));
 lineadd=i+1;
end

imasize= uint32(fix((double(bitpix)+7)/8)) * uint32(act_ysize+lineadd)* uint32(act_xsize);

%only for firewire add always some lines
%to ensure enough memory is allocated for the transfer
if(interface==INTERFACE_FIREWIRE)
  imas=uint32(fix((double(bitpix)+7)/8));
  imas= imas*uint32(act_ysize+lineadd)* uint32(act_xsize);
  i=floor(double(imas)/4096);
  i=i+1;
  i=i*4096;
  imasize=i;
  i=i-double(imas);
  xs=uint32(fix((double(bitpix)+7)/8));
  xs=xs*uint32(act_xsize);
  i=floor(i/double(xs));
  i=i+1;
  lineadd=lineadd+i;
 disp(['imasize is: ',int2str(imas),' aligned: ',int2str(imasize)]);
end
disp([int2str(lineadd),' additional lines must be allocated ']);

%Allocate data array for images
image_stack=ones(act_xsize,(act_ysize+lineadd),'uint16');

sBufNr=int16(-1);
im_ptr = libpointer('uint16Ptr',image_stack(:,:));
ev_ptr = libpointer('voidPtr');

%Allocate a SDK buffer and set address of buffer from data array image_stack
[errorCode,out_ptr,sBufNr]  = calllib('PCO_CAM_SDK','PCO_AllocateBuffer', out_ptr,sBufNr,imasize,im_ptr,ev_ptr);
if(errorCode)
 pco_errdisp('PCO_AllocateBuffer',errorCode);
 clear image_stack;
 commandwindow;
 return;
end

disp(['bufnr: ',int2str(sBufNr)]);

subfunc.fh_start_camera(out_ptr);

disp('get image');
tic;
%add the allocated buffer
[errorCode]  = calllib('PCO_CAM_SDK','PCO_AddBufferEx', out_ptr,0,0,sBufNr,act_xsize,act_ysize,bitpix);
if(errorCode)
 pco_errdisp('PCO_AddBufferEx',errorCode);
end

trigdone=int16(1);
[errorCode,out_ptr,trigdone]  = calllib('PCO_CAM_SDK','PCO_ForceTrigger',out_ptr,trigdone);
if(errorCode)
 pco_errdisp('PCO_ForceTrigger',errorCode);
else
 disp(['trigger done return: ',int2str(trigdone)]);
end

buflist=libstruct('PCO_Buflist');
buflist.sBufnr_1=uint16(sBufNr);
[errorCode,out_ptr,buflist]  = calllib('PCO_CAM_SDK','PCO_WaitforBuffer', out_ptr,1,buflist,1000);
if(errorCode)
 pco_errdisp('PCO_WaitforBuffer',errorCode);
else
 disp('PCO_WaitforBuffer done');
end

if(errorCode==PCO_NOERROR)
%disp(['statusdll: ',num2str(buflist_1.dwStatusDll,'%08X'),' statusdrv: ',num2str(buflist_1.dwStatusDrv,'%08X')]);
%test and display buffer
 if((bitand(buflist.dwStatusDll_1,hex2dec('00008000')))&&(buflist.dwStatusDrv_1==0))
  disp(['Event for buffer done, StatusDrv ',num2str(buflist.dwStatusDrv_1,'%08X')]);
%show image
  ima=im_ptr.Value;
  if(reduce_display_size~=0)
   [xs,ys]=size(ima);
   xmax=800;
   ymax=600;
   if((xs>xmax)&&(ys>ymax))
    ima=ima(1:xmax,1:ymax);
   elseif(xs>xmax)
    ima=ima(1:xmax,:);
   elseif(ys>ymax)
    ima=ima(:,1:ymax);
   end
  end
  ima=ima';
  m=max(max(ima(10:end-10,10:end-10)));
  imshow(ima,[0,m+100]);
  disp(['found max ',int2str(m)]);
 end
end

subfunc.fh_stop_camera(out_ptr);

%this will remove all pending buffers in the queue
[errorCode] = calllib('PCO_CAM_SDK', 'PCO_CancelImages', out_ptr);
if(errorCode)
 pco_errdisp('PCO_CancelImages',errorCode);
end

disp('Press "Enter" to close window and proceed')
pause();
close();

%set changed values back
subfunc.fh_set_triggermode(out_ptr,TRIGGER_MODE_AUTOTRIGGER);
%set saved RecoderSubmode
[errorCode] = calllib('PCO_CAM_SDK', 'PCO_SetRecorderSubmode', out_ptr,prev_rec_submode);
if(errorCode)
 pco_errdisp('PCO_SetRecorderSubmode',errorCode);
end

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

[errorCode]  = calllib('PCO_CAM_SDK','PCO_FreeBuffer',out_ptr,sBufNr);
if(errorCode)
 pco_errdisp('PCO_FreeBuffer',errorCode);
else
 disp('PCO_FreeBuffer done ');
end

subfunc.fh_stop_camera(out_ptr);

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 clearvars -except ME glvar;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  disp(txt);
  warning('off','backtrace')
  warning(msg)
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
 end
  clearvars;
  commandwindow;
  return;
 else
  clearvars -except ME;
  rethrow(ME)
 end
end

clearvars -except glvar;

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

clearvars glvar;
commandwindow;
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab