text:
function [preprocessSettingStruct, preprocessingSettingsAll] = getRegistrationSettings(obj,inputTitleStr,varargin)
	% Internal function to get user imaging preprocessing settings.
	% Biafra Ahanonu
	% started: 2016.05.31 [14:53:00]
	% inputs
		% inputTitleStr - string title to show user.
	% outputs
		% preprocessSettingStruct - structure containing user settings.
		% preprocessingSettingsAll - structure of the final pre-processing settings.

	% changelog
		% 2019.12.07 [16:23:36] - Added tooltips. Made options list easier on the eyes with improved spacing. Added a callback for user selecting the reference turboreg frame and or adding arbitrary values to a number of other options.
		% 2019.12.07 [17:46:03] - Change how settings are initialized, should make easier to maintain and add new settings.
		% 2019.12.08 [22:34:38] - Added additional tooltips and checks to make sure user inputs correct data type/size.
		% 2019.12.08 [22:49:48] - Allow users to input previous preprocessing settings.
		% 2020.04.18 [19:27:24] - Make sure any modifications from default are propagated upon multiple calls to getRegistrationSettings.
		% 2020.04.19 [15:20:48] - Fix how modified values stored so only store the modified value and not all the options at time of modification.
		% 2020.05.28 [21:07:27] - Added support for nargin=1.
		% 2020.10.21 [16:52:06] - Add support for user canceling the input.
		% 2020.10.24 [18:30:56] - Added support for calculating dropped frames if entire frame of a movie is a set value.
		% 2021.06.21 [14:27:26] - Switch to single page support.
		% 2021.06.21 [20:52:36] - Multi-column layout since most computers are widescreen now.
		% 2021.06.29 [12:52:36] - Added regRegionUseBtwnSessions support and made fix for loading prior settings if they do not contain new settings.
		% 2021.07.06 [09:42:09] - Switch to callbacks to close instead of using pause, allows for more flexibility going forward and use with MATLAB Online.
		% 2021.07.13 [17:55:30] - Added support for certain options (e.g. motion correction) to give a vector
		% 2021.07.22 [12:10:10] - Added movie detrend options.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.12.01 [19:21:00] - Handle "User selects specific value" going into infinite loop. Allow cancel option to prevent loop.
		% 2022.06.28 [15:28:13] - Add additional black background support.
		% 2022.07.10 [17:21:21] - Added largeMovieLoad setting.
	% TODO
		% DONE: Allow user to input prior settings, indicate those changed from default by orange or similar color.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Struct: Empty = no previous settings, else a structure of preprocessingSettingsAll from previous run.
	options.inputSettings = [];
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	if nargin==1
		disp('No title string given, using blank');
		inputTitleStr = '';
	end

	% Initialize constants
	nWorkersDefault = java.lang.Runtime.getRuntime().availableProcessors-1;
	userSelectStr = 'User selects specific value.';
	userSelectVal = 0;
	defaultTooltips = 'NO TIPS FOR YOU';

	vectorAllowList = {'motionCorrectionRefFrame'};

	% Options that allow users to manually set the value
	userOptionsAllowManualSelect = {...
	'numTurboregIterations',...
	'turboregNumFramesSubset',...
	'pxToCrop',...
	'motionCorrectionRefFrame',...
	'SmoothX',...
	'SmoothY',...
	'filterBeforeRegFreqLow',...
	'filterBeforeRegFreqHigh',...
	'medianFilterSize',...
	'downsampleFactorTime',...
	'downsampleFactorSpace',...
	'stripSize',...
	'stripfreqLowExclude',...
	'stripfreqHighExclude',...
	'filterBeforeRegImagejFFTLarge',...
	'filterBeforeRegImagejFFTSmall',...
	'loadMovieInEqualParts',...
	'nParallelWorkers',...
	'calcDroppedFramesFromMovie'...
	};

	% Create structure with options.
	tS = [];

	% ===================================
	tS.MAIN______________ = [];
		tS.MAIN______________.val = {{'====================='}};
		tS.MAIN______________.str = {{'====================='}};
		tS.MAIN______________.tooltip =  {{'====================='}};
	tS.checkConcurrentAnalysis = [];
		tS.checkConcurrentAnalysis.val = {{0,1}};
		tS.checkConcurrentAnalysis.str = {{'DO NOT check for analysis across workstations','DO check for analysis across workstations'}};
		tS.checkConcurrentAnalysis.tooltip =  {{'Select to allow multiple workstations to analyze several of the same imaging sessions WITHOUT duplicating analysis.'}};
	tS.resetParallelPool = [];
		tS.resetParallelPool.val = {{0,1}};
		tS.resetParallelPool.str = {{'DO NOT restart parallel pool after each folder is analyzed','DO restart parallel pool after each folder is analyzed'}};
		tS.resetParallelPool.tooltip =  {{['Parallel workers sometimes start to use a lot of memory.' 10 'Restart them to reset to lower memory usage.']}};
	tS.analyzeFromDisk = [];
		tS.analyzeFromDisk.val = {{0,1}};
		tS.analyzeFromDisk.str = {{'Load entire movie into RAM','Use disk, load part of movie into RAM'}};
		tS.analyzeFromDisk.tooltip =  {{defaultTooltips}};
	tS.calcDroppedFramesFromMovie = [];
		tS.calcDroppedFramesFromMovie.val = {{[],0,1,NaN}};
		tS.calcDroppedFramesFromMovie.str = {{'Calculate dropped frames from metadata.','All 0 frame is a dropped frame.','All 1 frame is a dropped frame.','All NaN frame is a dropped frame.'}};
		tS.calcDroppedFramesFromMovie.tooltip =  {{'If metadata does not contain dropped frame information, set this if specific frames are set to a single value to indicate dropped frames.'}};

	% ===================================
	tS.REGISTRATION______________ = [];
		tS.REGISTRATION______________.val = {{'====================='}};
		tS.REGISTRATION______________.str = {{'====================='}};
		tS.REGISTRATION______________.tooltip =  {{'====================='}};
	tS.mcMethod = [];
		tS.mcMethod.val = {{'turboreg','normcorre'}};
		tS.mcMethod.str = {{'Turboreg motion correction','NoRMCorre (non-rigid) motion correction'}};
		tS.mcMethod.tooltip =  {{'Method to use for motion correction.'}};
	tS.parallel = [];
		tS.parallel.val = {{1,0}};
		tS.parallel.str = {{'parallel processing','NO parallel processing'}};
		tS.parallel.tooltip =  {{'Use parallel processing during motion correction.'}};
	tS.registrationFxn = [];
		tS.registrationFxn.val = {{'transfturboreg','imwarp','imtransform'}};
		tS.registrationFxn.str = {{'transfturboreg','imwarp','imtransform'}};
		tS.registrationFxn.tooltip =  {{'Keep default unless have issues with transfturboreg (e.g. random frames look dim) then use MATLAB registration functions (imwarp or imtransform).'}};
	tS.turboregRotation = [];
		tS.turboregRotation.val = {{0,1}};
		tS.turboregRotation.str = {{'DO NOT turboreg rotation','DO turboreg rotation'}};
		tS.turboregRotation.tooltip =  {{'Unless you have rotation in your movie, leave OFF. Produces better results.'}};
	tS.RegisType = [];
		tS.RegisType.val = {{1,3,4,2}};
		tS.RegisType.str = {{'affine (1, parallelism maintained, no skew or rotation)','projective (3, parallelism not guaranteed, rotation allowed)','affine (4, parallelism maintained, skew allowed)','affine (2, parallelism maintained, rotation allowed)'}};
		tS.RegisType.tooltip =  {{'"affine" is good for most rigid cases, use "projective" if movie has complex motion'}};
	tS.numTurboregIterations = [];
		tS.numTurboregIterations.val = {{1,userSelectVal,2,3,4,5}};
		tS.numTurboregIterations.str = {{1,userSelectStr,2,3,4,5}};
		tS.numTurboregIterations.tooltip =  {{defaultTooltips}};
	tS.turboregNumFramesSubset = [];
		tS.turboregNumFramesSubset.val = {{5000,userSelectVal,1,5,10,50,100,250,500,1000,2000,3000,5000,10000,15000}};
		tS.turboregNumFramesSubset.str = {{5000,userSelectStr,1,5,10,50,100,250,500,1000,2000,3000,5000,10000,15000}};
		tS.turboregNumFramesSubset.tooltip =  {{'To save memory, motion correction is run on the indicated # of frames from the data.'}};
	tS.pxToCrop = [];
		tS.pxToCrop.val = {{14,userSelectVal,15,16,17,18,19,20,21,22,23,24,25}};
		tS.pxToCrop.str = {{14,userSelectStr,15,16,17,18,19,20,21,22,23,24,25}};
		tS.pxToCrop.tooltip =  {{'Maximum # of pixels to convert to NaNs around the border after motion correction.'}};
	tS.motionCorrectionRefFrame = [];
		tS.motionCorrectionRefFrame.val = {{100,userSelectVal,1,10,100,1000}};
		tS.motionCorrectionRefFrame.str = {{100,userSelectStr,1,10,100,1000}};
		tS.motionCorrectionRefFrame.tooltip =  {{['The reference frame for motion correction.' 10 'Avoid selecting 1st couple of movie frames or one that has issues.']}};
	tS.regRegionUseBtwnSessions = [];
		% tS.regRegionUseBtwnSessions.val = {{-1,0,-2,0}};
		tS.regRegionUseBtwnSessions.val = {{1,2,3,4}};
		tS.regRegionUseBtwnSessions.str = {{'NO | do not duplicate area coords across multiple folders','YES | duplicate area coords across multiple folders','YES | duplicate area coords if subject (animal) the same','YES | duplicate area coords across ALL folders'}};
		tS.regRegionUseBtwnSessions.tooltip =  {{['Between session motion correction area coordinates (area used to get registration translation coordinates).']}};

	% ===================================
	tS.REGISTRATION_NORMALIZATION______________ = [];
		tS.REGISTRATION_NORMALIZATION______________.val = {{'====================='}};
		tS.REGISTRATION_NORMALIZATION______________.str = {{'====================='}};
		tS.REGISTRATION_NORMALIZATION______________.tooltip =  {{'====================='}};
	tS.normalizeMeanSubtract = [];
		tS.normalizeMeanSubtract.val = {{1,0}};
		tS.normalizeMeanSubtract.str = {{'normalize movie before turboreg','do not normalize movie before turboreg'}};
		tS.normalizeMeanSubtract.tooltip = {{['Critical for most movies, normalizes each frame to improve motion correction.' 10 'DO NOT CHANGE unless you know what you are doing.']}};
	tS.normalizeMeanSubtractNormalize = [];
		tS.normalizeMeanSubtractNormalize.val = {{1,0}};
		tS.normalizeMeanSubtractNormalize.str = {{'subtract mean per frame','do not subtract mean per frame'}};
		tS.normalizeMeanSubtractNormalize.tooltip = {{'This subtracts the mean of each frame. Keep enabled is nearly all cases.'}};
	tS.normalizeComplementMatrix = [];
		tS.normalizeComplementMatrix.val = {{1,0}};
		tS.normalizeComplementMatrix.str = {{'invert movie before turboreg','DO NOT invert movie before turboreg'}};
		tS.normalizeComplementMatrix.tooltip = {{'Inverts all movie values to give more weight to dark features (e.g. blood vessels).'}};
	tS.normalizeType = [];
		tS.normalizeType.val = {{'bandpass','divideByLowpass','imagejFFT','highpass','matlabDisk'}};
		tS.normalizeType.str = {{'bandpass','divideByLowpass','imagejFFT','highpass','matlabDisk'}};
		tS.normalizeType.tooltip = {{['Spatial filtering applied before getting spatial translation coordinates.' 10 'Try "matlabDisk" if default does not work.']}};
	tS.normalizeFreqLow = [];
		tS.normalizeFreqLow.val = {{70,10,20,30,40,50,60,70,80,90}};
		tS.normalizeFreqLow.str = {{70,10,20,30,40,50,60,70,80,90}};
		tS.normalizeFreqLow.tooltip = {{'For any "normalizeType" FFT options selected, the lower frequency.'}};
	tS.normalizeFreqHigh = [];
		tS.normalizeFreqHigh.val = {{100,80,90,100,110}};
		tS.normalizeFreqHigh.str = {{100,80,90,100,110}};
		tS.normalizeFreqHigh.tooltip = {{'For any "normalizeType" FFT options selected, the higher frequency.'}};
	tS.normalizeBandpassType = [];
		tS.normalizeBandpassType.val = {{'bandpass','lowpass','highpass'}};
		tS.normalizeBandpassType.str = {{'bandpass','lowpass','highpass'}};
		tS.normalizeBandpassType.tooltip = {{'If "normalizeType" is an FFT, the type'}};
	tS.normalizeBandpassMask = [];
		tS.normalizeBandpassMask.val = {{'gaussian','binary'}};
		tS.normalizeBandpassMask.str = {{'gaussian','binary'}};
		tS.normalizeBandpassMask.tooltip = {{['If "normalizeType" is an FFT, how "normalizeBandpassType" is applied.' 10 '"binary" option will produce ringing.']}};
	tS.SmoothX = [];
		tS.SmoothX.val = {{10,userSelectVal,1,5,10,20,30,40,50,60,70,80,90}};
		tS.SmoothX.str = {{10,userSelectStr,1,5,10,20,30,40,50,60,70,80,90}};
		tS.SmoothX.tooltip = {{'Turboreg''s x smoothing size in pixels.'}};
	tS.SmoothY = [];
		tS.SmoothY.val = {{10,userSelectVal,1,5,10,20,30,40,50,60,70,80,90}};
		tS.SmoothY.str = {{10,userSelectStr,1,5,10,20,30,40,50,60,70,80,90}};
		tS.SmoothY.tooltip = {{'Turboreg''s y smoothing size in pixels.'}};
	tS.zapMean = [];
		tS.zapMean.val = {{0,1}};
		tS.zapMean.str = {{0,1}};
		tS.zapMean.tooltip = {{'Turboreg''s function to remove the mean. Leave disabled since this is already done outside turboreg.'}};

	% ===================================
	tS.MOVIE_NORMALIZATION______________ = [];
		tS.MOVIE_NORMALIZATION______________.val = {{'====================='}};
		tS.MOVIE_NORMALIZATION______________.str = {{'====================='}};
		tS.MOVIE_NORMALIZATION______________.tooltip =  {{'====================='}};
	tS.filterBeforeRegister = [];
		tS.filterBeforeRegister.val = {{[],'divideByLowpass','imagejFFT','bandpass'}};
		tS.filterBeforeRegister.str = {{'NO filtering before registering','matlab divide by lowpass before registering','imageJ divide by lowpass (requires Miji!)','matlab bandpass before registering'}};
		tS.filterBeforeRegister.tooltip = {{'IMPORTANT: type of spatial filtering applied after getting spatial transformation coordinates but BEFORE actually motion correcting the movie.'}};
	tS.saveBeforeFilterRegister = [];
		tS.saveBeforeFilterRegister.val = {{0,1}};
		tS.saveBeforeFilterRegister.str = {{'DO NOT SAVE registered movie sans spatial filtering','DO SAVE registered movie sans spatial filtering'}};
		tS.saveBeforeFilterRegister.tooltip = {{['Save registered movie before spatial filtering applied.' 10 'This is a good way to check the correlation of all frames to mean frame to show improvement with motion correction.']}};
	tS.saveFilterBeforeRegister = [];
		tS.saveFilterBeforeRegister.val = {{[],'save'}};
		tS.saveFilterBeforeRegister.str = {{'DO NOT save lowpass movie.','DO save lowpass movie.'}};
		tS.saveFilterBeforeRegister.tooltip = {{'Save a low-pass version of the movie, e.g. if looking at neuropil activity.'}};
	tS.filterBeforeRegFreqLow = [];
		tS.filterBeforeRegFreqLow.val = {{0,userSelectVal,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400,500,600}};
		tS.filterBeforeRegFreqLow.str = {{0,userSelectStr,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400,500,600}};
		tS.filterBeforeRegFreqLow.tooltip = {{['Value is in cycles per movie x-y dimensions.' 10 'If divideByLowpass, leave this alone.' 10 'If "matlab bandpass" then set to the lowest spatial frequency to keep (e.g. large objects).']}};
	tS.filterBeforeRegFreqHigh = [];
		tS.filterBeforeRegFreqHigh.val = {{10,userSelectVal,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400,500,600}};
		tS.filterBeforeRegFreqHigh.str = {{10,userSelectStr,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400,500,600}};
		tS.filterBeforeRegFreqHigh.tooltip = {{['Value is in cycles per movie x-y dimensions.' 10 'If divideByLowpass, set to the highest spatial frequency to exclude.' 10 'If "matlab bandpass" then set to the highest spatial frequency to keep (e.g. small ibjects).']}};
	tS.filterBeforeRegImagejFFTLarge = [];
		tS.filterBeforeRegImagejFFTLarge.val = {{10000,userSelectVal,100,500,1000,5000,8000}};
		tS.filterBeforeRegImagejFFTLarge.str = {{10000,userSelectStr,100,500,1000,5000,8000}};
		tS.filterBeforeRegImagejFFTLarge.tooltip = {{'Generally ignore since ImageJ isn''t needed for spatial filtering anymore'}};
	tS.filterBeforeRegImagejFFTSmall = [];
		tS.filterBeforeRegImagejFFTSmall.val = {{80,userSelectVal,10,20,30,40,50,60,70,90,100}};
		tS.filterBeforeRegImagejFFTSmall.str = {{80,userSelectStr,10,20,30,40,50,60,70,90,100}};
		tS.filterBeforeRegImagejFFTSmall.tooltip = {{'Generally ignore since ImageJ isn''t needed for spatial filtering anymore'}};
	tS.medianFilterSize = [];
		tS.medianFilterSize.val = {{3,userSelectVal,5,7,9,11,13,15}};
		tS.medianFilterSize.str = {{3,userSelectStr,5,7,9,11,13,15}};
		tS.medianFilterSize.tooltip = {{['The size in pixels of the median filter.' 10 'If "medianFilter" is selected for preprocessing.']}};

	% ===================================
	tS.MOVIE_DETREND______________ = [];
		tS.MOVIE_DETREND______________.val = {{'====================='}};
		tS.MOVIE_DETREND______________.str = {{'====================='}};
		tS.MOVIE_DETREND______________.tooltip =  {{'====================='}};
	tS.detrendDegree = [];
		tS.detrendDegree.val = {{1,userSelectVal,2,3,4,5,6}};
		tS.detrendDegree.str = {{1,userSelectStr,2,3,4,5,6}};
		tS.detrendDegree.tooltip = {{'For "detrend" function to compensate for photobleaching, the nth-degree polynomial trend to fit. Usually 1 (linear fit) if sufficient.'}};

	% ===================================
	tS.MOVIE_DOWNSAMPLING______________ = [];
		tS.MOVIE_DOWNSAMPLING______________.val = {{'====================='}};
		tS.MOVIE_DOWNSAMPLING______________.str = {{'====================='}};
		tS.MOVIE_DOWNSAMPLING______________.tooltip =  {{'====================='}};
	tS.downsampleFactorTime = [];
		tS.downsampleFactorTime.val = {{4,userSelectVal,1,2,4,6,8,10,20}};
		tS.downsampleFactorTime.str = {{4,userSelectStr,1,2,4,6,8,10,20}};
		tS.downsampleFactorTime.tooltip = {{'By what factor to downsample movie in time.'}};
	tS.downsampleFactorSpace = [];
		tS.downsampleFactorSpace.val = {{2,userSelectVal,1,2,4,6,8,10,20}};
		tS.downsampleFactorSpace.str = {{2,userSelectStr,1,2,4,6,8,10,20}};
		tS.downsampleFactorSpace.tooltip = {{'By what factor to downsample movie in space.'}};

	% ===================================
	tS.IO_and_MOVIE_IDENTIFICATION______________ = [];
		tS.IO_and_MOVIE_IDENTIFICATION______________.val = {{'====================='}};
		tS.IO_and_MOVIE_IDENTIFICATION______________.str = {{'====================='}};
		tS.IO_and_MOVIE_IDENTIFICATION______________.tooltip =  {{'====================='}};
	tS.inputDatasetName = [];
		tS.inputDatasetName.val = {{obj.inputDatasetName,'/1','/Movie','/movie','/images','/Data/Images','/data'}};
		tS.inputDatasetName.str = {{obj.inputDatasetName,'/1','/Movie','/movie','/images','/Data/Images','/data'}};
		tS.inputDatasetName.tooltip = {{'HDF5 dataset name of input data.'}};
	tS.outputDatasetName = [];
		tS.outputDatasetName.val = {{obj.outputDatasetName,'/1','/Movie','/movie','/images'}};
		tS.outputDatasetName.str = {{obj.outputDatasetName,'/1','/Movie','/movie','/images'}};
		tS.outputDatasetName.tooltip = {{'HDF5 dataset name of output data.'}};
	tS.fileFilterRegexp = [];
		tS.fileFilterRegexp.val = {{obj.fileFilterRegexpRaw,'concat_.*.h5','concatenated_.*.h5','crop.*.h5','recording.*.tif','concat.*.tif','dfstd_.*.h5','dfof_.*.h5','concat'}};
		tS.fileFilterRegexp.str = {{obj.fileFilterRegexpRaw,'concat_.*.h5','concatenated_.*.h5','crop.*.h5','recording.*.tif','concat.*.tif','dfstd_.*.h5','dfof_.*.h5','concat'}};
		tS.fileFilterRegexp.tooltip = {{'Regular expression used to find movie data files in each folder.'}};
	tS.processMoviesSeparately = [];
		tS.processMoviesSeparately.val = {{0,1}};
		tS.processMoviesSeparately.str = {{'No','Yes'}};
		tS.processMoviesSeparately.tooltip = {{['"No", all files in folder concatenated in natural sort order.' 10 'If "Yes" then if folder contains multiple movie files, they will be processed separately instead of concatenated together.']}};
	tS.loadMoviesFrameByFrame = [];
		tS.loadMoviesFrameByFrame.val = {{0,1}};
		tS.loadMoviesFrameByFrame.str = {{'No','Yes'}};
		tS.loadMoviesFrameByFrame.tooltip = {{['"No" allows movies to be loaded in one chunk.' 10 '"Yes" loads movies frame-by-frame, which can potentially improve memory usage.']}};
	tS.treatMoviesAsContinuousSwitch = [];
		tS.treatMoviesAsContinuousSwitch.val = {{1,0}};
		tS.treatMoviesAsContinuousSwitch.str = {{'Yes','No'}};
		tS.treatMoviesAsContinuousSwitch.tooltip = {{['"Yes" treats multiple movie files in a folder as one continuous movie.' 10 '"No" treats each movie file in a folder as separate movies.']}};
	tS.loadMovieInEqualParts = [];
		tS.loadMovieInEqualParts.val = {{0,userSelectVal,1,2,3,4,5,6,7,8,9,10}};
		tS.loadMovieInEqualParts.str = {{0,userSelectStr,1,2,3,4,5,6,7,8,9,10}};
		tS.loadMovieInEqualParts.tooltip = {{'Number of 50-frame segments evenly spaced in time to load.'}};
	tS.useParallel = [];
		tS.useParallel.val = {{1,0}};
		tS.useParallel.str = {{'Yes','No'}};
		tS.useParallel.tooltip = {{'Whether to use parallelization in processing steps.'}};
	tS.nParallelWorkers = [];
		% tS.nParallelWorkers.val = {[nWorkersDefault, userSelectVal, mat2cell([1:(nWorkersDefault*2)],1,ones(1,2*nWorkersDefault))]'};
		% tS.nParallelWorkers.str = {[nWorkersDefault, userSelectStr, mat2cell([1:(nWorkersDefault*2)],1,ones(1,2*nWorkersDefault))]'};
		tS.nParallelWorkers.val = {{nWorkersDefault, userSelectVal, nWorkersDefault*2}};
		tS.nParallelWorkers.str = {{nWorkersDefault, userSelectStr, nWorkersDefault*2}};
		tS.nParallelWorkers.tooltip = {{defaultTooltips}};
	tS.largeMovieLoad = [];
		tS.largeMovieLoad.val = {{0,1}};
		tS.largeMovieLoad.str = {{'No','Yes'}};
		tS.largeMovieLoad.tooltip = {{'Whether to load movie without pre-allocating matrix, useful for single large movie in each folder that needs to be processed.'}};

	% ===================================
	tS.STRIPE_REMOVAL______________ = [];
		tS.STRIPE_REMOVAL______________.val = {{'====================='}};
		tS.STRIPE_REMOVAL______________.str = {{'====================='}};
		tS.STRIPE_REMOVAL______________.tooltip =  {{'====================='}};
	tS.stripOrientationRemove = [];
		tS.stripOrientationRemove.val = {{'none','vertical','horizontal','both'}};
		tS.stripOrientationRemove.str = {{'none','vertical','horizontal','both'}};
		tS.stripOrientationRemove.tooltip = {{'Direction of stripes to remove from movie.'}};
	tS.stripSize = [];
		tS.stripSize.val = {{7,userSelectVal,1,3,5,7,9,11,13,15}};
		tS.stripSize.str = {{7,userSelectStr,1,3,5,7,9,11,13,15}};
		tS.stripSize.tooltip = {{['Width in pixels in FFT domain of filter.' 10 'Higher values allow more tolerance for strips that are not exactly vertical or horizontal.']}};
	tS.stripfreqLowExclude = [];
		tS.stripfreqLowExclude.val = {{20,userSelectVal,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400}};
		tS.stripfreqLowExclude.str = {{20,userSelectStr,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400}};
		tS.stripfreqLowExclude.tooltip = {{'Lowest frequency of stripe to exclude from strip filter.'}};
	tS.stripfreqHighExclude = [];
		tS.stripfreqHighExclude.val = {{20,userSelectVal,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400}};
		tS.stripfreqHighExclude.str = {{20,userSelectStr,1,2,3,4,5,7,10,15,20,25,30,35,40,45,50,60,70,80,100,150,200,250,300,350,400}};
		tS.stripfreqHighExclude.tooltip = {{'Highest frequency of stripe to exclude from strip filter.'}};
	tS.stripfreqBandpassType = [];
		tS.stripfreqBandpassType.val = {{'highpass','bandpass','lowpass'}};
		tS.stripfreqBandpassType.str = {{'highpass','bandpass','lowpass'}};
		tS.stripfreqBandpassType.tooltip = {{'Type of stripe bandpass to use.'}};

	if ~isempty(options.inputSettings)
		try
			nonDefaultProperties = {};
			propertyListTmp = fieldnames(tS);
			nProperties = size(propertyListTmp,1);
			for propertyNo = 1:nProperties
				property = char(propertyListTmp(propertyNo));

				% Check if value had been modified from default, if so then add that as the 1st, default option.
				% if isequal(options.inputSettings.(property).val{1},tS.(property).val{1})
				if isfield(options.inputSettings,property)
					if isfield(options.inputSettings.(property),'modified')
						nonDefaultProperties{end+1} = property;
						tS.(property).val = {[{options.inputSettings.(property).modified.val},tS.(property).val{1}]};
						tS.(property).str = {[options.inputSettings.(property).modified.str,tS.(property).str{1}]};

						% Ensure modification from default is passed onto future edits
						tS.(property).modified = [];
						% tS.(property).modified.val = tS.(property).val{1};
						% tS.(property).modified.str = tS.(property).str{1};
						tS.(property).modified.val = options.inputSettings.(property).modified.val;
						tS.(property).modified.str = options.inputSettings.(property).modified.str;
					end
				end
			end
			options.inputSettings
		catch err
			disp(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			disp(repmat('@',1,7))
		end
	else
		nonDefaultProperties = {};
	end

	% propertySettings = turboregSettingDefaults;
	preprocessingSettingsAll = tS;

	propertyList = fieldnames(preprocessingSettingsAll);
	nPropertiesToChange = size(propertyList,1);

	% add current property to the top of the list
	for propertyNo = 1:nPropertiesToChange
		property = char(propertyList(propertyNo));
		% propertyOptions = turboregSettingStr.(property);
		propertyOptions = preprocessingSettingsAll.(property).str{1};
		propertySettingsStr.(property) = propertyOptions;
		% propertySettingsStr.(property);
	end

	modX = 1.05;
	modX_text = 1.1;
	figNoDefault = 1337;
	uiListHandles = {};
	uiTextHandles = {};
	uiXIncrement = 0.025/modX;
	uiXOffset = 0.02/modX;
	uiYOffset = 0.90;
	uiYOffsetOrig = uiYOffset;
	uiTxtSize = 0.4;
	uiBoxSize = 0.55;
	uiFontSize = 9;
	nGuiSubsets = 3;
	% subsetList = round(linspace(1,nPropertiesToChange,nGuiSubsets));
	subsetList = [1 35 nPropertiesToChange];
	propNoSwitch = 36;
	% subsetList
	% for subsetNo = 1:nGuiSubsets
	nGuiSubsetsTrue = (nGuiSubsets-1);

	nGuiSubsets = 2;
	nGuiSubsetsTrue = (nGuiSubsets-1);
	% nGuiSubsetsTrue = 1;
	figure(figNoDefault);
	set(gcf,'color',[0 0 0]);
	for thisSet = 1:nGuiSubsetsTrue
		% 1:nPropertiesToChange
		% subsetStartTime = tic;
		subsetStartIdx = subsetList(thisSet);
		subsetEndIdx = subsetList(thisSet+1);

		subsetStartIdx = 1;
		subsetEndIdx = nPropertiesToChange;

		if thisSet==nGuiSubsetsTrue
			propertySubsetList = subsetStartIdx:(subsetEndIdx);
		else
			propertySubsetList = subsetStartIdx:(subsetEndIdx-1);
		end

		[figHandle figNo] = openFigure(figNoDefault, '');
		clf
		uicontrol('Style','Text','String',inputTitleStr,'Units','normalized','Position',[uiXOffset uiYOffset+0.05 0.8 0.05],'BackgroundColor','black','ForegroundColor','white','HorizontalAlignment','Left','FontSize',uiFontSize);
		uicontrol('Style','Text','String',sprintf('Options page %d/%d: Mouse over each options for tips. To continue, press enter. Orange = previous non-default settings.\n>>>On older MATLAB versions, select the command window before pressing enter.',thisSet,nGuiSubsets-1),'Units','normalized','Position',[uiXOffset uiYOffset+0.02 uiTxtSize+uiBoxSize 0.05],'BackgroundColor','black','ForegroundColor','white','HorizontalAlignment','Left','FontSize',uiFontSize);

		propNoDisp = 1;
		for propertyNo = propertySubsetList
			property = char(propertyList(propertyNo));
			% propertyTooltip = turboregSettingTooltips.(property);
			propertyTooltip = char(preprocessingSettingsAll.(property).tooltip{1});

			if propertyNo==propNoSwitch
				propNoDisp = 1;
				uiYOffset = uiYOffsetOrig;
			end
			if propertyNo>=propNoSwitch
				colNo = 49/100;
				uiTxtSize = 0.4/2;
				uiBoxSize = 0.55/2;
			else
				colNo = 0/100;
				uiTxtSize = 0.4/2;
				uiBoxSize = 0.55/2;
			end

			% disp([num2str(propertyNo) ' | ' property])
			if isempty(intersect(propertyNo,[1,propNoSwitch]))
				if isempty(regexp(property,'______________'))
					spaceMod = 0.00;
					FontWeightH = 'normal';
				else
					spaceMod = 0.02;
					uiYOffset = uiYOffset-spaceMod;
					FontWeightH = 'bold';
				end
			else
				FontWeightH = 'bold';
			end

			textPos = [uiXOffset+colNo uiYOffset-uiXIncrement*propNoDisp+modX_text*(0.027/modX) uiTxtSize 0.0225/modX];
			uiTextHandles{propertyNo} = uicontrol('Style','text','String',[property '' 10],'Units','normalized',...
				'Position',textPos,...
				'BackgroundColor',[0.9 0.9 0.9],'ForegroundColor','black','HorizontalAlignment','Left','FontSize',uiFontSize,'ToolTip',propertyTooltip,'FontWeight',FontWeightH);

			listPos = [uiXOffset+uiTxtSize+colNo uiYOffset-uiXIncrement*propNoDisp uiBoxSize 0.05];
			uiListHandles{propertyNo} = uicontrol('Style', 'popup','String', propertySettingsStr.(property),'Units','normalized',...
				'Position', listPos,'BackgroundColor','black','ForegroundColor','white',...
				'Callback',@subfxnInterfaceCallback,'FontSize',uiFontSize,'Tag',property);

			% jEdit = findjobj(uiTextHandles{propertyNo});
			% lineColor = java.awt.Color(1,0,0);  % =red
			% thickness = 3;  % pixels
			% roundedCorners = true;
			% newBorder = javax.swing.border.LineBorder(lineColor,thickness,roundedCorners);
			% jEdit.Border = newBorder;
			% jEdit.repaint;  % redraw the modified control
			% uiTextHandles{propertyNo}.Enable = 'Inactive';
			% optionCallback = ['set(uiListHandles{propertyNo}, ''Backgroundcolor'', ''g'')'];
			% uiListHandles{propertyNo} = uicontrol('Style', 'popup','String', propertySettingsStr.(property),'Units','normalized','Position', [uiXOffset+uiTxtSize uiYOffset-uiXIncrement*propNoDisp uiBoxSize 0.05],'Callback',@(hObject,callbackdata){set(hObject, 'Backgroundcolor', [208,229,180]/255);},'FontSize',uiFontSize);
			% ,'ToolTip',propertyTooltip

			% If property is non-default, set to orange to alert user.
			if any(ismember(nonDefaultProperties,property))
				set(uiListHandles{propertyNo},'Backgroundcolor',[254 216 177]/255/2);
			end

			propNoDisp = propNoDisp+1;
		end

		startMethodHandle = uicontrol('style','pushbutton','Units', 'normalized','position',[80 94 20 2]/100,'FontSize',9,'string','Click when done.','BackgroundColor',[153 255 153]/255,'callback',@subfxn_closeOptions);
		emptyBox = uicontrol('Style','Text','String',[''],'Units','normalized','Position',[1 1 1 1]/100,'BackgroundColor','black','ForegroundColor','white','HorizontalAlignment','Left','FontWeight','bold','FontSize',7);
		set(gcf,'KeyPressFcn', @subfxn_closeOptionsKey);
		ciapkg.view.changeFont('none','fontColor','w')
		% waitfor(gcf,'Tag');
		waitfor(emptyBox);
		% pause
		uiYOffset = 0.90;

		for propertyNo = propertySubsetList
			property = char(propertyList(propertyNo));
			uiListHandleData = get(uiListHandles{propertyNo});
			if isempty(regexp(property,'______________'))
				preprocessSettingStruct.(property) = preprocessingSettingsAll.(property).val{1}{uiListHandleData.Value};
			else
				preprocessSettingStruct.(property) = preprocessingSettingsAll.(property).val{1};
			end
			% preprocessSettingStruct.(property) = turboregSettingDefaults.(property){uiListHandleData.Value};
		end
	end
	close(1337)

	% ensure rotation setting matches appropriate registration type
	if preprocessSettingStruct.turboregRotation==1
		preprocessSettingStruct.RegisType = 3;
	end

	preprocessSettingStruct.refCropFrame = preprocessSettingStruct.motionCorrectionRefFrame;

	% if preprocessSettingStruct.refCropFrame==0
	% 	movieSettings = inputdlg({...
	% 			'frame to reference to: '...
	% 		},...
	% 		'view movie settings',1,...
	% 		{...
	% 			'100'...
	% 		}...
	% 	);
	% 	preprocessSettingStruct.refCropFrame = str2num(movieSettings{1});
	% end
	function subfxnInterfaceCallback(hObject,callbackdata)
		set(hObject, 'Backgroundcolor', [208,229,180]/255/2);

		% De-select the current option, allows user to press enter to continue.
		set(hObject, 'Enable', 'off');
		drawnow;
		set(hObject, 'Enable', 'on');

		% If user asks to manually specify motion correction frame
		if any(strcmp(userOptionsAllowManualSelect,get(hObject,'Tag')))==1
			thisProperty = userOptionsAllowManualSelect{strcmp(userOptionsAllowManualSelect,get(hObject,'Tag'))};
			gValue = get(hObject,'Value');
			gString = get(hObject,'String');
			gTooltip = char(get(hObject,'ToolTip'));
			gProperty = get(hObject,'Tag');
			if strcmp(thisProperty,'motionCorrectionRefFrame')==1
				defVal = num2str(obj.motionCorrectionRefFrame);
			else
				defVal = char(gString{1});
			end
			if strcmp(gString{gValue},userSelectStr)==1
				inputCheck = 'Americans love a winner.';
				lengthTooLongFlag = 1;
				numericCheckFlag = 1;
				nanCheckFlag = 0;
				updateValueFlag = 0;
				loopBreak = 0;
				while (nanCheckFlag==1||lengthTooLongFlag==1)&loopBreak==0
					movieSettings = inputdlg({...
							gTooltip...
						},...
						[thisProperty ' settings'],[1 70],...
						{...
							defVal...
						}...
					);

					if isempty(movieSettings)
						inputCheck = [];
						updateValueFlag = 0;
						loopBreak = 1;
						% uiwait(ciapkg.overloaded.msgbox('Re-enter a value, do not CANCEL the input dialog.'));
					else
						movieSettings = movieSettings{1};
						inputCheck = movieSettings;

						% Check users has input correct (no str or NaN), else ask again.
						if length(str2num(inputCheck))~=1
						    lengthTooLongFlag = 1;
						else
							lengthTooLongFlag = 0;
						end
						if isempty(str2num(inputCheck))
							numericCheckFlag = 0;
						else
							numericCheckFlag = 1;
						end
						if isnan(str2num(inputCheck))
							nanCheckFlag = 1;
						else
							nanCheckFlag = 0;
						end
						if nanCheckFlag==1||numericCheckFlag==0||lengthTooLongFlag==1
							if any(ismember(thisProperty,vectorAllowList))
								loopBreak = 1;
								updateValueFlag = 1;
							else
								uiwait(ciapkg.overloaded.msgbox('Please input a SINGLE numeric value (no strings or vectors).'));
							end
						else
							loopBreak = 1;
							updateValueFlag = 1;
						end
					end
				end
				% Reset to default value if user cancels
				if updateValueFlag==0
					movieSettings = defVal;
					updateValueFlag = 1;
				end
				if updateValueFlag==1
					gStringNew = [movieSettings;gString];
					set(hObject,'String',gStringNew);
					set(hObject,'Value',1);
					% preprocessingSettingsAll.(gProperty).val = [str2num(movieSettings{1}),preprocessingSettingsAll.(gProperty).val];
					preprocessingSettingsAll.(gProperty).val = {[str2num(movieSettings),preprocessingSettingsAll.(gProperty).val{1}]};
					preprocessingSettingsAll.(gProperty).str = {gStringNew};
					% turboregSettingDefaults.(thisProperty) = [str2num(movieSettings{1}),turboregSettingDefaults.(thisProperty)];
					% turboregSettingStr.(thisProperty) = gStringNew;
				end
			end
		end

		% Indicate value was modified for later input of previous settings.
		gProperty = get(hObject,'Tag');
		gValue = get(hObject,'Value');
		gString = get(hObject,'String');
		preprocessingSettingsAll.(gProperty).modified = [];
		preprocessingSettingsAll.(gProperty).modified.val = preprocessingSettingsAll.(gProperty).val{1}{gValue};
		preprocessingSettingsAll.(gProperty).modified.str = gString{gValue};
		% preprocessingSettingsAll.(gProperty).modified
		% preprocessingSettingsAll.(gProperty)
	end
	function subfxn_closeOptions(src,event)
		delete(emptyBox);
	end
	function subfxn_closeOptionsKey(src,event)
		if event.Key=="return"
			delete(emptyBox)
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
close;
clear;

load models.mat;
scale_rat = 10;
NN_Number = 5;

%[ img1_desc_types, detectedPts1, locations1,pts1 ] = create_mesh_sift_features(vertex1, faces1, scale_rat);
[ img1_desc_types, ~, ~,pts1 ] = create_mesh_spin_features(vertex1, faces1,5,scale_rat);

%[ img2_desc_types, detectedPts2, locations2,pts2] = create_mesh_sift_features(vertex2, faces2, scale_rat);
[ img2_desc_types, ~, ~,pts2] = create_mesh_spin_features(vertex2, faces2,5,scale_rat);

DescrData = pdist2(img1_desc_types,img2_desc_types);

[~,NN_Data]=sort(DescrData,2,'ascend');
NN_Data=NN_Data(:,1:NN_Number);

author: Sable

repo: mcbench-benchmarks


text:
%ANGDIFF Difference of two angles
%
% D = ANGDIFF(TH1, TH2) returns the difference between angles TH1 and TH2 on
% the circle.  The result is in the interval [-pi pi).  If TH1 is a column
% vector, and TH2 a scalar then return a column vector where TH2 is modulo
% subtracted from the corresponding elements of TH1.
%
% D = ANGDIFF(TH) returns the equivalent angle to TH in the interval [-pi pi).
%

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function d = angdiff(th1, th2)

    if nargin < 2
% THIS IS A BAD IDEA, WHERE IS IT USED?
%         if length(th1) > 1
%             d = th1(1) - th1(2);
%         else
%             d = th1;
%         end
        d = th1;
    else
        d = th1 - th2;
    end

    d = mod(d+pi, 2*pi) - pi;

% Simplistic version of the code, easy to see what it does, but slow...
%
% for very negative angles keep adding 2pi
%     while true
%         k = find(d < -pi);
%         if isempty(k)
%             break;
%         end
%         d(k) = d(k) + 2*pi;
%     end
%
%     % for very positive angles keep subtracting 2pi
%     while true
%         k = find(d > pi);
%         if isempty(k)
%             break;
%         end
%         d(k) = d(k) - 2*pi;
%     end

author: Allopart

repo: rbpf-gmapping


text:
function val = littleendian

% LITTLEENDIAN returns 1 (true) on a little endian machine, e.g. with an
% Intel or AMD, or 0 (false) otherwise
%
% Example
%   if (littleendian)
%     % do something, e.g. swap some bytes
%    end
%
% See also BIGENDIAN, SWAPBYTES, TYPECAST

% Copyrigth (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

val = (typecast(uint8([0 1]), 'uint16')==256);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [upsilonvp, upsilon] = lfmvpComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode)

% LFMVPCOMPUTEUPSILONMATRIX Upsilon matrix vel. pos. with t1, t2 limits
% FORMAT
% DESC computes a portion of the LFM kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% SEEALSO : lfmComputeUpsilonMatrix.F, lfmComputeH3.m
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

sigma = sqrt(sigma2);
gridt1 = repmat(t1, 1, length(t2));
gridt2 = repmat(t2', length(t1), 1);
timeGrid = gridt1 - gridt2;

if mode==0
    if nargout > 1
        upsilon = lfmComputeUpsilonMatrix(gamma, sigma2, t1, t2);
        upsilonvp = -gamma*upsilon + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2);
    else
        upsilonvp = -gamma*lfmComputeUpsilonMatrix(gamma, sigma2, t1, t2) ...
            + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2);
    end
else
    if nargout > 1
        upsilon = lfmComputeUpsilonMatrix(gamma, sigma2, t1, t2);
        upsilonvp = gamma*upsilon - (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2) ...
            + (2/(sqrt(pi)*sigma))*exp(-gamma*t1)*(exp(-(t2.^2)/sigma2)).';
    else
        upsilonvp = gamma*lfmComputeUpsilonMatrix(gamma, sigma2, t1, t2) ...
            - (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2) ...
            + (2/(sqrt(pi)*sigma))*exp(-gamma*t1)*(exp(-(t2.^2)/sigma2)).';
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [H]=freqresp(b,a,w)
% Frequency response function from difference equation
% [H] = freqresp(b,a,w)
% H = frequency response array evaluated a w frequencies
% b = numerator coefficient array
% a = denominator coefficient array (a(1)=1)
% w = frequency location array
m = 0:length(b)-1;
l = 0:length(a)-1;
num = b*exp(-j*m'*w);
den = a*exp(-j*l'*w);
H = num./den;


text:
function output = vbmc_output(vp,optimState,msg,stats,idx_best,vbmc_version)
%VBMC_OUTPUT Create OUTPUT struct for VBMC.

output.function = func2str(optimState.fun);
if all(isinf(optimState.LB)) && all(isinf(optimState.UB))
    output.problemtype = 'unconstrained';
else
    output.problemtype = 'boundconstraints';
end
output.iterations = optimState.iter;
output.funccount = optimState.funccount;
output.bestiter = idx_best;
output.trainsetsize = stats.Neff(idx_best);
output.components = vp.K;
output.rindex = stats.rindex(idx_best);
if stats.stable(idx_best)
    output.convergencestatus = 'probable';
else
    output.convergencestatus = 'no';
end
output.overhead = NaN;
output.rngstate = rng;
output.algorithm = 'Variational Bayesian Monte Carlo';
output.version = vbmc_version;
output.message = msg;

output.elbo = vp.stats.elbo;
output.elbo_sd = vp.stats.elbo_sd;

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function s=sprintsi(x,d,w)
%SPRINTSI Print X with SI multiplier S=(X,D,W)
% D is number of decimal places (+ve) or significant digits (-ve) [default=-3]
% |W| is total width including multiplier
% if W<=0 then trailing 0's will be eliminated
%
% Example: sprintsi(2345,-2) gives '2.3 k'

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

if nargin<3 w=0; end;
if nargin<2 d=-3; end;
f='afpnum kMGT';
e=max(-18,min(12,floor(log10(abs(x)))));
k=floor(e/3);
dp=max([0 d 3*k-d-e-1]);
if w<=0 & dp
   w=abs(w);
   dp=max(find([1 mod(mod(round(x*10^(dp-3*k)),10^dp),10.^(dp:-1:1))]))-1;
end
if(k)
   s=sprintf(sprintf('%%%d.%df %c',max(w-2,0),dp,f(k+7)),x*1e-3^k);
else
   s=sprintf(sprintf('%%%d.%df ',max(w-1,0),dp),x*1e-3^k);
end

author: jtkim-kaist

repo: VAD

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


text:
function ll = robThreeDynamicsLogLikelihood(model)

% ROBTHREEDYNAMICSLOGLIKELIHOOD Give the log likelihood of the robot three dynamics part.

% FGPLVM

ll = 0;
for i = 1:size(model.diffX)-1
  covMat = model.lambda*model.diffX(i, :)'*model.diffX(i, :) + ...
      eye(2)*model.sigma2;
  invCovMat = inv(covMat);
  ll = ll -0.5* model.diffX(i+1, :)*invCovMat*model.diffX(i+1, :)';
end


text:
function doZeroRange(solver,varargin)

zrm = getClass(varargin,'zeroRangeMethod');

if isempty(zrm), zrm = zeroRangeMethod(solver.pf,varargin{:}); end

isZero = zrm.checkZeroRange(solver.S3G);

solver.S3G = subGrid(solver.S3G,~isZero);

end

author: mtex-toolbox

repo: mtex

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


text:
function [g,a,info] = dtwfb2filterbank( dualwt, varargin)
%DTWFB2FILTERBANK DTWFB equivalent non-iterated filterbank
%   Usage: [g,a] = dtwfb2filterbank(dualwt)
%          [g,a,info] = dtwfb2filterbank(...)
%
%   Input parameters:
%      dualwt : Dual-tree wavelet filterbank specification.
%
%   Output parameters:
%      g      : Cell array of filters.
%      a      : Downsampling rate for each channel.
%      info   : Additional information.
%
%   `[g,a] = dtwfb2filterbank(dualwt)` constructs a set of filters *g* and
%   subsampling factors *a* of a non-iterated filterbank, which is
%   equivalent to the dual-tree wavelet filterbank defined by *dualwt*.
%   The returned parameters can be used directly in |filterbank| and other
%   routines. The format of *dualwt* is the same as in |dtwfb| and
%   |dtwfbreal|.
%
%   The function internally calls |dtwfbinit| and passes *dualwt* and all
%   additional parameters to it.
%
%   `[g,a,info] = dtwfb2filterbank(...)` additionally outputs a *info*
%   struct containing equivalent filterbanks of individual real-valued
%   trees as fields *info.g1* and *info.g2*.
%
%   Additional parameters:
%   ----------------------
%
%   `'real'`
%      By default, the function returns a filtebank equivalent to |dtwfb|.
%      The filters can be restricted to cover only the positive frequencies
%      and to be equivivalent to |dtwfbreal| by passing a `'real'` flag.
%
%   `'freq'`(default),`'nat'`
%     The filters are ordered to produce subbands in the same order as
%     |dtwfb| or |dtwfbreal| with the same flag.
%
%   Examples:
%   ---------
%
%   The following two examples create a multirate identity filterbank
%   using a duel-tree of depth 3:::
%
%      [g,a] = dtwfb2filterbank({'qshift3',3},'real');
%      filterbankfreqz(g,a,1024,'plot','linabs');
%
%   In the second example, the filterbank is identical to the full
%   wavelet tree:::
%
%      [g,a] = dtwfb2filterbank({'qshift3',3,'full'},'real');
%      filterbankfreqz(g,a,1024,'plot','linabs');
%
%   See also: dtwfbinit

complainif_notenoughargs(nargin,1,'DTWFB2FILTERBANK');

% Search for the 'real' flag
do_real = ~isempty(varargin(strcmp('real',varargin)));
if do_real
    %Remove the 'real' flag from varargin
    varargin(strcmp('real',varargin)) = [];
end

if ~isempty(varargin(strcmp('complex',varargin)))
    %Remove the 'complex' flag from varargin
    %It is not used elsewhere anyway
    varargin(strcmp('complex',varargin)) = [];
end

% Initialize the dual-tree
dtw = dtwfbinit({'strict',dualwt},varargin{:});

% Determine relation between the tree nodes
[wtPath, rangeLoc, rangeOut] = treeBFranges(dtw);
slice = ~cellfun(@isempty,rangeOut); % Limit to nodes with unconnected outputs
wtPath   = wtPath(slice);
rangeLoc = rangeLoc(slice);
rangeOut = rangeOut(slice);

% Multirate identity filters of the first tree
[g1,a] = nodesMultid(wtPath,rangeLoc,rangeOut,dtw);

% Multirate identity filters of the second tree
dtw.nodes = dtw.dualnodes;
g2 = nodesMultid(wtPath,rangeLoc,rangeOut,dtw);

if nargin>2
   % Return the filterbanks before doing the alignment
   info.g1 = cellfun(@(gEl) setfield(gEl,'h',gEl.h/2),g1,'UniformOutput',0);
   info.g2 = cellfun(@(gEl) setfield(gEl,'h',gEl.h/2),g2,'UniformOutput',0);
end

% Align filter offsets so they can be summed
for ii = 1:numel(g1)
   % Sanity checks
   assert(g1{ii}.offset<=0,sprintf('%s: Invalid wavelet filters.',upper(mfilename)));
   assert(g2{ii}.offset<=0,sprintf('%s: Invalid wavelet filters.',upper(mfilename)));

   % Insert zeros and update offsets
   offdiff = g1{ii}.offset-g2{ii}.offset;
   if offdiff>0
       g1{ii}.offset = g1{ii}.offset - offdiff;
       g1{ii}.h = [zeros(offdiff,1);g1{ii}.h(:)];
   elseif offdiff<0
       g2{ii}.offset = g2{ii}.offset + offdiff;
       g2{ii}.h = [zeros(-offdiff,1);g2{ii}.h(:)];
   end

   % Pad with zeros to a common length
   lendiff = numel(g1{ii}.h) - numel(g2{ii}.h);
   if lendiff~=0
       maxLen = max(numel(g1{ii}.h),numel(g2{ii}.h));
       g1{ii}.h = postpad(g1{ii}.h,maxLen);
       g2{ii}.h = postpad(g2{ii}.h,maxLen);
   end
end

% Filters covering the positive frequencies
g = cellfun(@(gEl,g2El) setfield(gEl,'h',(gEl.h+1i*g2El.h)/2),g1,g2,'UniformOutput',0);

% Mirror the filters when negative frequency filters are required too
if ~do_real
   gneg = cellfun(@(gEl,g2El) setfield(gEl,'h',(gEl.h-1i*g2El.h)/2),g1,g2,'UniformOutput',0);
   g = [g;gneg(end:-1:1)];
   a = [a;a(end:-1:1)];
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function out = mtimes( a,b )
%MTIMES for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% determine input types
inType = {class(a), class(b)};
if(strcmp(inType{1},'TRAFO') && strcmp(inType{2},'TRAFO'))
    [outA, idx] = flattenCellMatrix(a.data);
    [outB, idxB] = flattenCellMatrix(b.data);
    meta = a.meta;
    meta_b = b.meta;
    clear 'a' 'b'

    if(~isempty(idx) && ~isempty(idxB) && ~all(cellfun(@(x,y) isequal(x,y), idx, idxB)))
        error('TRAFO::mtimes: Nested cells must have the same size');
    end
    if(iscell(meta) && iscell(meta_b))
        if(~isempty(meta) && ~isempty(meta_b))
            if(iscell(meta{1}) && iscell(meta_b{1}))
                if(~all(cellfun(@(x,y) isequal(x,y), flattenCellMatrix(meta), flattenCellMatrix(meta_b))))
                    error('TRAFO::mtimes: Unequal reconstruction information');
                end
            else
                if(~all(cellfun(@(x,y) isequal(x,y), meta, meta_b)))
                    error('TRAFO::mtimes: Unequal reconstruction information');
                end
            end
        end
    else
        if ~isequal(meta, meta_b)
            error('TRAFO::mtimes: Unequal reconstruction information');
        end
    end

    out = cellfun(@(x,y) x*y, outA, outB, 'UniformOutput', false);

elseif(strcmp(inType{1},'TRAFO')) % assume b is double/uint/int
    [out, idx] = flattenCellMatrix(a.data);
    meta = a.meta;
    clear 'a'

    out = cellfun(@(x) x * b, out, 'UniformOutput', false);

elseif(strcmp(inType{2},'TRAFO'))
    [out, idx] = flattenCellMatrix(b.data);
    meta = b.meta;
    clear 'b'

    out = cellfun(@(x) a * x, out, 'UniformOutput', false);

else
    error('TRAFO::mtimes: Impossible constellation!');
end

out = reconFlatCellMatrix(out,idx);
out = TRAFO(out,meta); % return TRAFO object again without modifying input

end

author: thomaskuestner


text:
function relm_PaintSet(vRelm_LTest, vRelm_RTest)

figure

nLen = length(vRelm_LTest);
for nCnt = 1:nLen
  H = subplot(nLen, 2, (nCnt*2)-1);
  relm_PaintCumPlot(vRelm_LTest(nCnt), 'Log-Likelihood', H);
  H = subplot(nLen, 2, (nCnt*2));
  relm_PaintCumPlot(vRelm_RTest(nCnt), 'Log-Likelihood Ratio', H);
end

author: CelsoReyes


text:
function chebyshev_polynomial_test11 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST11 tests U_POLYNOMIAL_COEFFICIENTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST11\n' );
  fprintf ( 1, '  U_POLYNOMIAL_COEFFICIENTS determines the\n' );
  fprintf ( 1, '  polynomial coefficients for U(n,x).\n' );

  c = u_polynomial_coefficients ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(%d,x) = \n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( c(i+1,j+1) == 0.0 )

      elseif ( j == 0 )
        fprintf ( 1, '  %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function kern = rbfinfwhiteKernParamInit(kern)

% RBFINFWHITEKERNPARAMINIT The RBF-WHITE-INF kernel is a convolutional
% kernel obtained as the result of convolving a white noise input process
% with an RBF smoothing kernel with range between minus and plus infinity.
%
% Although it be used independently, it is mostly intended to be combined
% with other kernels in a compound kernel. For example, this is the kernel
% used as variational smoothing kernel for the DTC sparse GP approach.
%
% The parameters are sigma2, the process variance (kern.variance), and
% gamma, the inverse width (kern.inverseWidth). The inverse width controls
% how wide the basis functions are, the larger gamma, the smaller the basis
% functions are.
%
% It is very similar to the RBF-WHITE kernel.
%
% SEEALSO : cmpndKernParamInit, rbfwhiteKernParamInit
%
% FORMAT
% DESC initialises the RBF-WHITE kernel structure with some default
% parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : David Luengo, 2009

% KERN

if kern.inputDimension > 1
  error('RBF-WHITE kernel only valid for one-D input.')
end

kern.nParams = 2;
kern.inverseWidth = 1;
kern.variance = 1;

% Constrains parameters to be positive for optimisation.
kern.transforms.index = [1 2];
kern.transforms.type = optimiDefaultConstraint('positive');
kern.isStationary = true;


text:
%% Elasticity
%
%%


text:
function nifti = io_CSIload_nifti2(filename)
    arguments
        filename (1, :) char {mustBeFile}
    end
    matlabError = [];
    nifti = struct();
    try
        %open the file
        fileID  = fopen(filename, 'r');

        %sizeof_hdr: Needs to be 540
        header = fread(fileID, 540, '*int8')';

        nifti = readHeader(nifti, header);
        nifti = readExtensions(nifti, fileID);

        if(contains(nifti.datatype, 'complex', 'IgnoreCase', true))
            if(contains(nifti.datatype, '64'))
                dataType = 'float';
            elseif(contains(nifti.datatype, '128'))
                dataType = 'double';
            else
                fprintf('Datatype is %s bits long! Matlab does not contain this datatype. Using chars', bits{1});
                dataType = 'char';
            end
            data = fread(fileID, Inf, dataType);
            complexData = reshape(data, 2, []);
            data = complex(complexData(1, :), complexData(2, :));
        else
            data = fread(fileID, Inf, nifti.datatype);
        end
        data = reshape(data, nifti.dim(2:end));
        nifti.data = data;
    catch matlabError
        % catch here to catch errors in io. Then close file.
    end

    %close file
    fclose(fileID);

    % throw matlab error
    if ~isempty(matlabError)
        rethrow(matlabError);
    end
end

function dataType = getDataTypeFromInt(dataTypeInt)

    intTypes = [2, 4, 8, 16, 32, 64, 128, 256, 512, 768, 1024, 1280, 1536, 1792, 2048];
    dataLabels = ["UINT8", "INT16", "INT32", ...
        "FLOAT32", "COMPLEX64", "FLOAT64", ...
        "RGB24", "INT8", "UINT16", "UINT32", ...
        "INT64", "UINT64", "FLOAT128", ...
        "COMPLEX128", "COMPLEX256"];
    map = containers.Map(intTypes, dataLabels);
    dataType = map(dataTypeInt);
end

function nifti = readHeader(nifti, header)
    nifti.sizeof_hdr = typecast(header(1:4), 'int32');
    nifti.magic = cast(header(5:12), 'char');
    nifti.datatype =  getDataTypeFromInt(typecast(header(13:14)', 'int16'));
    nifti.bitpix =  typecast(header(15:16), 'int16');
    nifti.dim =  typecast(header(17:80), 'int64');
    nifti.intent_p1 =  typecast(header(81:88), 'double');
    nifti.intent_p2 =  typecast(header(89 :96 ), 'double');
    nifti.intent_p3 =  typecast(header(97 :104 ), 'double');
    nifti.pixdim = typecast(header(105 :168 ), 'double');
    nifti.vox_offset = typecast(header(169 :176 ), 'int64');
    nifti.scl_slope = typecast(header(177 :184 ), 'double');
    nifti.scl_inter = typecast(header(185 :192 ), 'double');
    nifti.cal_max = typecast(header(193 :200 ), 'double');
    nifti.cal_min = typecast(header(201 :208 ), 'double');
    nifti.slice_duration = typecast(header(209: 216), 'double');
    nifti.tofset = typecast(header(217 :224 ), 'double');
    nifti.slice_start = typecast(header(225 :232 ), 'int64');
    nifti.slice_end = typecast(header(233 :240 ), 'int64');
    nifti.descrip = cast(header(241 :320 ), 'char');
    nifti.aux_file = cast(header(321 :344 ), 'char');
    nifti.qform_cod = typecast(header(345 :348 ), 'int32');
    nifti.sform_cod = typecast(header(349 :352 ), 'int32');
    nifti.quatern_b = typecast(header(353 :360 ), 'double');
    nifti.quatern_c = typecast(header(361 :368 ), 'double');
    nifti.quatern_d = typecast(header(369 :376 ), 'double');
    nifti.qoffset_x = typecast(header(377 :384 ), 'double');
    nifti.qoffset_y = typecast(header(385 :392 ), 'double');
    nifti.qoffset_z = typecast(header(393 :400 ), 'double');
    nifti.srow_x = typecast(header(401 :432 ), 'double');
    nifti.srow_y = typecast(header(433 :464 ), 'double');
    nifti.srow_z = typecast(header(465 :496 ), 'double');
    nifti.slice_code = typecast(header(497 :500 ), 'int32');
    nifti.xyzt_units = typecast(header(501 :504 ), 'int32');
    nifti.intent_cod3 = typecast(header(505 :508 ), 'int32');
    nifti.intent_name = cast(header(509: 524), 'char');
    nifti.dim_info = cast(header(525 : 525), 'char');
end

function nifti = readExtensions(nifti, fileID)
    extension = fread(fileID, 4, 'int8=>int8');
    extensionLength = fread(fileID, 1, 'int32');
    nifti.extensionCode = fread(fileID, 1, 'int32');
    extensionJson = fread(fileID, extensionLength - 8, '*char')';
    nifti.extension = jsondecode(deblank(extensionJson));
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [Au, PHIu, Av, PHIv, w, TWOCIR]=ep2ap(SEMA, ECC, INC, PHA, plot_demo)
%
% Convert tidal ellipse parameters into amplitude and phase lag parameters.
% Its inverse is app2ep.m. Please refer to app2ep for the meaning of the
% inputs and outputs.
%
% Zhigang Xu
% Oct. 20, 2000
%
% Document:  tidal_ellipse.ps
%
if nargin < 5
     plot_demo=0;  % by default, no plot for the ellipse
end

   Wp = (1+ECC)/2 .*SEMA;
   Wm = (1-ECC)/2 .*SEMA;
   THETAp = INC-PHA;
   THETAm = INC+PHA;

   %convert degrees into radians
   THETAp = THETAp/180*pi;
   THETAm = THETAm/180*pi;

   %Calculate wp and wm.
   wp = Wp.*exp(i*THETAp);
   wm = Wm.*exp(i*THETAm);

   if nargout >= 5
      ndot=36;
      dot = 2*pi/ndot;
      ot = [0:dot:2*pi-dot];
      w = wp(:)*exp(i*ot)+wm(:)*exp(-i*ot);
      w=reshape(w, [size(wp) ndot]);
   end

   % Calculate cAu, cAv --- complex amplitude of u and v
   cAu = wp+conj(wm);
   cAv = -i*(wp-conj(wm));
   Au  = abs(cAu);
   Av  = abs(cAv);
   PHIu = -angle(cAu)*180/pi;
   PHIv = -angle(cAv)*180/pi;

   % flip angles in the range of [-180 0) to the range of [180 360).
   id = PHIu < 0; PHIu(id) = PHIu(id) + 360;
   id = PHIv < 0; PHIv(id) = PHIv(id) + 360;

  if any(plot_demo)
     plot_ell(SEMA,ECC,INC,PHA,plot_demo);
  end

  if nargout == 6
      TWOCIR=struct('Wp', Wp, 'THETAp', THETAp, 'wp', ...
            wp, 'Wm', Wm, 'THETAm', THETAm, 'wm', wm, 'ot', ot, 'dot', dot);
 end

%Authorship Copyright:
%
%    The author of this program retains the copyright of this program, while
% you are welcome to use and distribute this program as long as you credit
% the author properly and respect the program name itself. Particularly,
% you are expected to retain the original author's name in this original
% version of the program or any of its modified version that you might make.
% You are also expected not to essentially change the name of the programs
% except for adding possible extension for your own version you might create,
% e.g. app2ep_xx is acceptable.  Any suggestions are welcome and enjoy my
% program(s)!
%
%
%Author Info:
%_______________________________________________________________________
%  Zhigang Xu, Ph.D.
%  (pronounced as Tsi Gahng Hsu)
%  Research Scientist
%  Coastal Circulation
%  Bedford Institute of Oceanography
%  1 Challenge Dr.
%  P.O. Box 1006                    Phone  (902) 426-2307 (o)
%  Dartmouth, Nova Scotia           Fax    (902) 426-7827
%  CANADA B2Y 4A2                   email xuz@dfo-mpo.gc.ca
%_______________________________________________________________________
%
%Release Date: Nov. 2000


text:
function c=dctii(f,L,dim)
%DCTII  Discrete Consine Transform type II
%   Usage:  c=dctii(f);
%           c=dctii(f,L);
%           c=dctii(f,[],dim);
%           c=dctii(f,L,dim);
%
%   `dctii(f)` computes the discrete cosine transform of type II of the
%   input signal *f*. If *f* is multi-dimensional, the transformation is
%   applied along the first non-singleton dimension.
%
%   `dctii(f,L)` zero-pads or truncates *f* to length *L* before doing the
%   transformation.
%
%   `dctii(f,[],dim)` or `dctii(f,L,dim)` applies the transformation along
%   dimension *dim*.
%
%   The transform is real (output is real if input is real) and orthonormal.
%
%   This is the inverse of |dctiii|.
%
%   Let *f* be a signal of length *L*, let `c=dctii(f)` and define the
%   vector *w* of length *L* by
%
%   .. w = [1/sqrt(2) 1 1 1 1 ...]
%
%   .. math:: w\left(n\right)=\begin{cases}\frac{1}{\sqrt{2}} & \text{if }n=0\\1 & \text{otherwise}\end{cases}
%
%   Then
%
%   ..                     L-1
%     c(n+1) = sqrt(2/L) * sum w(n+1)*f(m+1)*cos(pi*n*(m+.5)/L)
%                          m=0
%
%   .. math:: c\left(n+1\right)=\sqrt{\frac{2}{L}}\sum_{m=0}^{L-1}w\left(n\right)f\left(m+1\right)\cos\left(\frac{\pi}{L} n\left(m+\frac{1}{2}\right)\right)
%
%   Examples:
%   ---------
%
%   The following figures show the first 4 basis functions of the DCTII of
%   length 20:::
%
%     % The dctiii is the adjoint of dctii.
%     F=dctiii(eye(20));
%
%     for ii=1:4
%       subplot(4,1,ii);
%       stem(F(:,ii));
%     end;
%
%   See also:  dctiii, dctiv, dstii
%
%   References: rayi90 wi94

%   AUTHOR: Peter L. Søndergaard
%   TESTING: TEST_PUREFREQ
%   REFERENCE: REF_DCTII

complainif_argnonotinrange(nargin,1,3,mfilename);

if nargin<3
  dim=[];
end;

if nargin<2
  L=[];
end;

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'DCTII');

if ~isempty(L)
  f=postpad(f,L);
end;
  c=comp_dct(f,2);
% c=zeros(L,W,assert_classname(f));
%
% m1=1/sqrt(2)*exp(-(0:L-1)*pi*i/(2*L)).';
% m1(1)=1;
%
% m2=1/sqrt(2)*exp((1:L-1)*pi*i/(2*L)).';
%
% s1=fft([f;flipud(f)]);
%
% % This could be done by a repmat instead.
% for w=1:W
%   c(:,w)=s1(1:L,w).*m1+[0;s1(2*L:-1:L+2,w).*m2];
% end;
%
% c=c/sqrt(L)/2;

if isreal(f)
  c=real(c);
end;

c=assert_sigreshape_post(c,dim,permutedsize,order);

% This is a slow, but convenient way of expressing the algorithm.
%R=1/sqrt(2)*[diag(exp((0:L-1)*pi*i/(2*L)));...
%	     zeros(1,L); ...
%	     [zeros(L-1,1),flipud(diag(exp(-(1:L-1)*pi*i/(2*L))))]];

%R(1,1)=1;

%c=R'*fft([f;flipud(f)])/sqrt(L)/2;


text:
function [M2,NM2]=compute_Hist_Metric2(img1,img2,index)

img1 = double(img1)/255;
img2 = double(img2)/255;
img1 = img1 - 0.01*(img1==1);
img2 = img2 - 0.01*(img2==1);

N = length(index);
sN = 20000;
gray_level = 64;%round(N^(1/(2+alpha)));
s1 = img1(index);
s2 = img2(index);
rand('state',1);
% if sampling
%     sindex = round(rand(sN,1)*(N-1))+1;
%     s1 = s1(sindex);
%     s2 = s2(sindex);
%     gray_level = round(sN^(1/(2+alpha)));
% end

bitLen = floor(log2(gray_level)-0.01) + 1;

sd1 = floor(s1*gray_level);
sd2 = floor(s2*gray_level);

h1 = getHist(sd1,bitLen);
h2 = getHist(sd2,bitLen);

h12 = getHist([sd1 sd2],bitLen);

e1 = -sum(h1.*log(h1));
e2 = -sum(h2.*log(h2));

e12 = -sum(h12.*log(h12));

M2 = e12 - e1- e2;
NM2 = M2/e12;

author: zhouyuanzxcv

repo: Hyperspectral


text:
function x = tridisolve ( a, b, c, d )

%*****************************************************************************80
%
%% TRIDISOLVE solves a tridiagonal system of linear equations.
%
%  Discussion:
%
%    We can describe an NxN tridiagonal matrix by vectors A, B, and C, where
%    A and C are of length N-1.  In that case, a linear system can be
%    represented as
%                        b(1) * x(1) + c(1) * x(2)   = d(1),
%      a(j-1) * x(j-1) + b(j) * x(j) + c(j) * x(j+1) = d(j), j = 2:n-1,
%      a(n-1) * x(n-1) + b(n) * x(n)                 = d(n)
%
%    This function produces the solution vector X.
%
%    This function is derived from Cleve Moler's Matlab suite.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Cleve Moler.
%
%  Parameters:
%
%    Input, real A(N-1), B(N), C(N-1), the matrix entries.
%
%    Input, real D(N), the right hand side.
%
%    Output, real X(N), the solution.
%
  x = d;
  n = length ( x );
  bi = zeros ( n, 1 );

  for j = 1 : n - 1
    bi(j) = 1.0 / b(j);
    mu = a(j) * bi(j);
    b(j+1) = b(j+1) - mu * c(j);
    x(j+1) = x(j+1) - mu * x(j);
  end

  x(n) = x(n) / b(n);
  for j = n - 1 : -1 : 1
    x(j) = ( x(j) - c(j) * x(j+1) ) * bi(j);
  end

  return
end


text:
function truncated_normal_ab_pdf_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_PDF_TEST tests TRUNCATED_NORMAL_AB_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_AB_PDF_TEST\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_AB_PDF evaluates the Truncated Normal PDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean = mu\n' );
  fprintf ( 1, '    standard deviation = sigma\n' );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval [a,b]\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                                                           Stored         Computed\n' );
  fprintf ( 1, '       X        Mu         S         A         B             PDF             PDF\n' );
  fprintf ( 1, '\n');

  n_data = 0;

  while ( 1 )

    [ n_data, mu, sigma, a, b, x, pdf1 ] ...
      = truncated_normal_ab_pdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    pdf2 = truncated_normal_ab_pdf ( x, mu, sigma, a, b );

    fprintf( 1, '  %8.1f  %8.1f  %8.1f  %8.1f  %8.1f  %14g  %14g\n', x, mu, sigma, a, b, pdf1, pdf2 );

  end

  return
end


text:
function value = r8_cin ( x )

%*****************************************************************************80
%
%% R8_CIN evaluates the alternate cosine integral Cin of an R8 argument.
%
%  Discussion:
%
%    CIN(X) = gamma + log(X)
%      + integral ( 0 <= T <= X ) ( cos ( T ) - 1 ) / T  dT
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the cosine integral Cin evaluated at X.
%
  persistent cincs
  persistent eul
  persistent ncin
  persistent xmin

  eul = 0.57721566490153286060651209008240;

  if ( isempty ( ncin ) )

    cincs = [ ...
      0.37074501750909688741654801228564992, ...
     -0.05893574896364446831956864397363697, ...
      0.00538189642113569124048745326203340, ...
     -0.00029860052841962135319594906563410, ...
      0.00001095572575321620077031054467306, ...
     -0.00000028405454877346630491727187731, ...
      0.00000000546973994875384912457861806, ...
     -0.00000000008124187461318157083277452, ...
      0.00000000000095868593117706609013181, ...
     -0.00000000000000920266004392351031377, ...
      0.00000000000000007325887999017895024, ...
     -0.00000000000000000049143726675842909, ...
      0.00000000000000000000281577746753902, ...
     -0.00000000000000000000001393986788501, ...
      0.00000000000000000000000006022485646, ...
     -0.00000000000000000000000000022904717, ...
      0.00000000000000000000000000000077273, ...
     -0.00000000000000000000000000000000233 ]';

    ncin = r8_inits ( cincs, 18, 0.1 * r8_mach ( 3 ) );
    xmin = sqrt ( r8_mach ( 1 ) );

  end

  absx = abs ( x );

  if ( absx <= xmin )
    value = 0.0;
  elseif ( absx <= 4.0 )
    value = r8_csevl ( ( x * x - 8.0 ) * 0.125, cincs, ncin ) * x * x;
  else
    [ f, g ] = r8_sifg ( absx );
    sinx = sin ( absx );
    value = - f * sinx + g * cos ( absx ) + log ( absx ) + eul;
  end

  return
end


text:
function [hog_data, valid_inds, vid_id] = Read_HOG_files_small(hog_files, hog_data_dir, num_samples)

    hog_data = [];
    vid_id = {};

    feats_filled = 0;

    curr_data_buff = [];

    for i=1:numel(hog_files)

        hog_file = [hog_data_dir, hog_files(i).name];

        fprintf('%d %s\n', i, hog_file);

        f = fopen(hog_file, 'r');

        curr_ind = 0;

        while(~feof(f))

            if(i == 1 && curr_ind == 0)
                num_cols = fread(f, 1, 'int32');
                if(isempty(num_cols))
                    break;
                end

                num_rows = fread(f, 1, 'int32');
                num_chan = fread(f, 1, 'int32');

                curr_ind = curr_ind + 1;

                % preallocate some space
                if(curr_ind == 1)
                    curr_data_buff = zeros(5000, 1 + num_rows * num_cols * num_chan);

                    num_feats =  1 + num_rows * num_cols * num_chan;
                end

                if(curr_ind > size(curr_data_buff,1))
                    curr_data_buff = cat(1, curr_data_buff, zeros(6000, num_rows * num_cols * num_chan));
                end
                feature_vec = fread(f, [1, 1 + num_rows * num_cols * num_chan], 'float32');
                curr_data_buff(curr_ind, :) = feature_vec;

            else

                % Reading in batches of 5000

                feature_vec = fread(f, [4 + num_rows * num_cols * num_chan, 5000], 'float32');
                feature_vec = feature_vec(4:end,:)';

                if(~isempty(feature_vec))
                    num_rows_read = size(feature_vec,1);

                    curr_data_buff(curr_ind+1:curr_ind+num_rows_read,:) = feature_vec;
                    %valid_data_buff =
                    curr_ind = curr_ind + size(feature_vec,1);
                end
            end

        end

        fclose(f);

        curr_data_small = curr_data_buff(1:curr_ind,:);
        vid_id_curr = cell(curr_ind,1);
        vid_id_curr(:) = {hog_files(i).name};

        % Keep up to 20 frames from the whole video (so that it is balanced
        % per dataset/video/participant)

        if(nargin > 2)
            num_instances = num_samples;
        else
            num_instances = 20;
        end

        increment = round(curr_ind / num_instances);
        if(increment == 0)
            increment = 1;
        end

        curr_data_small = curr_data_small(1:increment:end,:);
        vid_id_curr = vid_id_curr(1:increment:end,:);

        vid_id = cat(1, vid_id, vid_id_curr);

        % Assume same number of frames per video
        if(i==1)
            hog_data = zeros(10*numel(hog_files), num_feats);
        end

        if(size(hog_data,1) < feats_filled+size(curr_data_small,1))
           hog_data = cat(1, hog_data, zeros(feats_filled + size(curr_data_small,1) - size(hog_data,1), num_feats));
        end

        hog_data(feats_filled+1:feats_filled + size(curr_data_small,1),:) = curr_data_small;

        feats_filled = feats_filled + size(curr_data_small,1);

    end

    valid_inds = hog_data(1:feats_filled,1) > 0;
    hog_data = hog_data(1:feats_filled,2:end);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
% fig_error_J.m
% plot maximum normalized error vs J using QR method
% for min-max NUFFT with uniform scaling factors
% This is Fig. 3 in the 2003 IEEE T-SP paper.
% Jeff Fessler, The University of Michigan

if ~isvar('emax')
%	Nlist = 1;
	Nlist = 2^10;
	Mlist = [1.5 2 2.5 3 4 5];	% over-sampling factors
%	Jlist = [2:11];
	Jlist = [2:20];	% neighborhood size

	[JJ, MM, NN] = ndgrid(Jlist, Mlist, Nlist);

	emax = zeros(size(NN));

	alpha = [1]; beta = [];
	for ii=1:numel(NN)
		N = NN(ii);
		J = JJ(ii);
		M = MM(ii);
		printf('M=%d J=%d', M, J)
		K = M * N;
		gam = 2*pi/K;
		om = gam * [0:20]'/40;
		emax(ii) = max(nufft1_err_mm(om, N, J, K, 'qr', alpha, beta));
	end
end

%
% plot min-max errors
%
if 1
	semilogy(Jlist, emax(:,1), '-o', ...
		Jlist, emax(:,2), '-d', ...
		Jlist, emax(:,3), '-^', ...
		Jlist, emax(:,4), '-s', ...
		Jlist, emax(:,5), '-p', ...
		Jlist, emax(:,6), '-+')
%		Jlist(Jlist <= 99), emax(Jlist <= 99, 2), '-d', ...
%		Jlist(Jlist <= 98), emax(Jlist <= 98, 3), '-^', ...
%		Jlist(Jlist <= 96), emax(Jlist <= 96, 4), '-s', ...
%		Jlist(Jlist <= 95), emax(Jlist <= 95, 5), '-p')
%	axis([minmax(Jlist)'+[-0.2 0.2] 1e-5 1e-1])
%	axis tight
%	axisx([minmax(Jlist)'+[-0.2 0.2]])
	axis([minmax(Jlist)'+[-0.2 0.2] 1e-10 1e-0])
	xtick([2:3:20])
	ytick(10.^[-10:2:0])
%	text(2.2, 2e-4, '(for unit-norm signal)')
	xlabel J, ylabel 'E_{max}'
	title(sprintf('Maximum error for \\alpha = (%g)', alpha))
	leg = {};
	for ii=1:length(Mlist)
		leg{ii} = sprintf('K/N=%g', Mlist(ii));
	end
	legend(leg)

%	ir_savefig c 'fig_error_J'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function graph_demo (n)
%GRAPH_DEMO graph partitioning demo
%   graph_demo(n) constructs an set of n-by-n 2D grids, partitions them, and
%   plots them in one-second intervals.  n is optional; it defaults to 60.
%
%   Example:
%       graph_demo
%
%   See also DELSQ, NUMGRID, GPLOT, TREEPLOT

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 1)
    % construct a 60-by-60 grid
    n = 60 ;
end

figure (1)
clf

for regions = {'Square', 'C' 'Disc', 'Annulus', 'Heart', 'Butterfly', 'L'}

    % construct the grid
    region = regions {1} ;
    g = numgrid (region (1), n) ;
    x = repmat (0:n-1, n, 1) ;
    y = repmat (((n-1):-1:0)', 1, n)  ;
    A = delsq (g) ;
    x = x (find (g)) ;							    %#ok
    y = y (find (g)) ;							    %#ok

    % plot the original grid
    clf
    subplot (2,2,1)
    my_gplot (A, x, y)
    title (sprintf ('%s-shaped 2D grid', region)) ;
    axis equal
    axis off

    % bisect the graph
    s = bisect (A) ;
    [i j] = find (A) ;
    subplot (2,2,2)
    my_gplot (sparse (i, j, s(i) == s(j)), x, y) ;
    title ('node bisection') ;
    axis equal
    axis off

    % nested dissection
    nsmall = floor (size (A,1) / 2) ;
    defaults = 0 ;
    while (1)
        if (defaults)
            % use defaults
            [p cp cmember] = nesdis (A) ;
        else
            [p cp cmember] = nesdis (A, 'sym', nsmall) ;
        end

        % plot the components
        subplot (2,2,3)
        my_gplot (sparse (i, j, cmember(i) == cmember (j)), x, y) ;
        if (defaults)
            title ('nested dissection (defaults)') ;
        else
            title (sprintf ('nested dissection, nsmall %d', nsmall)) ;
        end
        axis equal
        axis off

        % plot the separator tree
        subplot (2,2,4)
        treeplot (cp, 'ko')
        title ('separator tree') ;
        axis equal
        axis off

        drawnow
        pause (0.1)

        if (defaults)
            break ;
        end

        nsmall = floor (nsmall / 2) ;
        if (nsmall < 20)
            defaults = 1 ;
            pause (0.2)
        end
    end
end

%-------------------------------------------------------------------------------

function my_gplot (A, x, y)
% my_gplot : like gplot, just a lot faster
[i, j] = find (A) ;
[ignore, p] = sort (max(i, j)) ;
i = i (p) ;
j = j (p) ;
nans = repmat (NaN, size (i)) ;
x = [ x(i) x(j) nans ]' ;
y = [ y(i) y(j) nans ]' ;
plot (x (:), y (:)) ;


text:
function w = rbfpak(net)
%RBFPAK	Combines all the parameters in an RBF network into one weights vector.
%
%	Description
%	W = RBFPAK(NET) takes a network data structure NET and combines the
%	component parameter matrices into a single row vector W.
%
%	See also
%	RBFUNPAK, RBF
%

%	Copyright (c) Ian T Nabney (1996-2001)

errstring = consist(net, 'rbf');
if ~errstring
  error(errstring);
end

w = [net.c(:)', net.wi, net.w2(:)', net.b2];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function truncated_normal_sparse_grid_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_SPARSE_GRID_TEST tests TRUNCATED_NORMAL_SPARSE_GRID library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_SPARSE_GRID_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TRUNCATED_NORMAL_SPARSE_GRID library.\n' );

  tno_order_test ( );
  tno_test ( );
  tno_persistence_test ( );
  tno_sparse_grid_test ( );

  tno_sparse_grid_size_test ( );
  tno_sparse_grid_test ( )
  tno_sparse_grid_write_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_SPARSE_GRID_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function scale = p26_scale ( neqn )

%*****************************************************************************80
%
%% P26_SCALE returns scale factors for problem p26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end


text:
function DVH_out = DVHwithShifts(planC,structNum,doseNum)
%function DVH_out = DVHwithShifts(planC,structNum,doseNum)
%
%APA, 11/09/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Number of Fractions and Trials
numFractions = 30;
numTrials    = 100;

% Systematic Shifts in cm (Inter-treatment)
XdispSys_Std = 0;
YdispSys_Std = 0;
ZdispSys_Std = 0;

% Random Shifts in cm (inter-fraction)
XdispRnd_Std = 1;
YdispRnd_Std = 1;
ZdispRnd_Std = 1;

indexS = planC{end};

%Distances are in cm
XpdfMean = 0;

YpdfMean = 0;

ZpdfMean = 0;

%Get scan associated with doseNum
scanNum = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID, planC);

if isempty(scanNum) %Assume dose is associated with this scan
    scanNum = getStructureAssociatedScan(structNum,planC);
end

%Get reference transformation matrix for doseNum
if ~isempty(scanNum) && isempty(planC{indexS.dose}(doseNum).transM)
    referenceTransM = planC{indexS.scan}(scanNum).transM;
else
    referenceTransM = planC{indexS.dose}(doseNum).transM;
end

%Store the reference doseArray
referenceDoseArray = planC{indexS.dose}(doseNum).doseArray;

%Try and get a binWidth from stateS.  If it doesnt exist, get it from
%the CERROptions file (allows this function to be called outside CERR)
global stateS;
if ~isempty(stateS) && isfield(stateS, 'optS') && isfield(stateS.optS, 'DVHBinWidth') && ~isempty(stateS.optS.DVHBinWidth)
    binWidth = stateS.optS.DVHBinWidth;
else
    optS = CERROptions;
    binWidth = optS.DVHBinWidth;
end

%Compute DVH at given dose
[dosesCurrentV, volsV] = getDVH(structNum, doseNum, planC);
[doseBinsCurrentV, volsHistCurrentV] = doseHist(dosesCurrentV, volsV, binWidth);

%Divide doseArray into fractions
planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray / numFractions;

%Get systematic errors
deltaX_systematic = randn(1, numTrials) * XdispSys_Std + XpdfMean;
deltaY_systematic = randn(1, numTrials) * YdispSys_Std + YpdfMean;
deltaZ_systematic = randn(1, numTrials) * ZdispSys_Std + ZpdfMean;

hWait = waitbar(0,'Computing plan robustness...');

if isempty(referenceTransM)
    referenceTransM_matrix = eye(4);
else
    referenceTransM_matrix = referenceTransM;
end

%Obtain DVH calculation points in blocks
optS    = planC{indexS.CERROptions};

ROIImageSize = [planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1  planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2];

deltaY = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;

%Get raster segments for structure.
[segmentsM, planC, isError] = getRasterSegments(structNum, planC);

if isempty(segmentsM)
    isError = 1;
end
numSegs = size(segmentsM,1);

%Relative sampling of ROI voxels in this place, compared to CT spacing.
%Set when rasterSegments are generated (usually on import).
sampleRate = optS.ROISampleRate;

%Sample the rows
indFullV =  1 : numSegs;
if sampleRate ~= 1
    rV = 1 : length(indFullV);
    rV([rem(rV+sampleRate-1,sampleRate)~=0]) = [];
    indFullV = rV;
end

%Block process to avoid swamping on large structures
if isfield(optS, 'DVHBlockSize') & ~isempty(optS.DVHBlockSize)
    DVHBlockSize = optS.DVHBlockSize;
else
    DVHBlockSize = 5000;
end

blocks = ceil(length(indFullV)/DVHBlockSize);

start = 1;

for b = 1 : blocks

  %Build the interpolation points matrix

  dummy = zeros(1,DVHBlockSize * ROIImageSize(1));
  x1V = dummy;
  y1V = dummy;
  z1V = dummy;
  volsSectionV =  dummy;

  if start+DVHBlockSize > length(indFullV)
    stop = length(indFullV);
  else
    stop = start + DVHBlockSize - 1;
  end

  indV = indFullV(start:stop);

  mark = 1;
  for i = indV

    tmpV = segmentsM(i,1:10);
    delta = tmpV(5) * sampleRate;
    xV = tmpV(3): delta : tmpV(4);
    len = length(xV);
    rangeV = ones(1,len);
    yV = tmpV(2) * rangeV;
    zV = tmpV(1) * rangeV;
    sliceThickness = tmpV(10);
    %v = delta^2 * sliceThickness;
    v = delta * (deltaY*sampleRate) * sliceThickness;
    x1V(mark : mark + len - 1) = xV;
    y1V(mark : mark + len - 1) = yV;
    z1V(mark : mark + len - 1) = zV;
    volsSectionV(mark : mark + len - 1) = v;
    mark = mark + len;

  end

  %cut unused matrix elements
  x1V = x1V(1:mark-1);
  y1V = y1V(1:mark-1);
  z1V = z1V(1:mark-1);
  volsSectionV = volsSectionV(1:mark-1);

  %Get transformation matrices for both dose and structure.
  transMDose    = getTransM('dose', doseNum, planC);
  transMStruct  = getTransM('struct', structNum, planC);

  %Forward transform the structure's coordinates.
  if ~isempty(transMStruct)
      [x1V, y1V, z1V] = applyTransM(transMStruct, x1V, y1V, z1V);
  end

  dvhCalcPtsC{b} = [x1V', y1V', z1V'];

  %Interpolate.
%   [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

%   dosesV = [dosesV, dosesSectionV];
%   volsV  = [volsV, volsSectionV];

  start = stop + 1;

end

%Loop over to generate DVH
for iTrial = 1:numTrials

    deltaX = deltaX_systematic(iTrial) + randn(1, numFractions) * XdispRnd_Std + XpdfMean;
    deltaY = deltaY_systematic(iTrial) + randn(1, numFractions) * YdispRnd_Std + YpdfMean;
    deltaZ = deltaZ_systematic(iTrial) + randn(1, numFractions) * ZdispRnd_Std + ZpdfMean;

    tmpDoseV = zeros(1,length(dosesCurrentV));

    for iFraction = 1:numFractions

        waitbar(((iTrial-1)*numFractions + iFraction)/(numTrials*numFractions),hWait)

        transM = referenceTransM_matrix;
        transM(1:3,4) = transM(1:3,4) + [deltaX(iFraction); deltaY(iFraction); deltaZ(iFraction)];

        %Apply the new transM to dose
        planC{indexS.dose}(doseNum).transM = transM;

        %Get doses and volumes of points in structure.
        %[dosesV, volsV] = getDVH(structNum, doseNum, planC);

        volsV  = [];
        dosesV = [];

        for b = 1 : blocks
            x1V = dvhCalcPtsC{b}(:,1);
            y1V = dvhCalcPtsC{b}(:,2);
            z1V = dvhCalcPtsC{b}(:,3);

            %Back transform the coordinates into the doses' coordinate system.
            [x1V, y1V, z1V] = applyTransM(inv(transM), x1V, y1V, z1V);

            %Interpolate.
            [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

            dosesV = [dosesV, dosesSectionV];
            volsV  = [volsV, volsSectionV];
        end

        tmpDoseV = tmpDoseV + dosesV;

    end

    DVHm(:,iTrial) = single(tmpDoseV(:));

    %Compute Mean and Std for binned dose
    %for iTrial=1:numTrials
        doseTrialV = DVHm(:,iTrial);
        [doseBinsV, volsHistV] = doseHist(doseTrialV, volsV, binWidth);
        doseBinsTmpV = doseBinsV;
        volsHistTmpV = volsHistV;
        if iTrial > 1
            AddToCurrent = ~ismember(doseBins{1},doseBinsV);
            AddToPrevious = ~ismember(doseBinsV,doseBins{1});
            doseBinsV = [doseBinsV doseBins{1}(AddToCurrent)];
            volsHistV = [volsHistV volsHist{1}(AddToCurrent)];
            for jTrial = 1:iTrial-1
                doseBins{jTrial} = [doseBins{jTrial} doseBinsTmpV(AddToPrevious)];
                volsHist{jTrial} = [volsHist{jTrial} volsHistTmpV(AddToPrevious)];
                [doseBins{jTrial},indSort] = sort(doseBins{jTrial});
                volsHist{jTrial} = volsHist{jTrial}(indSort);
            end
        end
        doseBins{iTrial} = doseBinsV;
        volsHist{iTrial} = volsHistV;
        [doseBins{iTrial},indSort] = sort(doseBins{iTrial});
        volsHist{iTrial} = volsHist{iTrial}(indSort);
    %end

end

close(hWait)

%Compute Std of binned volumes
blockSize = 100;
for iBlock = 1:ceil(length(volsHist{1})/blockSize)
    if iBlock == ceil(length(volsHist{1})/blockSize)
        indicesV = (iBlock-1)*blockSize+1:length(volsHist{1});
    else
        indicesV = (iBlock-1)*blockSize+1:iBlock*blockSize;
    end
    volsHistStdM = [];
    for iTrial = 1:length(volsHist)
        volsHistStdM(:,iTrial) = volsHist{iTrial}(indicesV);
    end
    volsHistStdV(indicesV) = std(volsHistStdM,0,2)';
end

%Reassign reference transM and doseArray to doseNum
planC{indexS.dose}(doseNum).transM = referenceTransM;
planC{indexS.dose}(doseNum).doseArray = referenceDoseArray;

%Compute Mean and Std Dev in blocks
numVoxels = length(volsV);
DVHBlockSize = 5000;
blocks = ceil(numVoxels/DVHBlockSize);
start = 1;

for b = 1 : blocks

  if start+DVHBlockSize > numVoxels
    stop = numVoxels;
  else
    stop = start + DVHBlockSize - 1;
  end

  DVH_block = DVHm(start:stop,:);

  meanDoseV(start:stop) = mean(DVH_block,2);
  stdDoseV(start:stop)  = std(DVH_block,0,2);

  start = stop + 1;

end

DVH_out = [meanDoseV(:)'; volsV(:)'];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%SerialLink.PAY Joint forces due to payload
%
% TAU = R.PAY(W, J) returns the generalised joint force/torques due to a
% payload wrench W (1x6) and where the manipulator Jacobian is J (6xN), and
% N is the number of robot joints.
%
% TAU = R.PAY(Q, W, F) as above but the Jacobian is calculated at pose Q
% (1xN) in the frame given by F which is '0' for world frame, 'e' for
% end-effector frame.
%
% Uses the formula TAU = J'W, where W is a wrench vector applied at the end
% effector, W = [Fx Fy Fz Mx My Mz]'.
%
% Trajectory operation::
%
% In the case Q is MxN or J is 6xNxM then TAU is MxN where each row is the
% generalised force/torque at the pose given by corresponding row of Q.
%
% Notes::
% - Wrench vector and Jacobian must be from the same reference frame.
% - Tool transforms are taken into consideration when F = 'e'.
% - Must have a constant wrench - no trajectory support for this yet.
%
% Author::
% Bryan Moutrie
%
% See also SerialLink.paycap, SerialLink.jacob0, SerialLink.jacobe.

% Copyright (C) Bryan Moutrie, 2013-2015
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.

function tauP = pay(robot, varargin)

    if length(varargin) == 2
        w = varargin{1};
        J = varargin{2};
        n = size(J,2);
    elseif length(varargin) == 3
        q = varargin{1};
        w = varargin{2};
        f = varargin{3};
        n = robot.n;
        J = zeros(6,n,size(q,1));
        switch f
            case '0'
                for i= 1: size(q,1)
                    J(:,:,i) = robot.jacob0(q(i,:));
                end
            case {'n', 'e'}
                for i= 1: size(q,1)
                    J(:,:,i) = robot.jacobe(q(i,:));
                end
        end
    end

    if ~isequal(size(w),[6 1]), error(pHRIWARE('error', 'inputSize')); end
    tauP = -reshape(J(:,:)'*w,n,[])';

end

author: petercorke

repo: robotics-toolbox-matlab


text:
% MACRO_METHODCALL: Script to insert at the beginning of all the brainstorm class functions
%
% DEPRECATED:  Use "eval(macro_method);" instead.

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

% Get current stack trace
dbStack = dbstack();
% Get the call just one level above
if (length(dbStack) > 1)
    strScript = sprintf('%s.m at line %d', dbStack(2).name, dbStack(2).line);
else
    strScript = 'Unknown';
end

% Display warning when this function is called
disp([10 '*******************************************************************************' ...
      10 '*** WARNING: In function: ' strScript ...
      10 '***          The process API changed, due to modifications in Matlab 2016b.' ...
      10 '***          Replace "macro_methodcall;" with "eval(macro_method);"' ...
      10 '***          at the beginning of all your Brainstorm functions.' ...
      10 '*******************************************************************************']);

% Matlab versions: 2006b-2016a
if bst_verlessthan(901)
    % No parameters: nothing to do
    if (nargin == 0)
    % Else : execute appropriate local function
    elseif ischar(varargin{1})
        if (nargout)
            [varargout{1:nargout}] = feval(str2func(varargin{1}), varargin{2:end});
        else
            feval(str2func(varargin{1}), varargin{2:end});
        end
    end

% Matlab versions: >= 2016b
% FiS Incompatibility fix (developed by The Mathworks, 2016)
else
    % Catch stupid warnings after 2016b (macro_methodcall.m not working anymore)
    warning('off', 'MATLAB:lang:GetCallerInfoWarning');
    % Builtin to capture the caller info, the structure contains nargin, nargout and local function handles in the caller.
    callerinfo = builtin('_GetCallerInfo');
    fcnNames = cellfun(@func2str,callerinfo.localFunctions,'UniformOutput',false);

    if(length(varargin) >= 1)
        ind1 = strcmp(fcnNames, varargin{1});
    end
    % FiS Incompatibility fix
    if (callerinfo.nargin == 0)
    % Else : execute appropriate local function
    elseif ischar(varargin{1})
        % FiS Incompatibility fix
        if (callerinfo.nargout)
            % FiS Incompatibility fix
            [varargout{1:callerinfo.nargout}] = feval(callerinfo.localFunctions{ind1}, varargin{2:end});
        else
            % FiS Incompatibility fix
            feval(callerinfo.localFunctions{ind1}, varargin{2:end});
        end
    end
end

author: brainstorm-tools


text:
function [ h, hmin, hmax ] = p09_stepsize ( option )

%*****************************************************************************80
%
%% P09_STEPSIZE returns step sizes for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =    0.300;
  hmin = 0.001;
  hmax = 0.600;

  return
end


text:
function y = smoothRamp(x,alpha,order)
%
% y = smoothRamp(x,alpha,order)
%
% FUNCTION:
%This function is a smooth version of the standard ramp function, which is:
%   ramp(x) = 0 ; (x<0)
%   ramp(x) = x ; (x>=0)
% Smoothing is done using a polynomial transition, the degree of smoothness
% at the transitions is set using order.
%
% INPUTS:
%   x = a vector of matrix of real numbers
%   alpha = smoothing parameter: -alpha<x<alpha will be smoothed
%   order = the order of the smoothing (0=no smoothing)
%           order must be in the set {0,1,..,8}
%
% OUTPUTS:
%   y(x<-alpha) = 0;
%   y(x>alpha) = x;
%   y(else) = smooth_transition
%
        y=x;
        c1 = x < -alpha;
        c2 = x > alpha;
switch order
    case 0 % No Smoothing
        y(x<0)=0;
    case 1
        y=x;
        y(c1)=0;
        p = [1/4,1/2,1/4];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 2
        y=x;
        y(c1)=0;
        p = [-1/16,0,3/8,1/2,3/16];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 3
        y=x;
        y(c1)=0;
        p = [1/32,0,-5/32,0,15/32,1/2,5/32];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 4
        y=x;
        y(c1)=0;
        p = [-5/256,0,7/64,0,-35/128,0,35/64,1/2,35/256];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 5
        y=x;
        y(c1)=0;
        p = [7/512,0,-45/512,0,63/256,0,-105/256,0,315/512,1/2,63/512];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 6
        y=x;
        y(c1)=0;
        p = [-21/2048,0,77/1024,0,-495/2048,0,231/512,0,-1155/2048,0,693/1024,1/2,231/2048];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 7
        y=x;
        y(c1)=0;
        p = [33/4096,0,-273/4096,0,1001/4096,0,-2145/4096,0,3003/4096,0,-3003/4096,0,3003/4096,1/2,429/4096];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    case 8
        y=x;
        y(c1)=0;
        p = [-429/65536,0,495/8192,0,-4095/16384,0,5005/8192,0,-32175/32768,0,9009/8192,0,-15015/16384,0,6435/8192,1/2,6435/65536];
        idx = ~c1&~c2;
        y(idx)=alpha*polyval(p,x(idx)/alpha);
    otherwise
        error('Order not supported');
end
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function range = p09_range ( )

%*****************************************************************************80
%
%% P09_RANGE returns an interval bounding the root for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = -  5.0;
  range(2) = + 16.0;

  return
end


text:
function name = p09_name ( )

%*****************************************************************************80
%
%% P09_NAME returns the name of problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'Patterson #7';

  return
end