text:
function straddling_graph(blind_time,frame_rate,pulse_distance,laser_energy,num_exposures_to_show,is_dbl_shutter,f1exp_cam)
if is_dbl_shutter % a double shutter camera like the pco panda and pixelfly capture 2 images for each frame.
	frame_rate=frame_rate*2;
end
period = 1/frame_rate*1000^2; % µs
cam_on = period-blind_time; % µs
pulse_length = laser_energy/100*pulse_distance; % µs
max_possible_pulse_distance = pulse_distance - (1/frame_rate*1000^2-cam_on); % µs
if pulse_length > max_possible_pulse_distance
	pulse_length = max_possible_pulse_distance;
end
duty_cycle=pulse_length/period;
if duty_cycle > 0.5
	pulse_length = period / 2;
	duty_cycle=pulse_length/period;
end

exposure_active_x=[];
exposure_active_y=[];
laser_active_x=[];
laser_active_y=[];

start_of_pulse1 = period-pulse_distance/2-pulse_length/2-blind_time/2;
end_of_pulse1 = period-pulse_distance/2-pulse_length/2-blind_time/2+pulse_length;

start_of_pulse2 = period-pulse_distance/2-pulse_length/2-blind_time/2+pulse_distance;
end_of_pulse2 = period-pulse_distance/2-pulse_length/2-blind_time/2+pulse_distance+pulse_length;

bot_cam = 0.45;
bot_las = 0.2;
amp = 0.2;

if is_dbl_shutter
	for i=0:2:num_exposures_to_show+4
		exposure_active_y= [exposure_active_y  [  bot_cam  bot_cam      bot_cam                            bot_cam+amp                    bot_cam+amp                        bot_cam                         bot_cam                      bot_cam+amp                 bot_cam+amp       bot_cam          ]];
		exposure_active_x =[exposure_active_x  [  0           0     start_of_pulse1-blind_time  start_of_pulse1-blind_time  start_of_pulse1+f1exp_cam-blind_time  start_of_pulse1+f1exp_cam-blind_time  start_of_pulse1+f1exp_cam-blind_time   start_of_pulse1+f1exp_cam-blind_time   end_of_pulse2+5000  end_of_pulse2+5000    ] +  period* i];
	end
else
	for i=0:num_exposures_to_show+4
		exposure_active_y=[exposure_active_y  [  bot_cam   bot_cam+amp       bot_cam+amp                      bot_cam               bot_cam    ]];
		exposure_active_x =[exposure_active_x [  0              0          period-blind_time            period-blind_time            period    ] +  period* i];
	end
end

for i=0:2:num_exposures_to_show+4
	laser_active_y= [laser_active_y  [  bot_las  bot_las      bot_las        bot_las+amp    bot_las+amp       bot_las         bot_las         bot_las+amp      bot_las+amp       bot_las          ]];
	laser_active_x =[laser_active_x  [  0           0     start_of_pulse1  start_of_pulse1  end_of_pulse1  end_of_pulse1  start_of_pulse2   start_of_pulse2   end_of_pulse2   end_of_pulse2    ] +  period* i];
end

straddling_figure=findobj('tag','straddling_figure');
if isempty(straddling_figure)
	hf = figure('numbertitle','off','MenuBar','figure','DockControls','off','Name','Camera exposure and pulse timing visualization','Toolbar','figure','CloseRequestFcn', @straddling_figure_CloseRequestFcn,'tag','straddling_figure','visible','on');
else
	hf = figure(straddling_figure);
	clf(hf)
end

plot(exposure_active_x,exposure_active_y,'linewidth',2)
hold on
plot(laser_active_x,laser_active_y,'linewidth',2)
axis tight
xlim([0 period*num_exposures_to_show])
%xlim([-period period*(num_cycles_to_show)+period])
ylim([0 1])
xlabel('time in µs')
ylabel([])
yticks([])

%remove unnecessary Toolbar stuff
% Get a handle to the standard plot toolbar.
tbh = findall(hf,'Type','uitoolbar');
% Get handles to each button we don't want on the standard toolbar and DELETE
ptPlotToolsOn    = findall(hf,'Tag','Plottools.PlottoolsOn');     delete(ptPlotToolsOn);
ptPlotToolsOff   = findall(hf,'Tag','Plottools.PlottoolsOff');    delete(ptPlotToolsOff);
ttInsertLegend   = findall(hf,'Tag','Annotation.InsertLegend');   delete(ttInsertLegend);
ttInsertColorbar = findall(hf,'Tag','Annotation.InsertColorbar'); delete(ttInsertColorbar);
ttLinking        = findall(hf,'Tag','DataManager.Linking');       delete(ttLinking);
ttRotate         = findall(hf,'Tag','Exploration.Rotate');        delete(ttRotate);
ttEditPlot       = findall(hf,'Tag','Standard.EditPlot');         delete(ttEditPlot);
ttEditPlot       = findall(hf,'Tag','Standard.PropertyInspector');         delete(ttEditPlot);
ptPrintFigure    = findall(hf,'Tag','Standard.PrintFigure');      delete(ptPrintFigure);
ptSaveFigure     = findall(hf,'Tag','Standard.SaveFigure');       delete(ptSaveFigure);
ptFileOpen       = findall(hf,'Tag','Standard.FileOpen');         delete(ptFileOpen);
ptNewFigure      = findall(hf,'Tag','Standard.NewFigure');        delete(ptNewFigure);

%Pulse_length
ha = annotation('doublearrow','Head1Style','plain','Head2Style','plain','Head1Length',5,'Head1Width',5,'Head2Length',5,'Head2Width',5);
ha.Parent=hf.CurrentAxes;
ha.X=[start_of_pulse1 end_of_pulse1];
ha.Y=[bot_las-0.01 bot_las-0.01];
plot ([start_of_pulse1 start_of_pulse1], [bot_las-0.01 bot_las],'LineStyle','--','Color',[1 0.7 0.7])
plot ([end_of_pulse1 end_of_pulse1], [bot_las-0.01 bot_las],'LineStyle','--','Color',[1 0.7 0.7])
text((start_of_pulse1+end_of_pulse1)/2,bot_las-0.01,{'Pulse length' [num2str(pulse_length) ' µs']},'Rotation',0,'HorizontalAlignment','center','VerticalAlignment','top','FontSize',8)

%Pulse distance
ha = annotation('doublearrow','Head1Style','plain','Head2Style','plain','Head1Length',5,'Head1Width',5,'Head2Length',5,'Head2Width',5);
ha.Parent=hf.CurrentAxes;
ha.X=[(start_of_pulse1+end_of_pulse1)/2 (start_of_pulse2+end_of_pulse2)/2];
ha.Y=[bot_las-0.1 bot_las-0.1];
plot ([(start_of_pulse1+end_of_pulse1)/2 (start_of_pulse1+end_of_pulse1)/2], [bot_las-0.1 bot_las+amp-0.01],'LineStyle','--','Color',[1 0.7 0.7])
plot ([(start_of_pulse2+end_of_pulse2)/2 (start_of_pulse2+end_of_pulse2)/2], [bot_las-0.1 bot_las+amp-0.01],'LineStyle','--','Color',[1 0.7 0.7])
text(((start_of_pulse1+end_of_pulse1)/2 + (start_of_pulse2+end_of_pulse2)/2)/2,bot_las-0.1,{'Pulse distance' [num2str(pulse_distance) ' µs']},'Rotation',0,'HorizontalAlignment','center','VerticalAlignment','top','FontSize',8)

%frame rate annotation
ha = annotation('doublearrow','Head1Style','plain','Head2Style','plain','Head1Length',5,'Head1Width',5,'Head2Length',5,'Head2Width',5);
ha.Parent=hf.CurrentAxes;
if is_dbl_shutter
	%frame rate double shutter cameras
	ha.X=[start_of_pulse1-blind_time  start_of_pulse1-blind_time+period*2];
	ha.Y=[bot_cam-0.01 bot_cam-0.01];
	plot ([start_of_pulse1-blind_time start_of_pulse1-blind_time], [bot_cam bot_cam-0.01],'LineStyle','--','Color',[0.7 0.7 1])
	plot ([start_of_pulse1-blind_time+period*2 start_of_pulse1-blind_time+period*2], [bot_cam bot_cam-0.01],'LineStyle','--','Color',[0.7 0.7 1])
	text(period*1.5,bot_cam-0.01,{'Double frame rate' [num2str(frame_rate/2) ' Hz'] [num2str(period*2) ' µs']},'Rotation',0,'HorizontalAlignment','center','VerticalAlignment','top','FontSize',8)
else
	%frame rate non-double shutter cameras
	ha.X=[period period*2];
	ha.Y=[bot_cam-0.01 bot_cam-0.01];
	plot ([period period], [bot_cam bot_cam-0.01],'LineStyle','--','Color',[0.7 0.7 1])
	plot ([period*2 period*2], [bot_cam bot_cam-0.01],'LineStyle','--','Color',[0.7 0.7 1])
	text(period*1.5,bot_cam-0.01,{'Frame rate' [num2str(frame_rate) ' Hz'] [num2str(period) ' µs']},'Rotation',0,'HorizontalAlignment','center','VerticalAlignment','top','FontSize',8)
end

%info top left
text(period/10,0.99,{['Laser duty cycle: ' num2str(round(duty_cycle*100,2)) ' %'] ['PIV data rate: ' num2str(frame_rate/2) ' Hz']},'HorizontalAlignment','left','VerticalAlignment','top')
%info top right
text(num_exposures_to_show*period,0.99,{' Use zoom buttons to see the details '},'HorizontalAlignment','right','VerticalAlignment','top')

margin=period/20*0;
smallmargin=margin/2*0;
for i=0:2:num_exposures_to_show-1

	ha=annotation('textbox',[0,0,0,0],'String',['Image Pair ' num2str(i/2+1)],'HorizontalAlignment','center','VerticalAlignment','top','BackgroundColor','k','FaceAlpha',0.1);
	ha.Parent=hf.CurrentAxes;
	ha.Position=[0+margin+  period* i           bot_cam+amp+0.04              period*2-margin*2                 0.175];

	ha=annotation('textbox',[0,0,0,0],'String',{'Image A' ['PIVlab_' sprintf('%3.3d',i/2) '_A.tif']},'HorizontalAlignment','center','VerticalAlignment','middle','BackgroundColor','k','FaceAlpha',0.1,'LineStyle','-','FontSize',8,'Interpreter','none');
	ha.Parent=hf.CurrentAxes;
	ha.Position=[0+margin+smallmargin + period* i              bot_cam+amp+0.05                period-margin*2-smallmargin*2               0.1];

	ha=annotation('textbox',[0,0,0,0],'String',{'Image B' ['PIVlab_' sprintf('%3.3d',i/2) '_B.tif']},'HorizontalAlignment','center','VerticalAlignment','middle','BackgroundColor','k','FaceAlpha',0.1,'LineStyle','-','FontSize',8,'Interpreter','none');
	ha.Parent=hf.CurrentAxes;
	ha.Position=[period+margin+smallmargin + period* i                    bot_cam+amp+0.05              period-margin*2-smallmargin*2              0.1];

end

legend({'Camera Exposure', 'Laser Active'},'Location','southeast')
set(gca,'InnerPosition',[0.01 0.1 0.98 0.9])
%grid on
%set(gca,'YGrid','off')
%set(gcf,'handlevisibility','callback')
end

function straddling_figure_CloseRequestFcn(hObject, ~, ~)
try
	delete(hObject);
catch
	delete(gcf);
end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function y=savitzkyGolayFilt(x,N,DN,F,W,DIM)
%savitzkyGolayFilt Savitzky-Golay Filtering.
%   savitzkyGolayFilt(X,N,DN,F) filters the signal X using a Savitzky-Golay
%   (polynomial) filter.  The polynomial order, N, must be less than the
%   frame size, F, and F must be odd.  DN specifies the differentiation
%   order (DN=0 is smoothing). For a DN higher than zero, you'll have to
%   scale the output by 1/T^DN to acquire the DNth smoothed derivative of
%   input X, where T is the sampling interval. The length of the input X
%   must be >= F.  If X is a matrix, the filtering is done on the columns
%   of X.
%
%   Note that if the polynomial order N equals F-1, no smoothing
%   will occur.
%
%   savitzkyGolayFilt(X,N,DN,F,W) specifies a weighting vector W with
%   length F containing real, positive valued weights employed during the
%   least-squares minimization. If not specified, or if specified as
%   empty, W defaults to an identity matrix.
%
%   savitzkyGolayFilt(X,N,DN,F,[],DIM) or savitzkyGolayFilt(X,N,DN,F,W,DIM)
%   operates along the dimension DIM.
%
%   See also savitzkyGolay, FILTER, sgolayfilt

%   References:
%     [1] Sophocles J. Orfanidis, INTRODUCTION TO SIGNAL PROCESSING,
%              Prentice-Hall, 1995, Chapter 8.

%   Author(s): R. Losada
%   Copyright 1988-2004 The MathWorks, Inc.
%   $Revision: 1.11.4.4 $  $Date: 2009/08/11 15:47:54 $

error(nargchk(4,6,nargin,'struct'));

% Check if the input arguments are valid
if round(F) ~= F, error(generatemsgid('MustBeInteger'),'Frame length must be an integer.'), end
if rem(F,2) ~= 1, error(generatemsgid('SignalErr'),'Frame length must be odd.'), end
if round(N) ~= N, error(generatemsgid('MustBeInteger'),'Polynomial order must be an integer.'), end
if N > F-1, error(generatemsgid('InvalidRange'),'The Polynomial order must be less than the frame length.'), end
if DN > N, error(generatemsgid('InvalidRange'),'The Differentiation order must be less than or equal to the Polynomial order.'), end

if nargin < 5 || isempty(W)
   % No weighting matrix, make W an identity
   W = ones(F,1);
else
   % Check for right length of W
   if length(W) ~= F, error(generatemsgid('InvalidDimensions'),'The weight vector must be of the same length as the frame length.'),end
   % Check to see if all elements are positive
   if min(W) <= 0, error(generatemsgid('InvalidRange'),'All the elements of the weight vector must be greater than zero.'), end
end

if nargin < 6, DIM = []; end

% Compute the projection matrix B
pp = fix(-F./2):fix(F./2);
B = savitzkyGolay(pp,N,DN,pp,W);

if ~isempty(DIM) && DIM > ndims(x)
	error(generatemsgid('InvalidDimensions'),'Dimension specified exceeds the dimensions of X.')
end

% Reshape X into the right dimension.
if isempty(DIM)
	% Work along the first non-singleton dimension
	[x, nshifts] = shiftdim(x);
else
	% Put DIM in the first dimension (this matches the order
	% that the built-in filter function uses)
	perm = [DIM,1:DIM-1,DIM+1:ndims(x)];
	x = permute(x,perm);
end

if size(x,1) < F, error(generatemsgid('InvalidDimensions'),'The length of the input must be >= frame length.'), end

% Preallocate output
y = zeros(size(x));

% Compute the transient on (note, this is different than in sgolayfilt,
% they had an optimization leaving out some transposes that is only valid
% for DN==0)
y(1:(F+1)/2-1,:) = fliplr(B(:,(F-1)/2+2:end)).'*flipud(x(1:F,:));

% Compute the steady state output
ytemp = filter(B(:,(F-1)./2+1),1,x);
y((F+1)/2:end-(F+1)/2+1,:) = ytemp(F:end,:);

% Compute the transient off
y(end-(F+1)/2+2:end,:) = fliplr(B(:,1:(F-1)/2)).'*flipud(x(end-(F-1):end,:));

% Convert Y to the original shape of X
if isempty(DIM)
	y = shiftdim(y, -nshifts);
else
	y = ipermute(y,perm);
end

author: Sable

repo: mcbench-benchmarks


text:
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%

%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  table = zeros(m,n);

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [vw, status, forceSave] = saveROI(vw, ROI, local, forceSave)
%
% [view, status, forceSave] = saveROI([view=cur view], [ROI='selected'], [local], [forceSave=0])
%
% Saves ROI to a file.
%
% INPUTS:
%   vw: mrVista view. Finds current view if omitted.
%
%   ROI: ROI structure, or index into the ROI structure in the view.
%   If this the string 'selected' the view's selected ROI to be saved.
%   [This is the default behavior if the ROI argument is omitted.]
%
% local: flag to save ROI in a local directory to the view (e.g.,
% 'Volume/ROIs'), or a shared directory (e.g., relative to the
% anatomy directory -- use "setpref('VISTA', 'defaultROIPath')" to
% set this directory.)  [For inplanes, defaults to 1, save locally.
% for other views, defaults to 0, saves in shared directory.]
%
% forceSave: if set to 0, if the ROI file exists, will ask the user if they
% want to save over the existing file. If 1, will save over existing files
% without asking. Useful for scripting. [Defaults to 0].
%
% OUTPUTS:
%   vw    : Modified view (e.g, Inplane ...).
%   status: 0 if the ROI saving was aborted. This is used when calling this
%           repeatedly (see 'saveAllROIs'), if the file exists and the
%           user selects 'Cancel'.
%   forceSave: modified forceSave flag, if the file exists, and the user
%           wants to bypass the dialog for future files. (again, see
%           'saveAllROIs'.)
%
% Example:  saveROI(INPLANE{1});
%
% See also:
%    The roi management is a mess.  This should be roiSave, but there is
%    another roiSave.  All the roi management routines have different names
%    and different behavior.  Sigh.
%
% Vistasoft Team
%
% djh, 1/24/98
% gmb, 4/25/98 added dialog box
% ARW 10/06/05 : Added default ROI path option
% (see also roiDir)
% ras 07/06: added forceSave flag.
% ras 10/06: ROI selection defaults to 'selected'; can enter ROI indices.
% ras 11/06: offers option to save over all, returning modified forceSave
% flag. (Very helpful when saving a lot of ROIs that you want to update.)

if notDefined('vw'),        vw = getCurView;        end
if notDefined('forceSave'), forceSave = 0;          end
if notDefined('ROI'),       ROI = 'selected';       end
if notDefined('local'),     local = isequal(vw.viewType, 'Inplane'); end

status = 0;
verbose = prefsVerboseCheck;

% disambiguate the ROI specification
if isnumeric(ROI), ROI = tc_roiStruct(vw, ROI); end

if ischar(ROI) && isequal(lower(ROI), 'selected')
    ROI = vw.ROIs( vw.selectedROI );
end

if isfield(ROI, 'roiVertInds')
    ROI = rmfield(ROI,'roiVertInds');
end

if isempty(ROI.coords) && forceSave==0
    q = sprintf('Save empty ROI %s?', ROI.name);
    confirm = questdlg(q, mfilename, 'Yes', 'No', 'No');
    if ~isequal(confirm, 'Yes'), disp('ROI not saved.'); return; end
end

% The method of saving the ROI starts here
% First we make sure we add a .mat extension. Matlab does this most of
% the times but not always, for example if your name contains a '.'.
% It still saves the file but without '.mat' and consequently cannot load it.
[~,~,ext]=fileparts(ROI.name);
if ~strcmp(ext,'.mat'), ROIname = [ROI.name '.mat'];
else                    ROIname = ROI.name;     end;

pathStr = fullfile(roiDir(vw,local), ROIname);

if check4File(pathStr) && forceSave==0
    q = sprintf('ROI %s already exists.  Overwrite?', ROI.name);
    saveFlag = questdlg(q, 'Save ROI', 'Yes', 'No', 'Yes To All', 'No');
else
    saveFlag = 'Yes';
end

switch saveFlag
    case 'Yes'
        if verbose,
            fprintf('Saving ROI "%s" in %s. \n', ROI.name, fileparts(pathStr));
        end
        save(pathStr,'ROI');
        status = 1;

    case 'No'
        if verbose,  fprintf('ROI %s not saved.', ROI.name);    end
        status = 1;     % not saved, but keep going if saving many ROIs

    case 'Yes To All'
        if verbose,     disp('Force-Saving all ROIs...');       end
        fprintf('Saving ROI "%s" in %s. \n', ROI.name, fileparts(pathStr));
        forceSave = 1;
        status = 1;

    case 'Cancel'
        if verbose,  disp('ROI Saving Aborted.');               end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% stein9 0-1 integer program from MIPLIB 2.0 by George Nemhauser et al.,
% solving a nine item Steiner triple system
% See, e.g., http://miplib.zib.de/miplib2/miplib/stein9.mps.gz

clear;

% the (dense) 0-1 matrix with the constraints
% the first twelve constraints are of the form <sum of three variables greater equal one>
% the last constraints gives a trivial bound on the objective and is actually redundant
matrix = [
0, 1, 1, 1, 0, 0, 0, 0, 0;...
1, 0, 1, 0, 1, 0, 0, 0, 0;...
1, 1, 0, 0, 0, 1, 0, 0, 0;...
0, 0, 0, 0, 1, 1, 1, 0, 0;...
0, 0, 0, 1, 0, 1, 0, 1, 0;...
0, 0, 0, 1, 1, 0, 0, 0, 1;...
1, 0, 0, 0, 0, 0, 0, 1, 1;...
0, 1, 0, 0, 0, 0, 1, 0, 1;...
0, 0, 1, 0, 0, 0, 1, 1, 0;...
1, 0, 0, 1, 0, 0, 1, 0, 0;...
0, 1, 0, 0, 1, 0, 0, 1, 0;...
0, 0, 1, 0, 0, 1, 0, 0, 1;...
1, 1, 1, 1, 1, 1, 1, 1, 1;...
];

% all constraints are of the type a_1*x_1+...+a9*x_9 >= b
% In SCIP, constraints are represented as left hand side <= linear sum <= right hand side
% Hence, the left hand sides are all finite, the right hand sides are plus infinity
lhs = [1,1,1,1,1,1,1,1,1,1,1,1,4]';
rhs = [1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20,1e+20]';

% the variables are binary, hence, the lower bounds are zero, the upper bounds are one
lb = [0,0,0,0,0,0,0,0,0]';
ub = [1,1,1,1,1,1,1,1,1]';

% all variables are binary decision variables (you put an item into the system or not)
% the goal is to minimize the number of used variables/items
vartype = ['b','b','b','b','b','b','b','b','b'];
obj = [1,1,1,1,1,1,1,1,1]';
objsense = 'min';

% call SCIP
[bestsol, objval] = scip(matrix, lhs, rhs, obj, lb, ub, vartype, objsense);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function L = log_marg_prob_node_case(CPD, y, x)
% LOG_MARG_PROB_NODE_CASE Compute prod_m log P(x(i,m)| x(pi_i,m)) for node i (tabular)
% L = log_marg_prob_node_case(CPD, self_ev, parent_ev)
%
% This is a slightly optimised version of log_marg_prob_node.
% We assume we have exactly 1 case, i.e., y is a scalar and x is a vector (not a cell array).

sz = CPD.sizes;
nparents = length(sz)-1;

% We assume the CPTs are already set to the mean of the posterior (due to update_params)

switch nparents
 case 0, p = CPD.CPT(y);
 case 1, p = CPD.CPT(x(1), y);
 case 2, p = CPD.CPT(x(1), x(2), y);
 case 3, p = CPD.CPT(x(1), x(2), x(3), y);
 otherwise,
  ind = subv2ind(sz, [x y]);
  p = CPD.CPT(ind);
end
L = log(p);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% ASTROTIK by Francesco Santilli
% R2BP (Restricted Two Bodies Problem)
% Compute the Minimum Orbit Intersection Distance (MOID)
%
% Usage: [d,f] = moid(orbitA,orbitB)
%
% where: orbit = [p e i o w] = 3d orbit elements
%        orbit = [p e w] = 2d orbit elements
%           p = semi-latus rectum [L] (p>0)
%           e = eccentricity [-] (0<=e<1)
%           i = inclinaion [rad]
%           o = raan [rad]
%           w = argument of perifocus [rad]
%        d = MOID [L]
%        f = [fA fB] = true anomalies [rad]

function [d,f] = moid(orbitA,orbitB)

    if ~(nargin == 2)
        error('Wrong number of input arguments.')
    end

    DA = check(orbitA,1);
    DB = check(orbitB,1);

    if DA ~= DB
        error('Wrong size of input arguments.')
    end

    if ~(DA==3 || DA==5)
        error('Wrong size of input arguments.')
    end
    d3 = (DA==5);

    % rotation matrix
    if d3
        [pA,eA,iA,oA,wA] = take(orbitA);
        [pB,eB,iB,oB,wB] = take(orbitB);
        RA = rotation([oA iA wA]);
        RB = rotation([oB iB wB]);
        RA = RA(1:2,:);
        RB = RB(1:2,:);
    else
        [pA,eA,wA] = take(orbitA);
        [pB,eB,wB] = take(orbitB);
        cwA = cos(wA);
        swA = sin(wA);
        cwB = cos(wB);
        swB = sin(wB);
        RA = [ cwA swA
             -swA cwA];
        RB = [ cwB swB
              -swB cwB];
    end

    if pA<=0 || pB<=0
        error('p must be a stricly positive value.')
    end

    if eA<0 || eB<0 || eA>=1 || eB>=1
        error('e must be in the range [0,1).')
    end

    L2 = pA^2+pB^2;
    f0 = [+1 +1 -1 -1
          +1 -1 +1 -1]*pi/2;

    f = zeros(2,4);
    d = zeros(1,4);
    ef = zeros(1,4);

    opt = optimset('LargeScale','on','Display','off','GradObj','on',...
                   'TolX',eps,'TolFun',0,... % 'Hessian','on',
                   'MaxFunEvals',Inf,'MaxIter',Inf);
    for k = 1:4
        try
            [f(:,k),d(k),ef(k)] = fminunc(@fun,f0(:,k),opt);
        catch
            f(:,k) = [NaN NaN]';
            d(k) = NaN;
            ef(k) = 0;
        end
    end

    bad = (ef==0);
    f(:,bad) = [];
    d(:,bad) = [];
    if isempty(f)
        error('No solution found.')
    end

    d = sqrt(d*L2);
    f = mod(f+pi,2*pi)-pi;
    [d,k] = min(d);
    f = f(:,k(1))';

    function [F,FF,FFF] = fun(f)

        fA = f(1);
        fB = f(2);

        cA = cos(fA);
        cB = cos(fB);
        sA = sin(fA);
        sB = sin(fB);

        rA = pA/(1+eA*cA);
        rB = pB/(1+eB*cB);

        A = rA*[cA sA]*RA;
        B = rB*[cB sB]*RB;

        D = A-B;
        F = D*D' / L2;

        if nargout > 1

            kA = rA^2/pA;
            kB = rB^2/pB;

            GA = [-sA eA+cA]*RA;
            GB = [-sB eB+cB]*RB;

            FA =  2*kA*D*GA';
            FB = -2*kB*D*GB';
            FF = [FA; FB] / L2;

            if nargout > 2

                kAA = 2*kA^2/rA*eA*sA;
                kBB = 2*kB^2/rB*eB*sB;

                GAA = -[cA sA]*RA;
                GBB = -[cB sB]*RB;

                FAA =  2*(kAA*D*GA' + kA*D*GAA' + kA^2*GA*GA');
                FBB = -2*(kBB*D*GB' + kB*D*GBB' + kB^2*GB*GB');
                FAB = -2*kA*kB*GA*GB';

                FFF = [FAA FAB
                       FAB FBB] / L2;

            end
        end

        %fprintf('%.20f %.20f -> %.20f\n',f',F);

    end

end


text:
classdef DescriptorExtractor < handle
    %DESCRIPTOREXTRACTOR  Common interface of 2D image Descriptor Extractors
    %
    % Class for computing descriptors for image keypoints.
    %
    % Extractors of keypoint descriptors in OpenCV have wrappers with a common
    % interface that enables you to easily switch between different algorithms
    % solving the same problem. This section is devoted to computing
    % descriptors represented as vectors in a multidimensional space. All
    % objects that implement the vector descriptor extractors inherit the
    % DescriptorExtractor interface.
    %
    % In this interface, a keypoint descriptor can be represented as a dense,
    % fixed-dimension vector of a basic type. Most descriptors follow this
    % pattern as it simplifies computing distances between descriptors.
    % Therefore, a collection of descriptors is represented as a matrix,
    % where each row is a keypoint descriptor.
    %
    % ## Example
    %
    %     detector = cv.FeatureDetector('ORB');
    %     keypoints = detector.detect(img);
    %     extractor = cv.DescriptorExtractor('ORB');
    %     descriptors = extractor.compute(img, keypoints);
    %
    % See also: cv.FeatureDetector, cv.DescriptorMatcher,
    %  cv.BOWImgDescriptorExtractor, extractFeatures, binaryFeatures
    %

    properties (SetAccess = private)
        % Object ID
        id
        % Type of the extractor
        Type
    end

    methods
        function this = DescriptorExtractor(extractorType, varargin)
            %DESCRIPTOREXTRACTOR  Creates a descriptor extractor by name
            %
            %     extractor = cv.DescriptorExtractor(type)
            %     extractor = cv.DescriptorExtractor(type, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __type__ The following extractor types are supported:
            %   * __BRISK__ see cv.BRISK
            %   * __ORB__ see cv.ORB (default)
            %   * __KAZE__ see cv.KAZE
            %   * __AKAZE__ see cv.AKAZE
            %   * __SIFT__ see cv.SIFT (requires `xfeatures2d` module)
            %   * __SURF__ see cv.SURF (requires `xfeatures2d` module)
            %   * __FREAK__ see cv.FREAK (requires `xfeatures2d` module)
            %   * __BriefDescriptorExtractor__ see cv.BriefDescriptorExtractor
            %     (requires `xfeatures2d` module)
            %   * __LUCID__ see cv.LUCID (requires `xfeatures2d` module)
            %   * __LATCH__ see cv.LATCH (requires `xfeatures2d` module)
            %   * __DAISY__ see cv.DAISY (requires `xfeatures2d` module)
            %   * __VGG__ see cv.VGG (requires `xfeatures2d` module)
            %   * __BoostDesc__ see cv.BoostDesc (requires `xfeatures2d` module)
            %
            % ## Options
            % Refer to the constructors of each descriptor extractor for a
            % list of supported options.
            %
            % See also: cv.DescriptorExtractor.compute
            %
            if nargin < 1, extractorType = 'ORB'; end
            this.Type = extractorType;
            this.id = DescriptorExtractor_(0, 'new', extractorType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     extractor.delete()
            %
            % See also: cv.DescriptorExtractor
            %
            if isempty(this.id), return; end
            DescriptorExtractor_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = extractor.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = DescriptorExtractor_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     extractor.clear()
            %
            % See also: cv.DescriptorExtractor.empty,
            %  cv.DescriptorExtractor.load
            %
            DescriptorExtractor_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = extractor.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.DescriptorExtractor.clear,
            %  cv.DescriptorExtractor.load
            %
            b = DescriptorExtractor_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     extractor.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.DescriptorExtractor.load
            %
            DescriptorExtractor_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     extractor.load(fname)
            %     extractor.load(str, 'FromString',true)
            %     extractor.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.DescriptorExtractor.save
            %
            DescriptorExtractor_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = extractor.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.DescriptorExtractor.save,
            %  cv.DescriptorExtractor.load
            %
            name = DescriptorExtractor_(this.id, 'getDefaultName');
        end
    end

    %% Features2D
    methods
        function ntype = defaultNorm(this)
            %DEFAULTNORM  Returns the default norm type
            %
            %     ntype = extractor.defaultNorm()
            %
            % ## Output
            % * __ntype__ Norm type. One of `cv::NormTypes`:
            %   * __Inf__
            %   * __L1__
            %   * __L2__
            %   * __L2Sqr__
            %   * __Hamming__
            %   * __Hamming2__
            %
            % See also: cv.DescriptorExtractor.compute, cv.DescriptorMatcher
            %
            ntype = DescriptorExtractor_(this.id, 'defaultNorm');
        end

        function sz = descriptorSize(this)
            %DESCRIPTORSIZE  Returns the descriptor size in bytes
            %
            %     sz = extractor.descriptorSize()
            %
            % ## Output
            % * __sz__ Descriptor size.
            %
            % See also: cv.DescriptorExtractor.descriptorType,
            %  cv.DescriptorExtractor.compute
            %
            sz = DescriptorExtractor_(this.id, 'descriptorSize');
        end

        function dtype = descriptorType(this)
            %DESCRIPTORTYPE  Returns the descriptor type
            %
            %     dtype = extractor.descriptorType()
            %
            % ## Output
            % * __dtype__ Descriptor type, one of numeric MATLAB class names.
            %
            % See also: cv.DescriptorExtractor.descriptorSize,
            %  cv.DescriptorExtractor.compute
            %
            dtype = DescriptorExtractor_(this.id, 'descriptorType');
        end

        function [descriptors, keypoints] = compute(this, img, keypoints)
            %COMPUTE  Computes the descriptors for a set of keypoints detected in an image or image set
            %
            %     [descriptors, keypoints] = extractor.compute(img, keypoints)
            %     [descriptors, keypoints] = extractor.compute(imgs, keypoints)
            %
            % ## Input
            % * __img__ Image (first variant).
            % * __imgs__ Image set (second variant), cell array of images.
            % * __keypoints__ Input collection of keypoints. Keypoints for
            %   which a descriptor cannot be computed are removed and the
            %   remaining ones may be reordered. Sometimes new keypoints can
            %   be added, for example: cv.SIFT duplicates keypoint with
            %   several dominant orientations (for each orientation). In the
            %   first variant, this is a struct-array of detected keypoints.
            %   the second variant, it is a cell-array, where `keypoints{i}`
            %   is a set of keypoints detected in `images{i}` (a struct-array
            %   like before). Each keypoint is a struct with the following
            %   fields:
            %   * __pt__ coordinates of the keypoint `[x,y]`
            %   * __size__ diameter of the meaningful keypoint neighborhood
            %   * __angle__ computed orientation of the keypoint (-1 if not
            %     applicable); it's in [0,360) degrees and measured relative
            %     to image coordinate system (y-axis is directed downward),
            %     i.e in clockwise.
            %   * __response__ the response by which the most strong keypoints
            %     have been selected. Can be used for further sorting or
            %     subsampling.
            %   * __octave__ octave (pyramid layer) from which the keypoint
            %     has been extracted.
            %   * **class_id** object class (if the keypoints need to be
            %     clustered by an object they belong to).
            %
            % ## Output
            % * __descriptors__ Computed descriptors. In the second variant of
            %   the method `descriptors{i}` are descriptors computed for a
            %   `keypoints{i}`. Row `j` in `descriptors` (or `descriptors{i}`)
            %   is the descriptor for `j`-th keypoint.
            % * __keypoints__ Optional output with possibly updated keypoints.
            %
            % See also: cv.DescriptorExtractor.DescriptorExtractor
            %
            [descriptors, keypoints] = DescriptorExtractor_(this.id, 'compute', img, keypoints);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function calcAllFusedTextures_HN(pathData,pathText,namePT,nameCT,nameROI,featType,CTweight_mat,scale_mat,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% function calcAllFusedTextures_HN(pathData,pathText,namePT,nameCT,nameROI,featType,CTweight_mat,scale_mat,Ng_mat)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes FUSED texture features for all patients, for all
% different combinations of the following texture extraction parameters:
% - CT weight: Weight given to MRI wavelet low-pass sub-bands in the
%               PET/CT fusion process.
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for FUSED scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM  imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/10.7937/K9/xxxxxxxxxxxxxxxxxxx, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathData: Full path to the HN sData files directory.
%              --> Ex: '/myProject/WORKSPACE/DATA'
% 2. pathNonText: Full path to the HN non texture features directory.
%              --> Ex: '/myProject/WORKSPACE/FEATURES/NON_TEXTURES'
% 3. namePT: Cell of strings of all PET sData files to read
%            --> Ex: {'HGJ_001_PT.PTscan.mat';'HGJ_022_PT.PTscan.mat'}
% 4. namePT: Cell of strings of all CT sData files to read
%            --> Ex: {'HGJ_001_CT.CTscan.mat';'HGJ_022_CT.CTscan.mat'}
% 5. nameROI: Cell of strings specifying the ROI names to analyze for the
%             patients defined by "namePT" and "nameCT"
%             --> Ex: {'GTV';'GTV-P'}
% 6. featType: Either 'GTVp' for primary GTV, or 'GTVtot' for primaty GTV +
%              nodal GTVs
% 7. CTweight_mat: Array vector specifying the different CT weights to test
%                  --> Ex: [1/4,1/3,1/2,2/3,3/4]
% 8. scale_mat: Array vector specifying the different 'Scale' values to test.
%               --> Ex: [1,2,3,4,5]
% 9. Ng_mat: Array vector specifying the different 'Ng' values to test.
%            --> Ex: [8,16,32,64]
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath=pwd;

% INITIALIZATION
nPatient = numel(namePT);

% COMPUTATION
fprintf('\n')
for i = 1:nPatient
    cd(pathData)
    sDataPT = load(namePT{i}); sDataPT = struct2cell(sDataPT); sDataPT = sDataPT{1};
    sDataCT = load(nameCT{i}); sDataCT = struct2cell(sDataCT); sDataCT = sDataCT{1};
    tStart = tic;
    fprintf(['\n*********************** COMPUTING TEXTURES: FUSED %s and %s ***********************'],namePT{i},nameCT{i})
    [textures] = calcPatientFusText_HN(sDataPT,sDataCT,nameROI{i},CTweight_mat,scale_mat,algo_cell,Ng_mat);
    ind = strfind(namePT{i},'_'); namePatient = namePT{i}(1:ind(2)-1);
    cd(pathText), save([namePatient,'_PTCT_',featType,'_text'],'textures')
    time = toc(tStart);
    fprintf('TOTAL TIME: %.2f seconds\n',time)
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function value = mono_total_enum ( m, n )

%*****************************************************************************80
%
%% MONO_TOTAL_ENUM enumerates monomials in M dimensions of degree equal to N.
%
%  Discussion:
%
%    For M = 3, we have the following values:
%
%    N  VALUE
%
%    0    1
%    1    3
%    2    6
%    3   10
%    4   15
%    5   21
%
%    In particular, VALUE(3,3) = 10 because we have the 10 monomials:
%
%      x^3, x^2y, x^2z, xy^2, xyz, xz^3, y^3, y^2z, yz^2, z^3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the maximum degree.
%
%    Output, integer VALUE, the number of monomials in D variables,
%    of total degree N.
%
  value = nchoosek ( n + m - 1, n );

  return
end


text:
function [index_set, abundances, reconstruct, err] = AAM_no_subset( X, L, K )
% AAM_NO_SUBSET Alternating angle minimization implementation of the MESMA
% algorithm, called by AAM
%
% This function executes the AAM algorithm on a set of libraries. For every
% pixel in x, an endmember set is constructed containing a single endmember
% from each library. The indices and abundances are returned. The algorithm
% employed is iterative alternating angle minimization.
%
% Input:  X: the N data points to unmix in d dimensions, (d,N)
%         L: spectral library, cell array of p elements of size (d,~)
%         K: Optional number of iterations. Default 3
% Output: index_set identifying the endmembers from each library, (p,N)
%         abundances with respect to these endmembers, (p,N)
%         reconstruct contains the reconstructed spectral (d,N)
%         err contains the reconstruction error (Euclidean distance) (1,N)
%
%
% Rob Heylen, 2016, University of Antwerp.

% Turn this on to activate sanity checks on the libraries and input values.
% Turn this off if you are sure there are no doubles in the libraries, and
% you do not want to check if pixels are contained in the libraries.
library_check=1;

% Initializations
if nargin==2
    K=3;
end
[d,numpx]=size(X);
p=numel(L);
for i=1:p
    N(i)=size(L{i},2);
end
flag=0;
index_set=zeros(p,numpx);
abundances=zeros(p,numpx);
F=zeros(d,p);
I=ones(p,1);
reconstruct=zeros(d,numpx);
err=zeros(1,numpx);

% Main loop over all pixels
for px=1:numpx
    x=X(:,px);

    % Check if x is a library member. If so, we can finish immediately
    if library_check
        for i=1:p
            %if sum(sum(abs(L{i}-x*ones(1,N(i))))==0)>0
            if numel(find(~sum(abs(L{i}-x*ones(1,N(i))))))>0
                I=ones(p,1);
                I(i)=find(sum(abs(L{i}-x*ones(1,N(i))))==0);
                index_set(:,px)=I;
                flag=1;
                break;
            end
        end
    end
    if flag
        flag=0;
        continue;
    end

    % Create random initial endmember set
    for i=1:p
        I(i)=ceil(rand*N(i));
        F(:,i)=L{i}(:,I(i));
    end

    % Iterate K times
    for it=1:K
        % Alternating angle optimization
        for i=1:p
            % Calculate angles
            Fi=F(:,[1:i-1 i+1:p]);       % Pivot
            Gi=[Fi x];                   % Plane through pivot and pixel
            E1=plane_project2(L{i},Fi);  % Project library onto pivot
            E2=plane_project2(L{i},Gi);  % Project library onto plane
            p1=sqrt(sum((E1-L{i}).^2));  % Distances from library to pivot
            p2=sqrt(sum((E2-L{i}).^2));  % Distances from library to plane
            ang=asin(p2./p1);            % Resulting angles

            % Find angles that should be inverted
            mask=(x-plane_project2(x,Fi))'*(E2-E1)<0;
            ang(mask)=pi-ang(mask);

            % Identify minimal angle, update index set
            [~,I(i)]=min(ang);

            % Update endmember set
            F(:,i)=L{i}(:,I(i));
        end
    end

    % Update index_set with obtained indices
    index_set(:,px)=I;
end

% Unmixing phase
E=zeros(d,p);
go=0;
for px=1:numpx
    for i=1:p
        E(:,i)=L{i}(:,index_set(i,px));
    end
    % Plug in your favorite unmixing program here

    if go==0
        [at,opt]=FCLSU_fast2(X(:,px)',E);
        abundances(:,px)=at';
        go=1;
    else
        at=FCLSU_fast2(X(:,px)',E,opt);
        abundances(:,px)=at';
    end

    % Reconstruction
    reconstruct(:,px)=E*abundances(:,px);

    % Error
    err(px)=norm(reconstruct(:,px)-X(:,px));
end

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [b,a]=potsband(fs)
%POTSBAND Design filter for 300-3400 telephone bandwidth [B,A]=(FS)
%
%Input: FS=sample frequency in Hz
%
%Output: B/A is a discrete time bandpass filter with a passband gain of 1
%
%The filter meets the specifications of G.151 for any sample frequency
%and has a gain of -3dB at the passband edges.

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

szp=[0.19892796195357i; -0.48623571568937+0.86535995266875i];
szp=[[0; -0.97247143137874] szp conj(szp)];
% s-plane zeros and poles of high pass 3'rd order chebychev2 filter with -3dB at w=1
zl=2./(1-szp*tan(300*pi/fs))-1;
al=real(poly(zl(2,:)));
bl=real(poly(zl(1,:)));
sw=[1;-1;1;-1];
bl=bl*(al*sw)/(bl*sw);
zh=2./(szp/tan(3400*pi/fs)-1)+1;
ah=real(poly(zh(2,:)));
bh=real(poly(zh(1,:)));
bh=bh*sum(ah)/sum(bh);
b=conv(bh,bl);
a=conv(ah,al);

author: jtkim-kaist

repo: VAD

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


text:
function [labels] = predictor_linear_classifier(X,B_hat,mapping)
%PREDICTOR_LINEAR_CLASSIFIER
%
%
%

[~,~,labels]  = linear_classifier('predict',X,[],[],'classification',B_hat,mapping);

end

author: epfl-lasa


text:
%% 3D Plot Testing

%% LP2
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)

%% LP2 [x0, unsolved]
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','x0',[7.3;1.16;26.163],'options',opts)
plot(Opt)

%% LP2 [x0, unsolved, multi scale]
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','x0',[7.3;1.16;26.163],'options',opts)
plot(Opt,[5 13 0 4 24 31])

%% LP2 [NO x0, unsolved]
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','options',opts)
plot(Opt)

%% LP2 [NO x0, unsolved, multi scale]
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','options',opts)
plot(Opt,[5 13 0 4 24 31])

%% LP2 [NO x0, unsolved, multi solve + multiplot]
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
%Setup Options
opts = optiset('solver','clp','display','iter');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',[1 1.1 1.2],40,'bounds',[6 1 25]',[12 3 30]','options',opts)
multisolve(Opt)
multiplot(Opt)

%% QCQP1 [-0.4004]
clc
%Objective & Constraints
H = [33 6    0;
     6  22   11.5;
     0  11.5 11];
f = [-1;-2;-3];
A = [-1 1 1; 1 -3 1];
b = [5;4];
Q = eye(3);
l = [2;2;2];
r = 0.5;
lb = [0;0;0];
ub = [40;inf;inf];
%Build & Solve
Opt = opti('H',H,'f',f,'ineq',A,b,'qc',Q,l,r,'bounds',lb,ub,'options',optiset('display','iter','solver','ipopt'))
[x,fval,exitflag,info] = solve(Opt)

plot(Opt)

%% MIQP2 [-2.75]
clc
%Objective & Constraints
H = eye(3);
f = -[2 3 1]';
A = [1 1 1;3 -2 -3; 1 -3 2];
b = [1;1;1];
%Setup Options
opts = optiset('solver','auto');
%Build & Solve
Opt = opti('hess',H,'grad',f,'ineq',A,b,'int','CIC','options',opts)
try
    [x,fval,exitflag,info] = solve(Opt)
    plot(Opt)
catch
end

%% MIQCQP2 [-0.0942]
clc
%Objective & Constraints
H = [33 6    0;
     6  22   11.5;
     0  11.5 11];
f = [-1;-2;-3];
A = [-1 1 1; 1 -3 1];
b = [5;4];
Q = eye(3);
l = [0;0;0];
r = 1;
lb = [0;0;0];
ub = [40;inf;inf];
int = 'CCI';
%Setup Options
opts = optiset('solver','auto','display','iter');
% Solve
Opt = opti('H',H,'f',f,'ineq',A,b,'qc',Q,l,r,'bounds',lb,ub,'int',int,'options',opts)
[x,fval,exitflag,info] = solve(Opt)

plot(Opt)

%% SDP4 [10.1787]
clc
clear
%Objective
f = [1 1 1];
%Linear Constraints
lb = [10;0;0];
ub = [1000;1000;1000];
%SDP Constraints [x1 1 2; 1 x2 3; 2 3 100] >= 0
C = -[0 1 2; 1 0 3; 2 3 100];
A0 = [1 0 0; 0 0 0; 0 0 0];
A1 = [0 0 0; 0 1 0; 0 0 0];
A2 = zeros(3);
sdcone = [C A0 A1 A2];
%Setup Options
dopts = dsdpset('ptol',1e-8,'rho',3,'zbar',0);
opts = optiset('solver','dsdp','display','iter','solverOpts',dopts);
% Solve
Opt = opti('f',f,'bounds',lb,ub,'sdcone',sdcone,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)

%% NLP1 Hock & Schittkowski #71
clc
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [ prod(x);
               sum(x.^2)];
nljac = @(x) [ prod(x)./x';
                2*x' ];
nlrhs = [25 40]';
nle = [1 0]'; % (>=, ==)
%Setup Options
opts = optiset('solver','ipopt','warnings','all','display','iter','derivCheck','on');
%Build & Solve
Opt = opti('obj',obj,'grad',grad,'nlmix',nlcon,nlrhs,nle,'nljac',nljac,'bounds',lb,ub,'options',opts)
x0 = [1 5 5 1]';
[x,fval,exitflag,info]= solve(Opt,x0)
plot(Opt,2)

%% NLP1 Hock & Schittkowski #71 [multiplot]
clc
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [ prod(x);
               sum(x.^2)];
nljac = @(x) [ prod(x)./x';
                2*x' ];
nlrhs = [25 40]';
nle = [1 0]'; % (>=, ==)
%Setup Options
opts = optiset('solver','ipopt','warnings','all','display','iter','derivCheck','on');
%Build & Solve
Opt = opti('obj',obj,'grad',grad,'nlmix',nlcon,nlrhs,nle,'nljac',nljac,'bounds',lb,ub,'options',opts)
x0 = [1 5 5 1]';
[x,fval,exitflag,info]= multisolve(Opt,x0)
multiplot(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y = winconv(x,varargin)
%WINCONV   Discrere time convolution of a sequence with a window.
%   y = WINCONV(X) convolves the sequence X with a rectangular window. The
%   length of the window is same as that of X.
%
%   y = WINCONV(X,WINTYPE) convolves the sequence X with the window
%   specified by WINTYPE. WINTYPE can be either 'rect', or 'hamming', or
%   'hanning', or 'bartlett', or 'blackman'. The length of the window is
%   same as that of X.
%
%   y = WINCONV(X,WINTYPE,WINAMP) convolves the sequence X with the window
%   specified by WINTYPE. The amplitude of the window is set by WINAMP. So
%   WINAMP must be a real number/vector. The length of the window is same
%   as that of X. WINAMP could be a real constant or a vector with WINLEN
%   elements.
%
%   y = WINCONV(X,WINTYPE,WINAMP,WINLEN) convolves the sequence X with the
%   window specified by WINTYPE having amplitude WINAMP and length WINLEN.
%
%   See also CONV, FFT, RECTWIN, HAMMING, HANNING, BARTLETT, BLACKMAN.
%
%   Author: Nabin Sharma
%   Date: 2009/03/15

error(nargchk(1,4,nargin,'struct'));

len = length(varargin);
switch len
    case 0
        wintype = 'rectwin';
        A = 1;
        L = length(x);
    case 1
        if ischar(varargin{1})
            wintype = lower(varargin{1});
            A = 1;
            L = length(x);
        end
    case 2
        if ischar(varargin{1}) && isreal(varargin{2})
            wintype = lower(varargin{1});
            A = varargin{2};
            L = length(x);
        end
    case 3
        if ischar(varargin{1}) && isreal(varargin{2}) &&...
                isreal(varargin{3})
            wintype = lower(varargin{1});
            A = varargin{2};
            L = varargin{3};
        end
end

% generate the window
w1 = (window(str2func(wintype),L)).'; A = A(:).';
w = A.*w1;

% perform the convolution using FFT
NFFT = 2^(nextpow2(length(x)+L));
X = fft(x,NFFT); W = fft(w,NFFT);
Y = X.*W;
y = ifft(Y,NFFT);


text:
function sr=gabreassign(s,tgrad,fgrad,a)
%GABREASSIGN  Reassign time-frequency distribution
%   Usage:  sr = gabreassign(s,tgrad,fgrad,a);
%
%   `gabreassign(s,tgrad,fgrad,a)` reassigns the values of the positive
%   time-frequency distribution *s* using the phase gradient given by *fgrad*
%   and *tgrad*. The lattice is determined by the time shift *a* and the
%   number of channels deduced from the size of *s*.
%
%   *fgrad* and *tgrad* can be obtained by the routine |gabphasegrad|.
%
%   Examples:
%   ---------
%
%   The following example demonstrates how to manually create a
%   reassigned spectrogram. An easier way is to just call |resgram|:::
%
%     % Create reassigned vector field of the bat signal.
%     a=4; M=100;
%     [tgrad, fgrad, c] = gabphasegrad('dgt',bat,'gauss',a,M);
%
%     % Perform the actual reassignment
%     sr = gabreassign(abs(c).^2,tgrad,fgrad,a);
%
%     % Display it using plotdgt
%     plotdgt(sr,a,143000,50);
%
%   See also: resgram, gabphasegrad
%
%   References: aufl95

% AUTHOR: Peter L. Søndergaard, 2008.

thisname = upper(mfilename);
complainif_notenoughargs(nargin,4,thisname);
complainif_notposint(a,'a',thisname);

% Basic checks
if any(cellfun(@(el) isempty(el) || ~isnumeric(el),{s,tgrad,fgrad}))
    error('%s: s, tgrad, fgrad must be non-empty and numeric.',...
          upper(mfilename));
end

% Check if argument sizes are consistent
if ~isequal(size(s),size(tgrad),size(fgrad))
   error('%s: s, tgrad, fgrad must all have the same size.',...
          upper(mfilename));
end

% Check if any argument is not real
if any(cellfun(@(el) ~isreal(el),{tgrad,fgrad}))
   error('%s: tgrad, fgrad must be real.',...
          upper(mfilename));
end

% if any(s<0)
%     error('%s: s must contain positive numbers only.',...
%         upper(mfilename));
% end

sr=comp_gabreassign(s,tgrad,fgrad,a);

% The following code is currently not actived. It calculates the
% reassigment using anti-aliasing, but it make very little visual
% difference, and it is slower.
  %   [M,N,W]=size(s);
  %   L=N*a;
  %   b=L/M;

  %   freqpos=fftindex(M);
  %   tgrad=bsxfun(@plus,tgrad/b,freqpos);

  %   timepos=fftindex(N);
  %   fgrad=bsxfun(@plus,fgrad/a,timepos.');

  %   tgrad=round(tgrad);
  %   fgrad=round(fgrad);

  %   tgrad=mod(tgrad,M);
  %   fgrad=mod(fgrad,N);

  %   sr=zeros(M,N,W);

  %   fk=mod(floor(tgrad),M)+1;
  %   ck=mod(ceil(tgrad),M)+1;
  %   fn=mod(floor(fgrad),N)+1;
  %   cn=mod(ceil(fgrad),N)+1;

  %   alpha = fgrad-floor(fgrad);
  %   beta  = tgrad-floor(tgrad);
  %   m1 =(1-alpha).*(1-beta).*s;
  %   m2 =(1-alpha).*beta.*s;
  %   m3 =alpha.*(1-beta).*s;
  %   m4 =alpha.*beta.*s;
  %   for ii=1:M
  %     for jj=1:N
  %       sr(fk(ii,jj),fn(ii,jj))=sr(fk(ii,jj),fn(ii,jj))+m1(ii,jj);
  %       sr(ck(ii,jj),fn(ii,jj))=sr(ck(ii,jj),fn(ii,jj))+m2(ii,jj);
  %       sr(fk(ii,jj),cn(ii,jj))=sr(fk(ii,jj),cn(ii,jj))+m3(ii,jj);
  %       sr(ck(ii,jj),cn(ii,jj))=sr(ck(ii,jj),cn(ii,jj))+m4(ii,jj);

  %     end;
  %   end;
  % end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function u = fem2d_bvp_serene ( nx, ny, a, c, f, x, y, show11 )

%*****************************************************************************80
%
%% FEM2D_BVP_SERENE solves boundary value problem on a rectangle.
%
%  Discussion:
%
%    The program uses the finite element method, with piecewise
%    serendipity basis functions to solve a 2D boundary value problem
%    over a rectangle.
%
%    The following differential equation is imposed inside the region:
%
%      - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y)
%
%    where a(x,y), c(x,y), and f(x,y) are given functions.
%
%    On the boundary, the solution is constrained to have the value 0.
%
%    The finite element method will use a regular grid of NX nodes in X, and
%    NY nodes in Y.  Both NX and NY must be odd.
%
%    The local element numbering is
%
%      3--2--1
%      |     |
%      4     8
%      |     |
%      5--6--7
%
%    The serendipity element mass matrix is a multiple of:
%
%       6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0
%      -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0
%       2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0
%      -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0
%       3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0
%      -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0
%       2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0
%      -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of X and Y grid values.
%    NX and NY must be odd and at least 3.
%
%    Input, function A(X,Y), evaluates a(x,y);
%
%    Input, function C(X,Y), evaluates c(x,y);
%
%    Input, function F(X,Y), evaluates f(x,y);
%
%    Input, real X(NX), Y(NY), the mesh points.
%
%    Input, integer SHOW11, is 1 to print out the element matrix
%    for the element in row 1, column 1.
%
%    Output, real U(MN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%

%
%  Quadrature definitions.
%
  quad_num = 3;
  abscissa(1) = -0.774596669241483377035853079956;
  abscissa(2) = 0.000000000000000000000000000000;
  abscissa(3) = 0.774596669241483377035853079956;
  weight(1) = 0.555555555555555555555555555556;
  weight(2) = 0.888888888888888888888888888889;
  weight(3) = 0.555555555555555555555555555556;
%
%  Make room for the matrix A and right hand side b.
%
  mn = fem2d_bvp_serene_node_num ( nx, ny );

  A = zeros ( mn, mn );
  b = zeros ( mn, 1 );
%
%  Compute the matrix entries by integrating over each element.
%
  ex_num = ( nx - 1 ) / 2;
  ey_num = ( ny - 1 ) / 2;

  for ey = 1 : ey_num

    s = 2 * ey - 1;
    mm = 2 * ey;
    n = 2 * ey + 1;

    ys = y(s);
    ym = y(mm);
    yn = y(n);

    yy(1) = y(n);
    yy(2) = y(n);
    yy(3) = y(n);
    yy(4) = y(mm);
    yy(5) = y(s);
    yy(6) = y(s);
    yy(7) = y(s);
    yy(8) = y(mm);

    for ex = 1 : ex_num

      w = 2 * ex - 1;
      cc = 2 * ex;
      e = 2 * ex + 1;

      xe = x(e);
      xc = x(cc);
      xw = x(w);

      xx(1) = x(e);
      xx(2) = x(cc);
      xx(3) = x(w);
      xx(4) = x(w);
      xx(5) = x(w);
      xx(6) = x(cc);
      xx(7) = x(e);
      xx(8) = x(e);
%
%  Node indices
%
%  3  2  1  wn  cn  en
%  4     8  wm      em
%  5  6  7  ws  cs  es
%
      node(1) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex + 1;
      node(2) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex;
      node(3) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex - 1;
      node(4) = ( 3 * ey - 1 ) * ey_num + 2 * ey +     ex - 1;
      node(5) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 3;
      node(6) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 2;
      node(7) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 1;
      node(8) = ( 3 * ey - 1 ) * ey_num + 2 * ey +     ex;

      if ( show11 )
        if ( ey == 1 && ex == 1 )
          ae = zeros(8,8);
          be = zeros(8,1);
        end
      end

      if ( 0 )
        fprintf ( 1, '  %2d  %2d  %2d  %2d  %2d  %2d  %2d  %2d\n', node(1:8) );
      end

      for qx = 1 : quad_num

        xq = ( ( 1.0 - abscissa(qx) ) * x(e)   ...
             + ( 1.0 + abscissa(qx) ) * x(w) ) ...
               / 2.0;

        for qy = 1 : quad_num

          yq = ( ( 1.0 - abscissa(qy) ) * y(n)   ...
               + ( 1.0 + abscissa(qy) ) * y(s) ) ...
                 / 2.0;

          wq = weight(qx) * ( x(e) - x(w) ) / 2.0 ...
             * weight(qy) * ( y(n) - y(s) ) / 2.0;

          v = basis_serene ( xq, yq, xw, ys, xe, yn, xx, yy );
          [ vx, vy ] = basisd_serene ( xq, yq, xw, ys, xe, yn, xx, yy );

          aq = a ( xq, yq );
          cq = c ( xq, yq );
          fq = f ( xq, yq );
%
%  Build the element matrix.
%
          if ( show11 )
            if ( ey == 1 && ex == 1 )
              for i = 1 : 8
                for j = 1 : 8
                  ae(i,j) = ae(i,j) + wq * ( vx(i) * aq * vx(j) ...
                                           + vy(i) * aq * vy(j) ...
                                           + v(i)  * cq * v(j) );
                end
                be(i) = be(i) + wq * ( v(i) * fq );
              end
            end
          end

          for i = 1 : 8
            ii = node(i);
            for j = 1 : 8
              jj = node(j);
              A(ii,jj) = A(ii,jj) + wq * ( vx(i) * aq * vx(j) ...
                                         + vy(i) * aq * vy(j) ...
                                         + v(i)  * cq * v(j) );
            end
            b(ii) = b(ii) + wq * ( v(i) * fq );
          end

        end
      end
%
%  Print a sample element matrix.
%
      if ( show11 )
        if ( ey == 1 && ex == 1 )
          scale = 0.5 * ae(1,3);
          fprintf ( 1, '\n' );
          fprintf ( 1, '  The Wathen elementary mass matrix:\n' );
          fprintf ( 1, '\n' );
          ae / scale
        end
      end

    end
  end
%
%  Where a node is on the boundary,
%  replace the finite element equation by a boundary condition.
%
  k = 0;

  for y = 1 : ny

    if ( mod ( y, 2 ) == 1 )
      xhi = nx;
    else
      xhi = 1 + ( nx - 1 ) / 2;
    end

    for x = 1 : 2 : xhi
      k = k + 1;
      if ( x == 1 | x == xhi | y == 1 | y == ny )
        A(k,1:mn) = 0.0;
        A(1:mn,k) = 0.0;
        A(k,k) = 1.0;
        b(k) = 0.0;
      end
    end

  end

  if ( 0 )
    spy ( A );
    pause
  end
%
%  Solve the linear system.
%
  u = A \ b;

  return
end


text:
% prtOptionsRvPlot contains available options for plotting RVs within the PRT
%
% Properties
%   nSamplesPerDim - 1x3 int array - number of samples to use for plotting
%       classifier confidence as a function of dimensionality [1D 2D 3D]
%       Default=[100 40 20]
%   colorMapFunction - a MATLAB colormap function
%
% Methods
%   Obj = prtOptionsRvPlot('paramName',paramVal,...) - Contructor
%

classdef prtOptionsRvPlot
    % Internal function
    % xxx Need Help xxx
    properties
        nSamplesPerDim = [100 40 20]; % Number of samples to use for plotting
        colorMapFunction = @(n)hot(n); % Two class colormap function handle
        nColorMapSamples = 256;
    end

    methods
        function obj = prtOptionsRvPlot(varargin)
            obj = prtUtilAssignStringValuePairs(obj,varargin{:});
        end
    end

%     methods
%         function obj = set.nSamplesPerDim(obj,val)
%             assert(length(val)==3,'largestMatrixSize must be a positive integer');
%
%             if val < 1e3
%                 warning('prt:prtOptionsComputation','Although valid, this value for largestMatrixSize is small. Consider using a larger value.')
%             end
%         end
%     end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [rootVals,exitCode]=polyRootsMultiDim(polyCoeffMats,maxDegIncreases,useMotzkinNull)
%%POLYROOTSMULTIDIM Find the roots of a system of simultaneous multivariate
%                   polynomials. Only the affine roots are found (which are
%                   generally the only ones desired), not the roots at
%                   infinity. The algorithm used transforms the
%                   multivariate root finding problem into a generalized
%                   eigenvalue problem involving the nullspace of a
%                   Macaulay matrix. Due to finite precision errors and the
%                   combinational complexity of the algorithm, it is
%                   generally best if the maximum number of variables is
%                   <=3 and the maximum degree of the polynomials is <=3.
%                   Sparser polynomials are much easier than dense ones
%                   (where all possible monomial coefficients are nonzero).
%                   When dealing with dense polynomials, it is generally
%                   best if the maximum number of variables is 2, where the
%                   maximum degree can get as high as sixth order without
%                   notable problems.
%
%INPUTS: polyCoeffMats An nX1 or 1Xn cell array containing n hypermatrices
%               of polynomial coefficients for a multivariate polynomial. A
%               hypermatrix of the coefficients for a multivariate
%               polynomial contains elements arranged such that
%               coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%               x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%               number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs). Note that this ordering is the reverse
%               that used in the 1D polyval function that is built into
%               Matlab. The number of elements for each index in coeffs is
%               the maximum order of that dimension +1.
% maxDegIncreases An optional parameter specifying the maximum number
%               of degree increases of the Motzkin matrix to use when
%               trying to solve the multivariate polynomial system. If
%               maxDegIncreases is too small, then the polynomial root
%               solver will fail. This parameter is provided as an option
%               to help identify when solving large systems would be very
%               slow. The default if this parameter is omitted or an empty
%               matrix is passed is 10*n. If one wishes there to be no
%               limit, then one can set maxDegIncreases to Inf and degree
%               increases will be allowed until Matlab runs out of memory.
% useMotzkinNull A key step in the algorithm is the determination of
%               nullspaces. By defauly, the nullspace function is used.
%               However, if desired, the Motzkin algorithm described in [1]
%               can be used. The Motzkin algorithm is generally not as
%               numerically stable, but it can be useful if one wishes to
%               step through the code and compare the results to the values
%               in [1].
%
%OUTPUTS: rootVals An nXnumSol matrix of the numSol zeros of the polynomial
%                  system found, or an empty matrix if exitCode does not
%                  equal zero.
%         exitCode A parameter indicating how the algorithm terminated.
%                  Possible values:
%                  0 The algorithm was a success; the roots were found.
%                  1 The maximum number of degree increases elapsed.
%                  2 A finite precision error occured either causing the
%                    nullity of the Macaulay matrix to change after
%                    stabilization, or to decrease when the degree was
%                    increased.
%
%This function implements Algorithm 3 of [1], an affine null space based
%root finding method for multivariate poylnomials. The polynomials in
%polyCoeffMats are converted into term matrices. The degree negative
%lexicographic ordering in the paper, for a fixed degree, coincides with
%the ordering of the terms provided by unrankTComposition with
%firstElMostSig=true. Thus, ranking and unranking of compositions is used
%to keep track of the monomials in the columns of the Macaulay matrix.
%
%The algorithm requires that a "shift function" be chosen (g(x) in [1]).
%The shift function was arbitrarily chosen to be the sum_{i=1}^n i*xi,
%where xi is the ith unknown variable. That is, the sum of all first degree
%terms multiplied by their position.
%
%The use of the Motzkin algorithm for finding the nullspace computation
%(Sections 3.2.1, 3.2.2 and 6.2.5 of [1]) allows one to avoid computation
%of degree-augmented Macaulay matrices, and provides a basis for the
%nullspace havign the nice sparse structure as used in [1]. However,
%Section 6.2.3 demonstrated that the choice of basis for the nullspace does
%not change the final solution. Thus, the option to use the (possibly more
%stable) null command (which uses singular value decomposition) is offered,
%though it does require keeping track of the entire Macaulay matrix.
%
%Though the Macaulay matrix has a sparse quasi-Toeplitz structure, the
%sparsity was not taken advantage of in the implementation of this
%algorithm.
%
%Roots at infinity are not "true" solutions to the problem; they only arise
%when the system of polynomials has been homogenized (all monomial terms in
%each equation made to have the same total degree) through the introduction
%of an additional variable. Section 2.3 of [1] discusses what they are.
%This algorithm does not find them.
%
%EXAMPLE 1:
%An example is the third-order polynomial system on page 94 of [1]. The
%equations are
%0=x1^2+5*x1*x2+4*x2*x3-10
%0=x2^3+3*x1^2*x2-12
%0=x3^3+4*x1*x2*x3-8
%Thus, the solution is
% p=zeros(4,4,4);%First Equation
% p(2+1,0+1,0+1)=1;
% p(1+1,1+1,0+1)=5;
% p(0+1,1+1,1+1)=4;
% p(0+1,0+1,0+1)=-10;
%
% q=zeros(4,4,4);%Second Equation
% q(0+1,3+1,0+1)=1;
% q(2+1,1+1,0+1)=3;
% q(0+1,0+1,0+1)=-12;
%
% k=zeros(4,4,4);%Third Equation
% k(0+1,0+1,3+1)=1;
% k(1+1,1+1,1+1)=4;
% k(0+1,0+1,0+1)=-8;
%
% polyCoeffMats=cell(3,1);
% polyCoeffMats{1}=p;
% polyCoeffMats{2}=q;
% polyCoeffMats{3}=k;
% rootVals=polyRootsMultiDim(polyCoeffMats)
%Where 18 solutions are found. Two are real; the rest complex.
%
%EXAMPLE 2:
%This is the system of equations used as an example in Section 2.1. of [1].
%0=-x1^2+2*x1*x2+x2^2+5*x1-3*x2-4
%0=x1^2+2*x1*x2+x2^2-1
%Thus, the solution is
% p=zeros(2+1,2+1);
% p(0+1,0+1)=-4;
% p(2+1,0+1)=-1;
% p(1+1,1+1)=2;
% p(0+1,2+1)=1;
% p(1+1,0+1)=5;
% p(0+1,1+1)=-3;
%
% q=zeros(3+1,3+1);
% q(0+1,0+1)=-1;
% q(2+1,0+1)=1;
% q(1+1,1+1)=2;
% q(0+1,2+1)=1;
% polyCoeffMats=cell(2,1);
% polyCoeffMats{1}=p;
% polyCoeffMats{2}=q;
% rootVals=polyRootsMultiDim(polyCoeffMats)
%In this instance, all four roots of the system are real. The solutions are
%(4,-5), (1,0), (3,-2), (0,-1).
%
%REFERENCES:
%[1] P. Dreesen, "Back to the roots: Polynomial system solving using linear
%    algebra," Ph.D. dissertation, Katholieke Universiteit Leuven, Leuven,
%    Flanders, Belgium, Sep. 2013.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The number polynomials and the number of variables.
n=length(polyCoeffMats);

if(nargin<2||isempty(maxDegIncreases))
    maxDegIncreases=10*n;
end

if(nargin<3||isempty(useMotzkinNull))
    useMotzkinNull=false;
end

%Determine the number of variables. Not all polynomials have to have all of
%the same variables, so
termMats=cell(n,1);

%Holds the degrees of each polynomial.
d=zeros(n,1);

%Convert the polynomial hypermatrices to term matrices. The order of the
%monomials in the term matrix does not matter here (though the ordering of
%the monomials when added to the Motzkin matrix is important).
for curPoly=1:n
    termMat=multiDimPolyMat2Terms(polyCoeffMats{curPoly},2,n);

    %Normalize the equation so that the highest coefficient has magnitude
    %one. This helps reduce finite precision problems.
    maxVal=max(abs(termMat(1,:)));
    termMat(1,:)=termMat(1,:)/maxVal;

    termMats{curPoly}=termMat;

    %The total degree of the current polynomial.
    d(curPoly)=max(sum(termMat(2:end,:),1));
end

%The maximum total degree of all of the polynomials.
d0=max(d);

%This keeps track of the number of columns of the Macaulay matrix  before
%the first column containing a certain degree monomial term (the
%number of monomial terms before the first monomial term of the given
%degree. The array numBeforeDeg is extended each time another degree block
%is added to the Macaulay matrix. numBeforeDeg(i+1) is the number of items
%before degree i, starting with degree 0, where numBeforeDeg(0+1)=0.
numBeforeDeg=getNumElsBeforeDeg(n,d0);

%Construct the initial Macaulay matrix of degree d0, which is the highest
%total degree of all of the polynomials going in.
[N,numBeforeDeg]=buildInitialMacaulayMat(termMats,d,d0,numBeforeDeg);
p=size(N,1);
q=size(N,2);

%Compute the inital basis for the nullspace
if(useMotzkinNull)
    Z=MotzkinMatrix(N);
else
    Z=nullspace(N);
end
N=[];%N is no longer needed.

%We keep track of the nullity to determine when the size of the nullspace
%has stabilized. Only after it has stabilized do we consider whether we can
%identify where any solutions at infinity are (so as to eliminate them).
nullity=size(Z,2);

dG=[];
dCur=d0;

nullityStabilized=false;
for curDegIncrease=1:maxDegIncreases
    %Get the new rows that would have to be augmented.
    [NRows,numBeforeDeg]=newRows4MacaulayMatrix(p,numBeforeDeg,termMats,d,dCur);
    qNew=size(NRows,2);

    %Compute the augmented nullspace for the Macaulay matrix using
    %Motzkin's algorithm (Algorithm 4). This used the block method of
    %Section 6.2.5 to expand the nullspace.
    N1=NRows(:,1:q);
    N2=NRows(:,(q+1):qNew);
    if(useMotzkinNull)
        XY=MotzkinMatrix([N1*Z,N2]);
    else
        XY=nullspace([N1*Z,N2]);
    end
    numZPrev=size(Z,2);
    X=XY(1:numZPrev,:);
    Y=XY((numZPrev+1):end,:);
    Z=[Z*X;Y];

    nullityNew=size(Z,2);

    dCur=dCur+1;
    p=size(N,1);
    q=qNew;

    %If we have found the degree for which the nullity no longer increases.
    if(nullityNew==nullity)
        nullityStabilized=true;

        %Check whether we are at dG -the degree where we can detect the
        %affine basis set (and thus know that there are ma solutions).
        [isAtdG,ma,degOfGap]=checkFordG(Z,numBeforeDeg,dCur);

        if(isAtdG)
            dG=dCur;
            break;
        end
    elseif(nullityStabilized==true||nullityNew<nullity)
        %The nullity has stabilized and then become unstabilized or the
        %nullity decreased. Neither should happen and both are indicative
        %of finite precision problems.
        rootVals=[];
        exitCode=2;
        return;
    else
        nullity=nullityNew;
    end
end

if(isempty(dG))
    rootVals=[];
    exitCode=1;%The maximum number of iterations elapsed.
    return;
end

%S1 will select all monomials up to one less than the degree of the gap. We
%will use all monomials in the nullspace up to the degree of the gap.
k=numBeforeDeg(degOfGap+2);
S1=eye(numBeforeDeg(degOfGap+1),k);

%Perform column compression on Z to get W11 using Theorem 6.9
[~,~,Q]=svd(Z(1:k,:));
W=Z*Q;
W11=W(1:k,1:ma);

Sg=constructSg(n,degOfGap,numBeforeDeg);

%The generalized eigenvalue problem to solve is S1*W11*V11*D=Sg*W11*V11
%However, eig cannot solve that as it is a rectangular problem. We use the
%method of Section 6.2.2 to transform it into a square eigenvalue problem.
[V11,~]=eig(lsqminnorm(S1*W11,Sg*W11));

%Using Corollary 6.11 to extract the actual solutions.
Ka1=W11*V11;
Ka1=bsxfun(@rdivide,Ka1,Ka1(1,:));
rootVals=Ka1(2:(n+1),:);

exitCode=0;%Successful return

end

function Sg=constructSg(n,maxDeg,numBeforeDeg)
%This function constructs Sg where g is the shift function, chosen to be
%the sum of all first order monomial terms weighted by their index.
%That is x1+2*x2+3*x3+...+n*xn
%These values are defined through Proposition 6.3. Our chose of sg was
%arbitrary.

termMats{1}=[(1:n);eye(n)];

Sg=buildInitialMacaulayMat(termMats,1,maxDeg,numBeforeDeg);
end

function [dGHasBeenReached,ma,degOfGap]=checkFordG(Z,numBeforeDeg,maxDeg)
%This implements Corollary 6.12 to determine whether we have reached a
%point where we can identify the affine basis set.

    curRank=1;%The rank of the degree-0 block.

    for curDeg=1:maxDeg
        selRows=(1:numBeforeDeg(curDeg+2));
        newRank=rank(Z(selRows,:));
        if(newRank==curRank)
            dGHasBeenReached=true;
            ma=curRank;%The number of affine bases.
            degOfGap=curDeg;
            return;
        end
        curRank=newRank;
    end
    dGHasBeenReached=false;%No gap detected.
    ma=[];
    degOfGap=[];
end

function W=MotzkinRow(b,epsVal)
%Chapter 3.2.1. A subroutine in the Motzkin method for finding the null
%space of a matrix. This implements the row-subroutine used by the
%MotzkinMatrix function.

n=length(b);
W=zeros(n,n-1);

%Explicitly zero values that are too small in magnitude.
b(abs(b)<epsVal)=0;

%Step 1: Find the rightmost pivot. This is the rightmost nonzero value.
ip=find(b,1,'last');

%If there is no nonzero element.
if(isempty(ip))
    W=eye(n,n);
    return;
end

%Step 2: Rescale b
b=b/b(ip);

%Step 3
e=zeros(n,1);
for i=(ip+1):n
    e(i)=1;
    W(:,i-1)=e;
    e(i)=0;
end

%Step 4
for i=(ip-1):-1:1
    W(i,i)=1;
    W(ip,i)=-b(i);
end
end

function H=MotzkinMatrix(A)
%Chapter 3.2.2 to find the Motzkin (canonical) nullspace of a matrix.
%An example with A=[0,2,1,0,2,0;1,0,3,2,0,1;4,0,3,2,1,0] in the paper.

m=size(A,1);

%A value for declaring elements of b zero for purposes of the Motzkin row
%algorithm.
epsVal=max(size(A))*eps(norm(A));

b=A(1,:);

H=MotzkinRow(b,epsVal);
for i=2:m
   b=A(i,:)*H;
   H=H*MotzkinRow(b,epsVal);
end
end

function [MRows,numBeforeDeg]=newRows4MacaulayMatrix(p,numBeforeDeg,termMats,d,dCur)
%To enlarge the Macaulay matrix to one degree higher, new rows must be
%appended. This computes those rows. The Macaulay matrix is defined in
%Section 5.1.2. This function only computes the rows of a particular new
%degree. When increasing the degree of the Macaulay matrix, it is not
%necessary to respect the same row ordering as in a Macaulay matrix built
%to a particular degree from scratch.

s=length(termMats);%The number of polynomials
n=size(termMats{1},1)-1;%The number of variables

%The new maximum degree if these rows are appended.
d0=dCur+1;

%Update offsets with the current degree.
numCompositionsInDeg=binomial(dCur+1+n-1,n-1);
numBeforeDeg(end+1)=numBeforeDeg(end)+numCompositionsInDeg;

[pNew,qNew]=MacaulayMatrixSize(d0,d);

%Allocate space for the new rows.
MRows=zeros(pNew-p,qNew);

curRow=1;
for i=1:s%Going through all the polynomials
    termMat=termMats{i};%The current polynomial being manipulated.
    numTerms=size(termMat,2);

    degree=d0-d(i);

    %For each value n, all possible monomials are given by all
    %compositions of the degree into n parts (as n is the number of
    %variables). Of course, these compositions are defined such that
    %the minimum value in each element is 1 and we want it to be zero,
    %so we have to adjust for that.

    numMonomials=binomial(degree+n-1,n-1);
    for j=0:(numMonomials-1)
        curMonomial=unrankTComposition(j,n,degree+n,true)-1;

        %Add the polynomial multiplied by the current monomial.
        for curTerm=1:numTerms
            monomial2Add=curMonomial+termMat(2:end,curTerm);

            %The degree of the term to add.
            deg=sum(monomial2Add);
            offset=numBeforeDeg(deg+1);

            idx=rankTComposition(monomial2Add+1,true)+offset+1;
            MRows(curRow,idx)=termMat(1,curTerm);
        end
        curRow=curRow+1;
    end
end
end

function numBeforeDeg=getNumElsBeforeDeg(n,dMax)
%A composition of n into d parts is the number of ways of splitting the
%integer n into the sum of d nonzero values. Compositions can also be
%offset to split n into d values, which are >=0. The lexicographic rank of
%a composition (starting from 0) is the offset within graded lexicographic
%ordering of a value with a particular fixed degree. We are going to have
%to keep track of these offsets to be able to index vectors.

numBeforeDeg=zeros(dMax+1,1);
numBeforeDeg(0+1)=0;
numBeforeDeg(1+1)=1;
for degree=1:dMax
    %The number of compositions of degree into n parts, where the parts can
    %be zero. We are including d0 as the function checkFordG will need it
    %even though this function only needs it up to degree d0-1.
    numCompositionsInDeg=binomial(degree+n-1,n-1);

    numBeforeDeg(degree+2)=numBeforeDeg(degree+1)+numCompositionsInDeg;
end

end

function [M,numBeforeDeg]=buildInitialMacaulayMat(termMats,d,d0,numBeforeDeg)
%%BUILDINITIASLMACAULAYMAT Construct the minimum-sized Macaulay matrix for
%             the system given. The degree of the matrix is the maximum
%             total degree of all of the polynomials. The Macaulay matrix
%             is defined in Section 5.1.2. This function is also used as a
%             convenience function to construct Sg.

s=length(termMats);%The number of polynomials
n=size(termMats{1},1)-1;%The number of variables

%Initially a Macaulay matrix of degree d0 is constructed. The number of
%rows p and columns q of the matrix are given
[p,q]=MacaulayMatrixSize(d0,d);

%Allocate space for the d0th order Macaulay matrix.
M=zeros(p,q);

%Now, fill in the Macaulay matrix of order d0. Here, each polynomial is
%multiplied by all monomials of degree <=d0-d(i) and added. This includes
%the zeroth-degree, which contains
curRow=1;
for i=1:s%Going through all the polynomials
    termMat=termMats{i};%The current polynomial being manipulated.

    %First, add the zeroth-order term (the polynomial itself.
    numTerms=size(termMat,2);
    for curTerm=1:numTerms
        %The degree of the term to add.
        deg=sum(termMat(2:end,curTerm));
        offset=numBeforeDeg(deg+1);

        idx=rankTComposition(termMat(2:end,curTerm)+1,true)+offset+1;
        M(curRow,idx)=termMat(1,curTerm);
    end
    curRow=curRow+1;

    %Now, add all other terms associated with this polynomial.
    for degree=1:(d0-d(i))%All the orders for the current polynomial.
        %For each order, all possible monomials are given by all
        %compositions of degree into n parts (as n is the number of
        %variables). Of course, these compositions are defined such that
        %the minimum value in each element is 1 and we want it to be zero,
        %so we have to adjust for that

        numMonomials=binomial(degree+n-1,n-1);
        for j=0:(numMonomials-1)
            curMonomial=unrankTComposition(j,n,degree+n,true)-1;

            %Add the polynomial multiplied by the current monomial.
            for curTerm=1:numTerms
                monomial2Add=curMonomial+termMat(2:end,curTerm);

                %The degree of the term to add.
                deg=sum(monomial2Add);
                offset=numBeforeDeg(deg+1);

                idx=rankTComposition(monomial2Add+1,true)+offset+1;
                M(curRow,idx)=termMat(1,curTerm);
            end
            curRow=curRow+1;
        end
    end
end
end

function [p,q]=MacaulayMatrixSize(d0,d)
%%%MACAULAYMATRIXSIZE The size of a d0th order Macaulay matrix when the
%                   polynomials have degree d(i). This implements Lemma 5.8
%                   in Section 5.2.1 of [1].

%The number of polynomials.
n=length(d);

q=binomial(n+d0,d0);

p=0;
for i=1:n
   p=p+binomial(n+d0-d(i),d0-d(i));
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [procMask3M, planC] = post_process_larynx(strNum,paramS,planC)
% AI 03/26/2020
% Morphological post-processing for auto-segmentation of larynx.
%--------------------------------------------------------------------------

%Get auto-segemented mask
[label3M, planC] = getStrMask(strNum,planC);
slicesV = find(squeeze(sum(sum(double(label3M)))>0));

maskSiz = size(label3M,1);
scale = 512/maskSiz;

%Post-process
if ~isempty(slicesV)

    filtSize = floor(5/scale);
    filtSize = max(1,2*floor(filtSize/2)-1); %Nearest odd val.
    conn = 26;

    strMask3M = zeros(size(label3M,1),size(label3M,1),length(slicesV));
    sliceLabels3M = label3M(:,:,slicesV);

    %Fill holes
    sliceLabels3M = imclose(sliceLabels3M,strel('sphere',floor(3/scale)));

    %Remove islands
    for s = 1:size(sliceLabels3M,3)
        slcMask = sliceLabels3M(:,:,s);
        slcMask = bwareaopen(slcMask,floor(10/scale^2),8);
        strMaskM = zeros(size(slcMask));
        connCompS = bwconncomp(slcMask,8);
        ccSiz = cellfun(@numel,[connCompS.PixelIdxList]);
        sel = ccSiz==max(ccSiz);
        if ~ (isempty(slcMask(sel)) | max(ccSiz)< floor(20/scale^2))
            idx = connCompS.PixelIdxList{sel};
            strMaskM(idx) = 1;
        end
        sliceLabels3M(:,:,s) = strMaskM;
    end

    %Retain largest connected component
    connCompS = bwconncomp(sliceLabels3M,conn);
    ccSiz = cellfun(@numel,[connCompS.PixelIdxList]);
    sel = ccSiz==max(ccSiz);
    if ~ (isempty(sliceLabels3M(sel)) | max(ccSiz)< floor(1000/scale^2))
        idx = connCompS.PixelIdxList{sel};
        strMask3M(idx) = 1;
    end

    %Smooth
    if size(strMask3M,3)>1
        smoothedlabel3M = smooth3(double(strMask3M),'box',filtSize);
        strMask3M = smoothedlabel3M > 0.5;
    end

    label3M(:,:,slicesV) = strMask3M;

end

procMask3M = label3M;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function A = convertRadToDegrees(InA)
    A = InA*180/pi;
end


text:
function rotClusters = cluster_views(model_3d,inds,nClusters,vis)

    if(nargin<4)
        vis=0;
    end
    rots = model_3d.rots(inds);
    rots = cellfun(@(x)(x(:)'),rots,'UniformOutput',false);
    rots = vertcat(rots{:});
    T = clusterdata(rots,'distance',@riemannian_dist,'linkage','average','maxclust',nClusters);
    nClusters = length(unique(T));
    rotClusters = cell(nClusters,1);
    for i=1:nClusters
        rotClusters{i} = inds(T==i);
    end

    % Visualize
    if(vis)
        montage.W = 3000;
        montage.H = 3000;
        montage.sc = 0.5;
        for i=1:nClusters
            ImAll = vis_viewclusters(model_3d, rotClusters{i}, montage);
            figure;imshow(ImAll{1});
        end
    end

end

function d = riemannian_dist(R,Rmats)
    N = size(Rmats,1);
    R = reshape(R',[3 3]);
    Rmats = reshape(Rmats',[3 3 N]);
    d = zeros(N,1);
    parfor i=1:N
        d(i) = norm(log(eig(Rmats(:,:,i)'*R)));
    end
end

function ImAll = vis_viewclusters(model_3d,inds, montage)
    globals;
    pasdir = PASCAL_DIR;
    I = cell(length(inds),1);
    for i=1:length(inds)
        im = imread([pasdir model_3d.voc_image_id{inds(i)} '.jpg']);
        if(model_3d.flip(inds(i)))
            im(:,:,1) = fliplr(im(:,:,1));im(:,:,2) = fliplr(im(:,:,2));im(:,:,3) = fliplr(im(:,:,3));
        end
        imbbox   = model_3d.bbox(inds(i),:);
        I{i} = imcrop(im,imbbox);
    end
    [~,ImAll] = makeMontage(I,montage.W,montage.H,montage.sc);
end

repo: CategoryShapes


text:
function [Sl1, Sl2, Jac] = deformation_field_3d_analysis(Tr, pixel_resolution)
% calculates sliding amount jacobian and other parameters of transformation
% Tr should be given in physical units (mm)!

dT = cell(3,3);
sz = size(Tr);
[m1,m2,m3] = ndgrid(1:sz(1), 1:sz(2), 1:sz(3));
mm{1} = m1;
mm{2} = m2;
mm{3} = m3;
for i1 = 1 : 3
    for i2 = 1: 3
        tmp = Tr(:,:,:,i1)+mm{i1}*pixel_resolution(i1);
        dT{i1,i2} = DGradient(tmp, [], i2) / pixel_resolution(i2);
    end
end

% EV = zeros([3, sz(1), sz(2), sz(3)]);
Sl1 = zeros([sz(1), sz(2), sz(3)]);
Sl2 = zeros([sz(1), sz(2), sz(3)]);
Jac = zeros([sz(1), sz(2), sz(3)]);
for i3 = 1 : sz(3)
for i2 = 1 : sz(2)
for i1 = 1 : sz(1)

    mat = [dT{1,1}(i1,i2,i3), dT{1,2}(i1,i2,i3), dT{1,3}(i1,i2,i3);...
           dT{2,1}(i1,i2,i3), dT{2,2}(i1,i2,i3), dT{2,3}(i1,i2,i3);...
           dT{3,1}(i1,i2,i3), dT{3,2}(i1,i2,i3), dT{3,3}(i1,i2,i3); ];

    mat_sym = mat'*mat;

    v = eig(mat_sym);
    v = sort(v);
    v = sqrt(abs(v));
    Sl1(i1,i2,i3) = (v(3)-v(1))/2;
    Sl2(i1,i2,i3) = v(3)/(v(2) + v(1));
    Jac(i1,i2,i3) = det(mat);
end
end
end

end

repo: pTVreg


text:
% HIV system
% System identification: (Delay)DMDc

clear all, close all, clc
figpath = '../FIGURES/HIV/'; mkdir(figpath)
datapath = '../DATA/HIV/'; mkdir(datapath)
addpath('../utils');

SystemModel = 'HIV';

% Model type selected based on DATA_ENSEMBLE and  Ndelay
%% Generate Data
InputSignalType = 'prbs';
Nvar = 5;
Ndelay = 1;  % Choose 1 for DMDc, and 10 for DelayDMDc in paper
DATA_ENSEMBLE = 0; % Choose 0 to reproduce results in paper

%%
if DATA_ENSEMBLE == 1 % ONLY FOR Ndelay = 1 implemented
    ModelName = 'DMDc';
    getTrainingData_Ensemble

    %% Reshape
    X = x(1:end-1,:,:);         % Array of snapshots
    Xp = x(2:end,:,:);          % time-shifted Array
    U = u(1:end-1,:,:);         % Array of inputs
    M = size(X,1);
    n = size(X,2);

    % Reorganize data from array to time-state matrix
    X_tmp = zeros(Nvar,Nic*M);
    Xp_tmp = zeros(Nvar,Nic*M);
    U_tmp = zeros(1,Nic*M);
    for iIC = 1:Nic
        for i = 1:Nvar
            X_tmp(i,(iIC-1)*M+1:iIC*M) = X(:,i,iIC);
            Xp_tmp(i,(iIC-1)*M+1:iIC*M) = Xp(:,i,iIC);
        end
        U_tmp(1,(iIC-1)*M+1:iIC*M) = U(:,1,iIC);
    end
    X = X_tmp; Xp = Xp_tmp; U = U_tmp;
    clear X_tmp Xp_tmp U_tmp
    M = size(X,2);

    %% Construct data matrices // Mean-correction
    xmean = zeros(1,5); %mean(X,2)';%xref;
    X   = X - repmat(xmean',[1 M]);
    Xp  = Xp - repmat(xmean',[1 M]);

    r1 = size(X,1); r2 = size(Xp,1);
    [sysmodel_DMDc,Psi,Psi_p] = DelayDMDc_MV(zeros(5,1),U,r1,r2,dt,size(X,1),size(U,1),2,X,Xp);

    Nt = length(t)-1;

    %% Validation over training phase
    xDMDc_Valid = zeros(size(X));
    for iIC = 1:Nic
        [xDMDc,~] = lsim(sysmodel_DMDc,squeeze(u(1:Nt,1,iIC)),t(1:Nt),x0_ensemble(iIC,:));
        xDMDc = xDMDc + repmat(xmean,[Nt 1]);
        xDMDc_Valid(1:5,(iIC-1)*Nt+1:iIC*Nt) = xDMDc';
        disp(['PROGRESS: ',num2str(100*iIC/Nic),'%'])
    end
    %%
    figure, hold on, box on
    plot(X','-k','LineWidth',2)
    plot(xDMDc_Valid','--','LineWidth',2)
else

    getTrainingData
    %% DMDc: B = unknown  and with time delay coordinates
    Ndelay = 1; % 1 for DMDc, 10 for DelayDMDc in paper
    if Ndelay == 1
        ModelName = 'DMDc';
    elseif Ndelay>1
        ModelName = 'DelayDMDc';
    end

    % Construct data matrices
    Hu = getHankelMatrix_MV(u,Ndelay);
    xmean = xref; %zeros(1,Nvar);%xref;%mean(x);%xref; %mean(x);
    X   = x - repmat(xmean,[T 1]);
    Hx  = getHankelMatrix_MV(X,Ndelay);
    numOutputs = size(Hx,1); numInputs = size(Hu,1); numVar = 5;
    r1 = size(Hx,1); r2 = size(Hx,1);
    [sysmodel_DMDc,U,Up] = DelayDMDc_MV(Hx,Hu,size(Hx,1),size(Hx,1),dt,size(Hx,1),size(Hu,1),2);

    Nt = length(t)-Ndelay+1;
    %% Validation over training phase
    [xDMDc,~] = lsim(sysmodel_DMDc,Hu',tspan(1:Nt),Hx(:,1));
    xDMDc = xDMDc(:,end-Nvar+1:end);
    xDMDc = xDMDc + repmat(xmean,[Nt 1]);

    %% Show prediction over training stage
    clear ph
    figure,box on,
    ccolors = get(gca,'colororder');
    ccolors_valid = [ccolors(1,:)-[0 0.2 0.2];
        ccolors(2,:)-[0.1 0.2 0.09];
        ccolors(3,:)-[0.1 0.2 0.09];
        ccolors(4,:)-[0.1 0.1 0.2];
        ccolors(5,:)-[0.1 0.2 0.09]];
    for i = 1:Nvar
        ph(i) = semilogy(tspan,x(:,i),'-','Color',ccolors(i,:),'LineWidth',1); hold on
    end
    for i = 1:Nvar
        ph(Nvar+i) = semilogy(tspan(Ndelay:Nt+Ndelay-1),xDMDc(:,i),'--','Color',ccolors_valid(i,:),'LineWidth',2);
    end
    xlabel('Time')
    ylabel('xi')
    legend(ph([1,6]),'True',ModelName)
    set(gca,'LineWidth',1, 'FontSize',14)
    set(gcf,'Position',[100 100 300 200])
    set(gcf,'PaperPositionMode','auto')
    print('-depsc2', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'.eps']);

end
%% Save Data & Model
Model.name = ModelName;
Model.sys = sysmodel_DMDc;
Model.Ndelay = Ndelay;
Model.dt = dt;
Model.xmean = xmean;
save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'.mat']),'Model')

author: eurika-kaiser

repo: SINDY-MPC