text:
function info = fcnTest(varargin)

run matconvnet/matlab/vl_setupnn ;
addpath matconvnet/examples ;

% experiment and data paths
opts.expDir = 'data/fcn32s-voc11' ;
opts.dataDir = 'data/voc11' ;
opts.modelPath = 'data/fcn32-voc11/net-epoch-50.mat' ;
opts.modelFamily = 'matconvnet' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

% experiment setup
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat') ;
opts.vocEdition = '11' ;
opts.vocAdditionalSegmentations = true ;
opts.vocAdditionalSegmentationsMergeMode = 2 ;
opts.gpus = [] ;
opts = vl_argparse(opts, varargin) ;

resPath = fullfile(opts.expDir, 'results.mat') ;
if exist(resPath)
  info = load(resPath) ;
  return ;
end

if ~isempty(opts.gpus)
  gpuDevice(opts.gpus(1))
end

% -------------------------------------------------------------------------
% Setup data
% -------------------------------------------------------------------------

% Get PASCAL VOC 11/12 segmentation dataset plus Berkeley's additional
% segmentations
if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = vocSetup('dataDir', opts.dataDir, ...
    'edition', opts.vocEdition, ...
    'includeTest', false, ...
    'includeSegmentation', true, ...
    'includeDetection', false) ;
  if opts.vocAdditionalSegmentations
    imdb = vocSetupAdditionalSegmentations(...
      imdb, ...
      'dataDir', opts.dataDir, ...
      'mergeMode', opts.vocAdditionalSegmentationsMergeMode) ;
  end
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% Get validation subset
val = find(imdb.images.set == 2 & imdb.images.segmentation) ;

% Compare the validation set to the one used in the FCN paper
% valNames = sort(imdb.images.name(val)') ;
% valNames = textread('data/seg11valid.txt', '%s') ;
% valNames_ = textread('data/seg12valid-tvg.txt', '%s') ;
% assert(isequal(valNames, valNames_)) ;

% -------------------------------------------------------------------------
% Setup model
% -------------------------------------------------------------------------

switch opts.modelFamily
  case 'matconvnet'
    net = load(opts.modelPath) ;
    net = dagnn.DagNN.loadobj(net.net) ;
    net.mode = 'test' ;
    for name = {'objective', 'accuracy'}
      net.removeLayer(name) ;
    end
    net.meta.normalization.averageImage = reshape(net.meta.normalization.rgbMean,1,1,3) ;
    predVar = net.getVarIndex('prediction') ;
    inputVar = 'input' ;
    imageNeedsToBeMultiple = true ;

  case 'ModelZoo'
    net = dagnn.DagNN.loadobj(load(opts.modelPath)) ;
    net.mode = 'test' ;
    predVar = net.getVarIndex('upscore') ;
    inputVar = 'data' ;
    imageNeedsToBeMultiple = false ;

  case 'TVG'
    net = dagnn.DagNN.loadobj(load(opts.modelPath)) ;
    net.mode = 'test' ;
    predVar = net.getVarIndex('coarse') ;
    inputVar = 'data' ;
    imageNeedsToBeMultiple = false ;
end

if ~isempty(opts.gpus)
  gpuDevice(opts.gpus(1)) ;
  net.move('gpu') ;
end
net.mode = 'test' ;

% -------------------------------------------------------------------------
% Train
% -------------------------------------------------------------------------

numGpus = 0 ;
confusion = zeros(21) ;

for i = 1:numel(val)
  imId = val(i) ;
  name = imdb.images.name{imId} ;
  rgbPath = sprintf(imdb.paths.image, name) ;
  labelsPath = sprintf(imdb.paths.classSegmentation, name) ;

  % Load an image and gt segmentation
  rgb = vl_imreadjpeg({rgbPath}) ;
  rgb = rgb{1} ;
  anno = imread(labelsPath) ;
  lb = single(anno) ;
  lb = mod(lb + 1, 256) ; % 0 = ignore, 1 = bkg

  % Subtract the mean (color)
  im = bsxfun(@minus, single(rgb), net.meta.normalization.averageImage) ;

  % Soome networks requires the image to be a multiple of 32 pixels
  if imageNeedsToBeMultiple
    sz = [size(im,1), size(im,2)] ;
    sz_ = round(sz / 32)*32 ;
    im_ = imresize(im, sz_) ;
  else
    im_ = im ;
  end

  if ~isempty(opts.gpus)
    im_ = gpuArray(im_) ;
  end

  net.eval({inputVar, im_}) ;
  scores_ = gather(net.vars(predVar).value) ;
  [~,pred_] = max(scores_,[],3) ;

  if imageNeedsToBeMultiple
    pred = imresize(pred_, sz, 'method', 'nearest') ;
  else
    pred = pred_ ;
  end

  % Accumulate errors
  ok = lb > 0 ;
  confusion = confusion + accumarray([lb(ok),pred(ok)],1,[21 21]) ;

  % Plots
  if mod(i - 1,30) == 0 || i == numel(val)
    clear info ;
    [info.iu, info.miu, info.pacc, info.macc] = getAccuracies(confusion) ;
    fprintf('IU ') ;
    fprintf('%4.1f ', 100 * info.iu) ;
    fprintf('\n meanIU: %5.2f pixelAcc: %5.2f, meanAcc: %5.2f\n', ...
            100*info.miu, 100*info.pacc, 100*info.macc) ;

    figure(1) ; clf;
    imagesc(normalizeConfusion(confusion)) ;
    axis image ; set(gca,'ydir','normal') ;
    colormap(jet) ;
    drawnow ;

    % Print segmentation
    figure(100) ;clf ;
    displayImage(rgb/255, lb, pred) ;
    drawnow ;

    % Save segmentation
    imPath = fullfile(opts.expDir, [name '.png']) ;
    imwrite(pred,labelColors(),imPath,'png');
  end
end

% Save results
save(resPath, '-struct', 'info') ;

% -------------------------------------------------------------------------
function nconfusion = normalizeConfusion(confusion)
% -------------------------------------------------------------------------
% normalize confusion by row (each row contains a gt label)
nconfusion = bsxfun(@rdivide, double(confusion), double(sum(confusion,2))) ;

% -------------------------------------------------------------------------
function [IU, meanIU, pixelAccuracy, meanAccuracy] = getAccuracies(confusion)
% -------------------------------------------------------------------------
pos = sum(confusion,2) ;
res = sum(confusion,1)' ;
tp = diag(confusion) ;
IU = tp ./ max(1, pos + res - tp) ;
meanIU = mean(IU) ;
pixelAccuracy = sum(tp) / max(1,sum(confusion(:))) ;
meanAccuracy = mean(tp ./ max(1, pos)) ;

% -------------------------------------------------------------------------
function displayImage(im, lb, pred)
% -------------------------------------------------------------------------
subplot(2,2,1) ;
image(im) ;
axis image ;
title('source image') ;

subplot(2,2,2) ;
image(uint8(lb-1)) ;
axis image ;
title('ground truth')

cmap = labelColors() ;
subplot(2,2,3) ;
image(uint8(pred-1)) ;
axis image ;
title('predicted') ;

colormap(cmap) ;

% -------------------------------------------------------------------------
function cmap = labelColors()
% -------------------------------------------------------------------------
N=21;
cmap = zeros(N,3);
for i=1:N
  id = i-1; r=0;g=0;b=0;
  for j=0:7
    r = bitor(r, bitshift(bitget(id,1),7 - j));
    g = bitor(g, bitshift(bitget(id,2),7 - j));
    b = bitor(b, bitshift(bitget(id,3),7 - j));
    id = bitshift(id,-3);
  end
  cmap(i,1)=r; cmap(i,2)=g; cmap(i,3)=b;
end
cmap = cmap / 255;

author: vlfeat

repo: matconvnet-fcn


text:
function [f, logscale] = fwd(engine, fpast, ev, t)
% Forwards pass.

f.obslik = mk_hmm_obs_lik_vec(engine, ev);
transmat = engine.transprob;
f.past_alpha = fpast.alpha;
if engine.maximize
  Q = length(fpast.alpha);
  A = repmat(fpast.alpha, [1 Q]);
  m = max(transmat .* A, [], 1);
  [f.alpha, scale] = normalise(m(:) .* f.obslik);
else
  [f.alpha, scale] = normalise((transmat' * fpast.alpha) .* f.obslik);
end
logscale = log(scale);
%f.xi = normalise((fpast.alpha * obslik') .* transmat); % t-1,t
f.t = t;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function erase_string(str)
% Backspaces the number of chars of the length of str

    fprintf(1, repmat('\b', 1, length(str)));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

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

% Last Modified by GUIDE v2.5 20-Jul-2011 14:25:09

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

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

% BrowserCurationGUI(browserObj,curationCfg)
% browserObj: handle to OlyDat.Browser object
% curationCfg: struct array providing curation fields. curationCfg
% configures the ith column in the curation table. Note that the 1st column
% is reserved for the "excluded" column. curationCfg fields:
%   * Name: curation field name
%   * TableName: name for column header in table
%   * TableColWidth: width in pixels for col. 0 indicates auto-size.
function BrowserCurationGUI_OpeningFcn(hObject, eventdata, handles, varargin)

handles.output = hObject;
handles.browser = varargin{1};
handles.curationCfg = varargin{2};
handles.selectedCell = [];

% table metadata
handles.tableMetadataExcludedColIdx = 1;
handles.tableMetadataCurationField2Col = struct();
handles.tableMetadataCol2Field = {''}; % first col is empty for excluded col

NcurationCfg = numel(handles.curationCfg);
for cfgIdx = 1:NcurationCfg
    name = handles.curationCfg(cfgIdx).Name;
    handles.tableMetadataCurationField2Col.(name) = cfgIdx+1;
    handles.tableMetadataCol2Field{cfgIdx+1,1} = name;
end

% configure table based on tableMetadata
Ncol = NcurationCfg+1;
tfColEditable = false(1,Ncol);
tfColEditable(handles.tableMetadataExcludedColIdx) = true; % excluded col is editable; rest are not
colFmt = cell(1,Ncol);
colFmt{1} = 'logical'; % unused, see below
colNames = [{'Excluded'} {handles.curationCfg.TableName}];
colWidths = [{60} {handles.curationCfg.TableColWidth}];
for colIdx = 1:Ncol
    if colWidths{colIdx}==0
        colWidths{colIdx} = 'auto';
    end
end

% (*) AL 1/5/2012: including ColumnFormat in the set leads to a SEGV on
% maci64. This is consistent/reproduceable. However simply leaving out the
% ColumnFormat set seems to work fine. The ColumnFormat property will not
% have the right size (eg will not have Ncol number of elements), but
% in this situation MATLAB appears to just use auto-formatting for all
% columns and for now that is working fine.
set(handles.tblCuration,'ColumnEditable',tfColEditable,... % 'ColumnFormat',colFmt,... (*) see above
                        'ColumnName',colNames,...
                        'ColumnWidth',colWidths,...
                        'Data',cell(1,Ncol));

% context menu
uh = uicontextmenu;
uimenu(uh,'label','Enter Curation Info','Callback',@ctxtMnuEnterCurationInfo);
uimenu(uh,'label','Clear Curation Info','Callback',@ctxtMnuClearCurationInfo);
uimenu(uh,'label','Sort Ascending','Callback',@ctxtMnuSortAscending,'Separator','on');
uimenu(uh,'label','Sort Descending','Callback',@ctxtMnuSortDescending);
uimenu(uh,'label','Bring Excluded/Marked to Top','Callback',@ctxtMnuBringExcludedMarkedToTop);
set(handles.tblCuration,'UIContextMenu',uh);

% handles.tableMetadataCurationField2Col = struct('experiment_id',2,'experiment_name',3,'manual_pf',4,
% 'notes_curation',5,'flag_redo',6,'flag_review',7,'manual_curation_date',8,'manual_curator',9);
% handles.tableMetadataCol2Field = {'';'experiment_id';'experiment_name';'manual_pf';'notes_curation';
%     'flag_redo';'flag_review';'manual_curation_date';'manual_curator'};

guidata(hObject,handles);

% This gui is custom-resizeable
set(handles.uipanel1,'Units','Normalized');

% Explicitly call resizefcn b/c the orig layout is not precisely consistent
% with what is done in resize.
uipanel1_ResizeFcn(handles.uipanel1,[],handles);

centerfig(hObject);

function uipanel1_ResizeFcn(hObject,~,handles)
TBLCURATION_Y_TOP_SPACER = 8.0;
TBLCURATION_Y_BOTTOM_SPACER = 10.0;
TBLCURATION_X_SPACER = 8.0;
TBL_CURATION_BETWEEN_BUTTONS_X_SPACER = 3.0;
TBL_CURATION_AFTER_BUTTONS_X_SPACER = 11.0;

% get panel width/height in pixels
origUnits = get(hObject,'Units');
set(hObject,'Units','pixels');
pnlPos = get(hObject,'Position');
set(hObject,'Units',origUnits);
pnlWidth = pnlPos(3);
pnlHeight = pnlPos(4);

% get top posn of unexclude all button
pbUnexcludeAllPos = get(handles.pbUnexcludeAll,'Position');
pbTopPosInPixels = pbUnexcludeAllPos(2)+pbUnexcludeAllPos(4);

% set tbl position
tblPos(1) = TBLCURATION_X_SPACER;
tblPos(2) = pbTopPosInPixels + TBLCURATION_Y_BOTTOM_SPACER;
tblPos(3) = pnlWidth-2.0*TBLCURATION_X_SPACER;
tblPos(4) = pnlHeight-tblPos(2)-TBLCURATION_Y_TOP_SPACER;
set(handles.tblCuration,'Position',tblPos);

% set posn of UnexcludeAll, MarkAllCurated buttons
pbMarkAllUncuratedPos = get(handles.pbMarkAllUncuratedAsPass,'Position');
pbUnexcludeAllXPos = pnlWidth - pbUnexcludeAllPos(3) - pbMarkAllUncuratedPos(3) ...
    - TBL_CURATION_BETWEEN_BUTTONS_X_SPACER - TBL_CURATION_AFTER_BUTTONS_X_SPACER;
pbUnexcludeAllPos(1) = pbUnexcludeAllXPos;
pbMarkAllUncuratedXPos = pnlWidth - pbMarkAllUncuratedPos(3) - TBL_CURATION_AFTER_BUTTONS_X_SPACER;
pbMarkAllUncuratedPos(1) = pbMarkAllUncuratedXPos;
set(handles.pbUnexcludeAll,'Position',pbUnexcludeAllPos);
set(handles.pbMarkAllUncuratedAsPass,'Position',pbMarkAllUncuratedPos);

function figure1_ResizeFcn(~,~,~)
%none

function varargout = BrowserCurationGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% A note on exp selection. When rows of the CurationGUI are selected, the
% corresponding experiments are selected in the Browser. However, when
% experiments are selected in the Browser, the corresponding rows of the
% curationGUI are NOT selected. The primary reason for this is that the
% native ML tables do not permit programmatically setting their selection.

function tblCuration_CellSelectionCallback(hObject,eventdata,handles)
if ~isempty(eventdata.Indices)
    handles.selectedCell = eventdata.Indices;
    guidata(hObject,handles);
    zlclSelectExperimentsInBrowserBasedOnSelectedCellCache(handles);
else
    handles.selectedCell = [];
    guidata(hObject,handles);
end

function tblCuration_CellEditCallback(hObject, eventdata, handles) %#ok<*INUSL,*DEFNU>
% If there are multiple selections in the table, the upper-left-most will be used
rowIdx = eventdata.Indices(1,1);
colIdx = eventdata.Indices(1,2);
data = get(hObject,'Data');

switch colIdx
    case handles.tableMetadataExcludedColIdx
        % 'excluded' column

        % Make sure this experiment (row) is selected. This may be
        % unnecessary, ie if tblCuration_CellSelectionCallback is always
        % called first.
        expID = data{rowIdx,handles.tableMetadataCurationField2Col.experiment_id};
        handles.browser.selectExperiment(expID);

        tf = eventdata.NewData;
        handles.browser.setRemoveExp(tf,false);

        parentFigure = ancestor(hObject,'figure');
        uistack(parentFigure,'top');
    otherwise
        assert(false,'This no longer ever occurs');
end

% This is unused
% function s = zlclRow2CurationInfoStruct(row,handles)
% flds = handles.tableMetadataCol2Field;
% Nflds = numel(flds);
% assert(numel(row)==Nflds);
% s = struct();
% for c = 1:Nflds
%     if ~isempty(flds{c}) % excluded col will have an empty field
%         val = row{c};
%         if ischar(val)
%             % convert empty option ' ' to ''
%             val = strtrim(val);
%         end
%         s.(flds{c}) = val;
%     end
% end

%%% ACTIONS/HELPERS
function zlclSelectExperimentsInBrowserBasedOnSelectedCellCache(handles)
rowIdxs = handles.selectedCell(:,1);
data = get(handles.tblCuration,'Data');
eids = data(rowIdxs,handles.tableMetadataCurationField2Col.experiment_id);
handles.browser.selectExperiment(cell2mat(eids));

function zlclDisplaySelectCellMsgBox(msg)
msgbox(msg,'No cell selected');

function zlclActionEnterCurationMarks(handles)
if isempty(handles.selectedCell)
    zlclDisplaySelectCellMsgBox('Please select a cell to specify a row.');
else
    % select experiments in the Browser and open CurationInfoEntry.
    zlclSelectExperimentsInBrowserBasedOnSelectedCellCache(handles);
    handles.browser.curationOpenCurationInfoEntry();
end

function zlclActionClearCurationMarks(handles)
if isempty(handles.selectedCell)
    zlclDisplaySelectCellMsgBox('Please select a cell to specify a row.');
else
    zlclSelectExperimentsInBrowserBasedOnSelectedCellCache(handles);
    handles.browser.curationClearCurationInfoForCurrentExp();
end

function zlclSortByData(tableH,sortData,dir)
data = get(tableH,'Data');
assert(size(data,1)==numel(sortData(:)));

[~,idx] = sort(sortData(:));
switch dir
    case 'ascend'
        %none
    case 'descend'
        idx = idx(end:-1:1);
end

data = data(idx,:);
set(tableH,'Data',data);

function zlclSortHelper(handles,direction)
if isempty(handles.selectedCell)
    zlclDisplaySelectCellMsgBox('Select a cell to specify a column for sorting.');
else
    data = get(handles.tblCuration,'Data');
    colIdx = handles.selectedCell(1,2); % use column of first selected cell
    colData = data(:,colIdx);
    if isnumeric(colData{1}) || islogical(colData{1})
        colData = cell2mat(colData);
    end
    zlclSortByData(handles.tblCuration,colData,direction);

    % After setting the table data to a new value, MATLAB de-selects any
    % selected cells. Ideally, would like to have the original cell (in
    % terms of coordinates) selected, so that the user may eg 'sort
    % ascending', then immediately 'sort descending' if they got the wrong
    % direction. Unfortunately, this is not possible with the ML uitable
    % without dropping into Java; this is a potentialy maintenance
    % dependency so for now we don't do it. To sort one way and then the
    % other, the user will have to re-select a cell in the appropriate
    % column.
end

% At the moment, this checks for marking by checking the manual_curator
% field (hardcoded).
function zlclBringExcludedMarkedToTop(handles)
data = get(handles.tblCuration,'Data');
Nrows = size(data,1);

excludedColIdx = handles.tableMetadataExcludedColIdx;
if isfield(handles.tableMetadataCurationField2Col,'manual_curator')
  % todo: this should not be hardcoded against manual_curator
  manualCuratorColIdx = handles.tableMetadataCurationField2Col.manual_curator;
  tfManualCuration = true;
else
  tfManualCuration = false;
end

indicator = false(Nrows,1);
for c = 1:Nrows
    row = data(c,:);

    tfExcluded = row{excludedColIdx};
    if tfExcluded
        indicator(c) = true;
        continue;
    end

    if tfManualCuration && ~isempty(row{manualCuratorColIdx})
        indicator(c) = true;
    end
end
zlclSortByData(handles.tblCuration,indicator,'descend');

%%% BUTTONS
function pbEnterCurationMarks_Callback(hObject,eventdata,handles)
zlclActionEnterCurationMarks(handles);

function pbClearMarks_Callback(hObject,eventdata,handles)
zlclActionClearCurationMarks(handles);

function pbUnexcludeAll_Callback(hObject, eventdata, handles)
handles.browser.unexcludeAllExps();

function pbMarkAllUncuratedAsPass_Callback(hObject, eventdata, handles)
handles.browser.curationMarkAllUncuratedAsPass();

%%% CONTEXT MENU
function ctxtMnuEnterCurationInfo(src,~)
zlclActionEnterCurationMarks(guidata(src));

function ctxtMnuClearCurationInfo(src,~)
zlclActionClearCurationMarks(guidata(src));

function ctxtMnuSortAscending(src,~)
zlclSortHelper(guidata(src),'ascend');

function ctxtMnuSortDescending(src,~)
zlclSortHelper(guidata(src),'descend');

function ctxtMnuBringExcludedMarkedToTop(src,~)
zlclBringExcludedMarkedToTop(guidata(src));

%%% REGULAR MENU
function mnuFileLoad_Callback(hObject, eventdata, handles)
handles.browser.curationLoadCurationFile();

function mnuFileAppend_Callback(hObject, eventdata, handles)
handles.browser.curationAppendCurationDataToFile();

function mnuFileWriteNew_Callback(hObject, eventdata, handles)
handles.browser.curationWriteCurationDataToNewFile();

function mnuSortAscending_Callback(hObject,eventdata,handles)
zlclSortHelper(handles,'ascend');

function mnuSortDescending_Callback(hObject,eventdata,handles)
zlclSortHelper(handles,'descend');

function mnuBringExcludedMarkedToTop_Callback(hObject,eventdata,handles)
zlclBringExcludedMarkedToTop(handles);

author: kristinbranson

repo: JAABA


text:
% function [x,y,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);
% Input:
%    I = image
%    par = vector for 4 parameters
%      [number of filter orientations, number of scales, filter size, elongation]
%      To use default values, put 0.
%    threshold = threshold on edge strength
% Output:
%    [x,y,gx,gy] = locations and gradients of an ordered list of edgels
%       x,y could be horizontal or vertical or 45 between pixel sites
%       but it is guaranteed that there [floor(y) + (floor(x)-1)*nr]
%       is ordered and unique.  In other words, each edgel has a unique pixel id.
%    par = actual par used
%    threshold = actual threshold used
%    mag = edge magnitude
%    mage = phase map
%    g = gradient map at each pixel
%    [FIe,FIo] = odd and even filter outputs
%    mago = odd filter output of optimum orientation
%
% Stella X. Yu, 2001

function [x,y,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);

if nargin<3 | isempty(threshold),
    threshold = 0.2;
end

[r,c] = size(I);
def_par = [8,1,20,3];

% take care of parameters, any missing value is substituted by a default value
if nargin<2 | isempty(par),
   par = def_par;
end
par(end+1:4)=0;
par = par(:);
j = (par>0);
have_value = [ j, 1-j ];
j = 1; n_filter = have_value(j,:) * [par(j); def_par(j)];
j = 2; n_scale  = have_value(j,:) * [par(j); def_par(j)];
j = 3; winsz    = have_value(j,:) * [par(j); def_par(j)];
j = 4; enlong   = have_value(j,:) * [par(j); def_par(j)];

% always make filter size an odd number so that the results will not be skewed
j = winsz/2;
if not(j > fix(j) + 0.1),
    winsz = winsz + 1;
end

% filter the image with quadrature filters
FBo = make_filterbank_odd2(n_filter,n_scale,winsz,enlong);
FBe = make_filterbank_even2(n_filter,n_scale,winsz,enlong);
n = ceil(winsz/2);
f = [fliplr(I(:,2:n+1)), I, fliplr(I(:,c-n:c-1))];
f = [flipud(f(2:n+1,:)); f; flipud(f(r-n:r-1,:))];
FIo = fft_filt_2(f,FBo,1);
FIo = FIo(n+[1:r],n+[1:c],:);
FIe = fft_filt_2(f,FBe,1);
FIe = FIe(n+[1:r],n+[1:c],:);

% compute the orientation energy and recover a smooth edge map
% pick up the maximum energy across scale and orientation
% even filter's output: as it is the second derivative, zero cross localize the edge
% odd filter's output: orientation
mag = sqrt(sum(FIo.^2,3)+sum(FIe.^2,3));
mag_a = sqrt(FIo.^2+FIe.^2);
[tmp,max_id] = max(mag_a,[],3);
base_size = r * c;
id = [1:base_size]';
mage = reshape(FIe(id+(max_id(:)-1)*base_size),[r,c]);
mage = (mage>0) - (mage<0);

ori_incr=pi/n_filter; % to convert jshi's coords to conventional image xy
ori_offset=ori_incr/2;
theta = ori_offset+([1:n_filter]-1)*ori_incr; % orientation detectors
% [gx,gy] are image gradient in image xy coords, winner take all
mago = reshape(FIo(id+(max_id(:)-1)*base_size),[r,c]);
ori = theta(max_id);
ori = ori .* (mago>0) + (ori + pi).*(mago<0);
gy = mag .* cos(ori);
gx = -mag .* sin(ori);
g = cat(3,gx,gy);

% phase map: edges are where the phase changes
mag_th = max(mag(:)) * threshold;
eg = (mag>mag_th);
h = eg & [(mage(2:r,:) ~= mage(1:r-1,:)); zeros(1,c)];
v = eg & [(mage(:,2:c) ~= mage(:,1:c-1)), zeros(r,1)];
[y,x] = find(h | v);
k = y + (x-1) * r;
h = h(k);
v = v(k);
y = y + h * 0.5; % i
x = x + v * 0.5; % j
t = h + v * r;
gx = g(k) + g(k+t);
k = k + (r * c);
gy = g(k) + g(k+t);

author: jwyang


text:
function aa = fitnesstra2(pop,para)
 global angles
[pops,numvar]=size(pop);
for i=1:pops
    time=pop(i,5)+pop(i,6);
    kk=trajt(para,pop(i,:));
    angles= kk(1:2,[5,7,10,13,15,17,20,23,25,27,30,33,35]);
    tt=torque(kk);
    ft=ftorque(tt);
    poo=kk(1:2,:);
    fq=sum(sum(abs(diff(poo'))));
    pos=forkin(poo);xx=pos(1,:);yy=pos(2,:);
    x=diff(xx);
    y=diff(yy);
    dis=sqrt(x.^2+y.^2);
    fdis=sum(dis);
    a1=2.3;
    a2=1.8;
    a3=2;
    a4=1;
    f=a1*ft+a2*fq+a3*fdis+a4*time;
    fob=fobstacle2(xx,yy);
    aa(i)=fob/f;
end

author: Sable

repo: mcbench-benchmarks


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

%%
clear; close all; clc;

%% Syntax
% |DET=ndet(N1,N2,N3);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function plotIPDF(f,varargin)
% plot orientations into inverse pole figures
%
% Syntax
%   plotIPDF(ori,[r1,r2,r3])
%   plotIPDF(ori,[r1,r2,r3],'points',100)
%   plotIPDF(ori,[r1,r2,r3],'points','all')
%   plotIPDF(ori,[r1,r2,r3],'contourf')
%   plotIPDF(ori,[r1,r2,r3],'antipodal')
%   plotIPDF(ori,data,[r1,r2,r3])
%
% Input
%  ebsd - @EBSD
%  r   - @vector3d specimen directions
%
% Options
%  RESOLUTION - resolution of the plots
%  property   - user defined colorcoding
%
% Flags
%  antipodal - include <VectorsAxes.html antipodal symmetry>
%  complete  - plot entire (hemi)--sphere
%
% See also
% S2Grid/plot savefigure Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

[mtexFig,isNew] = newMtexFigure('datacursormode',@tooltip,varargin{:});

% maybe we should call this function with the option add2all
if ~isNew && ~check_option(varargin,'parent') && ...
    ((((ishold(mtexFig.gca) && nargin > 1 && isa(varargin{1},'vector3d') && length(varargin{1})>1))) || check_option(varargin,'add2all'))
  plot(ori,varargin{:},'add2all');
  return
end

% find inverse pole figure direction
r = [];
try r = getappdata(mtexFig.currentAxes,'inversePoleFigureDirection'); end
if isempty(r), r = varargin{1}; end
argin_check(r,'vector3d');

for ir = 1:length(r)

  if ir>1, mtexFig.nextAxis; end

  % the crystal directions
  h = f.orientation \ r(ir);

  if ~check_option(varargin,{'complete','noSymmetry'})
    h = h.project2FundamentalRegion;
  end

  %  plot
  [~,cax] = h.line('fundamentalRegion','doNotDraw',varargin{:});

  if isNew, mtexTitle(cax(1),char(r(ir),'LaTeX')); end

  setappdata(cax,'inversePoleFigureDirection',r(ir));
  set(cax,'tag','ipdf');
  setappdata(cax,'CS',f.CS);
  setappdata(cax,'SS',f.SS);

end

if isNew || check_option(varargin,'figSize')
  mtexFig.drawNow('figSize',getMTEXpref('figSize'),varargin{:});
end

% --------------- Tooltip function ------------------
function txt = tooltip(empt,eventdata) %#ok<INUSL>

pos = get(eventdata,'Position');
xp = pos(1); yp = pos(2);

rho = atan2(yp,xp);
rqr = xp^2 + yp^2;
theta = acos(1-rqr/2);

m = Miller(vector3d.byPolar(theta,rho),getappdata(gcf,'CS'));
m = round(m);
txt = char(m,'tolerance',3*degree,'commasep');

author: mtex-toolbox

repo: mtex

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


text:
function Results=weightedfit(data)
% This code fits makes a linear fit to a data set (using y =bx+a) where each data point
% has a different or constant standard deviation. Your data should have three or two columns.
% The first column should be the independent variable(x) and the second
% column should be the dependent variable(y). Column three should contain
% your standard deviations for each datapoint. In the situations where you
% do not specify a column three, the code assigns a weight of one to all
% data points and this corresponds to the regular linear fits.
%==========
% INPUTS
%==========
%data = 3 columns; column 1 = x, column2 = y and column 3 = standard dev.

%==========
%OUTPUTS
%==========
%Result.slope= b; Fitted slope
%Result.Intercept = a; Fitted intercept

%Coded by Ebo Ewusi-Annan (University of Florida, 2011)
%============
%REFERENCES
%===========
%1. Willam H. Press, Saul A. Teukolsky and Willan T. Vetterling (1997).
%Numerical Recipes in Fortran.
%2. Philip R. Bevington and D. Keith Robinson (2003). Data Reduction and
%Error Analysis for the Physical Sciences.

x= data(:,1);
y=data(:,2);
[s t]= size(data);
stdv=ones(s,1);
if t==3
stdv=data(:,3);
end
w = 1./stdv.^2;
S = sum(w);
Sx = sum(w.*x);
Sy = sum(w.*y);
Sxx= sum(w.*x.^2);
Sxy= sum(w.*x.*y);
Delta = S*Sxx - (Sx)^2;
a = (Sxx*Sy - Sx*Sxy)./Delta;
b = (S*Sxy - Sx*Sy)./Delta;
fprintf('\n slope=%f Int=%f \n',b,a)
Results.slope=b;
Results.Intercept= a;
y_fit = a + b*x;
  if t==2
    h=plot(x,y,'rs','MarkerFaceColor','r');
     title('Unweighted fit')
  else
     clf;set(gcf,'color','w');
     h=errorbar(x,y,stdv,'rs','MarkerFaceColor','r');
     title('Weighted fit')
  end
hold on; q=plot(x,y_fit,'b.--','linewidth',2);
xlabel('x (Column 1)')
ylabel('y (Column 2)')
legend([h(1),q(1)],'Data',sprintf('\nSlope=%f\nIntercept=%f\n',b,a),'location','Southeast')

end


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

%%
clear; close all; clc;

%% Syntax
% |[Fc_Q,Fc_T,Ft_Q,Ft_T,Ct_Q,Ct_T,Vt]=rhombicDodecahedronMesh(r,nCopies)|

%% Description
% Creates a rhombic dodecahedron mesh where r sets the radias and nCopies
% (a 1x3 vector) sets the number of copies in the x, y, and z direction.
% The output consists of:
%
% Fc_Q, Fc_T: the quadrilateral and triangular face cell arrays (1 cell
% entry per element).
%
% Ft_Q, Ft,T: the quadrilateral and triangular face arrays
%
% Ct_Q, Ct,T: color/label data for the face arrays
%
% Vt: the vertex array

%%
% Plot settings
fontSize=15;
faceAlpha1=0.8;

%% Creating a mesh of rhombic dodecahedra

r=0.5; %Radii, results in a width of 1
nCopies=[3 3 3]; %Number of offset copies

[E,V,C,F,CF]=rhombicDodecahedronMesh(r,nCopies);

%%
% Plotting results
cFigure; hold on;
title('A mesh of rhombicDodecahedra','FontSize',fontSize);
gpatch(F,V,CF,'k',faceAlpha1);
colormap(gjet);
axisGeom(gca,fontSize);
camlight('headlight');
drawnow;

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

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


text:
%PDF2EPS  Convert a pdf file to eps format using pdftops
%
% Examples:
%   pdf2eps source dest
%
% This function converts a pdf file to eps format.
%
% This function requires that you have pdftops, from the Xpdf suite of
% functions, installed on your system. This can be downloaded from:
% http://www.foolabs.com/xpdf
%
%IN:
%   source - filename of the source pdf file to convert. The filename is
%            assumed to already have the extension ".pdf".
%   dest - filename of the destination eps file. The filename is assumed to
%          already have the extension ".eps".

% Copyright (C) Oliver Woodford 2009-2010

% Thanks to Aldebaro Klautau for reporting a bug when saving to
% non-existant directories.

function pdf2eps(source, dest)
% Construct the options string for pdftops
options = ['-q -paper match -eps -level2 "' source '" "' dest '"'];
% Convert to eps using pdftops
[status, message] = pdftops(options);
% Check for error
if status
    % Report error
    if isempty(message)
        error('Unable to generate eps. Check destination directory is writable.');
    else
        error(message);
    end
end
% Fix the DSC error created by pdftops
fid = fopen(dest, 'r+');
if fid == -1
    % Cannot open the file
    return
end
fgetl(fid); % Get the first line
str = fgetl(fid); % Get the second line
if strcmp(str(1:min(13, end)), '% Produced by')
    fseek(fid, -numel(str)-1, 'cof');
    fwrite(fid, '%'); % Turn ' ' into '%'
end
fclose(fid);
end

author: goGPS-Project


text:
function p_ = ViewRanking(X,p,Lower,Upper)

[J,N]=size(X);
K=length(Lower);

% constrain probabilities to sum to one...
Aeq = ones(1,J);
beq=1;

% ...constrain the expectations...
V=X(:,Lower) - X(:,Upper);

A = V';
b = 0;

% ...compute posterior probabilities
p_ = EntropyProg(p,A,b,Aeq ,beq);


text:
function design = ula_1d(n, d, name)
%ULA_1D Generates a 1D ULA.
%Syntax:
%   design = ULA_1D(10, wavelength/2);
%   design = ULA_1D(10, 0.5, 'ULA with 10 sensors');
%Inputs:
%   n - Number of elements.
%   d - Inter-element spacing.
%   name - Custom name of the array. Default is 'ULA with n elements'.
%Outputs:
%   design - An array design struct.
if n <= 0
    error('n_sensor must be a positive integer.');
end
if d <= 0 || ~isreal(d)
    error('d must be a positive real number.');
end
if nargin <= 2
    name = sprintf('ULA with %d elements', n);
elseif ~ischar(name)
    error('Name must be a string.');
end
design.element_indices = (0:n-1);
design.element_positions = design.element_indices*d;
design.element_spacing = d;
design.element_count = n;
design.dim = 1;
design.type = 'ula';
design.name = name;
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function out =  bicubic_interpolation_cell(p, x, y)
v(1) = cubic_interpolation_cell(p(1,:), y);
v(2) = cubic_interpolation_cell(p(2,:), y);
    v(3) = cubic_interpolation_cell(p(3,:), y);
    v(4) = cubic_interpolation_cell(p(4,:), y);
   out =  cubic_interpolation_cell(v, x);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Yagiz Aksoy
% 2017

tbloc = fileparts(mfilename('fullpath'));
addpath(fullfile(tbloc, 'affinity'));
addpath(fullfile(tbloc, 'trimming'));
addpath(fullfile(tbloc, 'common'));
clear tbloc;

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
function pass = test_mean( )
% Test MEAN()

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

g = diskfun();
pass(1) = isempty(mean(g));

%%
g = diskfun(@(x,y) 0*x+1);
pass(2) = ( norm(mean(g)-.5) < tol ); %this direction includes measure on disk
pass(3) = ( norm(mean(g,2)-1)<tol );

g = diskfun(@(t,r) r.^2.*sin(2*t), 'polar');
pass(4) = ( norm(mean(g)-chebfun(@(t) .25*sin(2*t), [-pi pi], 'trig') ) < tol );
pass(5) = ( norm(mean(g,2)) < tol);

g = diskfun( @(x,y) exp(-10*x.^2-10*y.^2));
pass(6) = ( abs(sum(sum(g)) - sum2(g)) < tol );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function linpack_s_test19 ( )

%*****************************************************************************80
%
%% TEST19 tests SPPCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST19\n' );
  fprintf ( 1, '  For a positive definite symmetric packed matrix,\n' );
  fprintf ( 1, '  SPPCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the matrix A.
%
  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      if ( i == j - 1 )
        a(k) = -1.0;
      elseif ( i == j )
        a(k) = 2.0;
      else
        a(k) = 0.0;
      end
    end
  end
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition number.\n' );

  [ a, rcond, z, info ] = sppco ( a, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reciprocal condition number = %f\n', rcond );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Nystrom_method
clc;
clear;
close all;
%% the example fuction
datasize=200;
var_noise=0.1;
x=linspace(-10,10,datasize)';
    f=sin(x)./x;
noise=var_noise*randn(datasize,1);
y=f+noise;
x_raw=x;
plot(x_raw,f,'r-','LineWidth',2),hold on;
%% the training data samples
Inset=30;
dataindex=randsample(datasize,Inset);
xI=x(dataindex);
yI=y(dataindex);
plot(xI,yI,'b+'),hold on;
x(dataindex)=[];
y(dataindex)=[];
x=[xI;x];
y=[yI;y];
%% the covariance matrix,kernel;
K_mat=Kernel_func(x,x,0.7,1.2);
K_mm=K_mat(1:Inset,1:Inset);

%%
[vecmat_m,valmat_m]=eig(K_mm);
valmat_m=diag(valmat_m);
[valmat_m,ix]=sort(valmat_m,'descend');
vecmat_m=vecmat_m(:,ix);
valmat_n=datasize/Inset*valmat_m;
prin_val=sum(valmat_n>0);

for i=1:prin_val;
vecmat_n(:,i)=sqrt(Inset/datasize)/valmat_m(i)*K_mat(1:datasize,1:Inset)*...
    vecmat_m(:,i);
valmat_n(i)=datasize/Inset*valmat_m(i);
end

valmat_n=diag(valmat_n(1:prin_val));

new_kernel=inv(valmat_n)+vecmat_n'*1/var_noise*eye(datasize)*vecmat_n;
% rank(new_kernel)
L=chol(new_kernel,'lower');  %using cholesky decomposition to be more stable
%%test
testnum=500;
x_test=linspace(-15,15,testnum)';
K_test_train=Kernel_func(x_test,x,0.7,1.2);
K_test=Kernel_func(x_test,x_test,0.7,1.2);

b_temp=vecmat_n'*1/var_noise*eye(datasize)*y;
f_mean=K_test_train*(1/var_noise*eye(datasize)*y-1/var_noise*eye(datasize)*...
    vecmat_n*(L'\(L\b_temp)));
for i=1:testnum
K_ss=Kernel_func(x_test(i),x_test(i),0.7,1.2);
K_s=Kernel_func(x_test(i),x,0.7,1.2);
v=L\( vecmat_n'*1/var_noise*eye(datasize)*K_s');
  f_var(i)=K_ss-(K_s*1/var_noise*eye(datasize)*K_s'-v'*v);
%  f_var(i)=K_ss-K_s*(1/var_noise*eye(datasize)-1/var_noise*eye(datasize)*vecmat_n...
%      *inv(new_kernel)*vecmat_n'*1/var_noise*eye(datasize))*K_s';
end
f_var=f_var';

% in_pos=find(f_var>0);
% f_var=f_var(in_pos);
% x_test=x_test(in_pos);
% f_mean=f_mean(in_pos);

plot(x_test,f_mean,'b-','LineWidth',2),hold on;
%%plot the variance
x_area=[x_test;flipdim(x_test,1)];
ybound_std=[(f_mean+1.96*sqrt(f_var));flipdim(f_mean-1.96*sqrt(f_var),1)];
fill(x_area,ybound_std,[187 255 0]/255,'EdgeColor','none');
plot(x_raw,f,'r-','LineWidth',2),hold on;
plot(xI,yI,'b+'),hold on;
plot(x_test,f_mean,'b-','LineWidth',2),hold on;
%  alpha(0.6)
legend1=legend('target function','training points','preditive distribution','95% confidence interval');
 set(legend1,'Box','off','Color','none',...
     'Location','NorthEast');
xlabel('30 active training points')
set(gca,'XTick',[-15:3:15])
matlab2tikz( 'Nys_data30.tex' )
end
%%

author: epfl-lasa


text:
% TEST_WAVELET_2D Test case for WAVELET_2D
%
% See also
%   WAVELET_2D
classdef test_wavelet_2d < matlab.unittest.TestCase
    methods(Test)

        function testWithNoOptions(testcase)
            %% define
            x = rand(64, 64);
            filters = morlet_filter_bank_2d(size(x));
            %% with no options
            [x_phi, x_psi] = ...
                wavelet_2d(x, filters);
            %% check number of high pass filter
            J = filters.meta.J;
            L = filters.meta.L;
            Q = filters.meta.Q;
            expected = J*L*Q;
            actual = numel(x_psi.signal);
            %% assert
            testcase.assertEqual(expected, actual);

        end

         function testWithRandomOptions(testcase)
            for i = 1:32
                %%
                sz = 1 + floor(128*rand(1,2));
                %% define
                x = rand(sz);

                white_list={'x_resolution','oversampling'};
                type={'d','d'};
                values={{1,2,4,8,16},{0,1}};

                options_W = generate_random_options(white_list,type,values);

                filters = morlet_filter_bank_2d(sz);
                %% with no options
                [x_phi, x_psi] = ...
                    wavelet_2d(x, filters,options_W);
                %% check number of high pass filter
                J = filters.meta.J;
                L = filters.meta.L;
                Q = filters.meta.Q;
                expected = J*L*Q;
                actual = numel(x_psi.signal);
                %% assert
                testcase.assertEqual(expected, actual);

            end
         end

        function testWithRandomSize(testcase)
            for i = 1:32
                %%
                sz = 1 + floor(128*rand(1,2));
                %% define
                x = rand(sz);
                filters = morlet_filter_bank_2d(sz);
                %% with no options
                [x_phi, x_psi] = ...
                    wavelet_2d(x, filters);
                %% check number of high pass filter
                J = filters.meta.J;
                L = filters.meta.L;
                Q = filters.meta.Q;
                expected = J*L*Q;
                actual = numel(x_psi.signal);
                %% assert
                testcase.assertEqual(expected, actual);

            end

        end
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function h=sgsf_2d(x,y,px,py,flag_coupling)
%Function:
%       2-D Savitzky-Golay smoothing filter (i.e., the polynomial smoothing
%       filter, or the least-squares smoothing filter)
%       See Ref. [1] for details on the 1-D Savitzky-Golay smoothing filter.
%       One can also refer to the following URL where a program of
%       1-D Savitzky-Golay smoothing (and differentiation) filter is given:
%       http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4038&objectType=file
%       See Ref. [2] and [3] for details on the 2-D Savitzky-Golay smoothing filter.
%Usage:
%       h=sgsf_2d(x,y,px,py,flag_coupling)
%       x    = x data point, e.g., -3:3
%       y    = y data point, e.g., -2:2
%       px    =x polynomial order       default=1
%       py    =y polynomial order       default=1
%       flag_coupling  = with or without the consideration of the coupling terms, between x and y. default=0
%Example:
%       sgsf_2d(-3:3,-3:3,2,2)
%       sgsf_2d(-3:3,-3:3,2,2,1)
%Author:
%       Jianwen Luo <luojw@bme.tsinghua.edu.cn, luojw@ieee.org> 2003-12-15
%       Department of Biomedical Engineering, Department of Electrical Engineering
%       Tsinghua University, Beijing 100084, P. R. China
%Reference
%[1]A. Savitzky and M. J. E. Golay, "Smoothing and Differentiation of Data by Simplified Least Squares Procedures,"
%   Analytical Chemistry, vol. 36, pp. 1627-1639, 1964.
%[2]K. L. Ratzlaff and J. T. Johnson, "Computation of Two-Dimensional Polynomial Least-Squares Convolution Smoothing Integers,"
%   Analytical Chemistry, vol. 61, pp. 1303-1305, 1989.
%[3]J. E. Kuo, H. Wang, and S. Pickup, "Multidimensional Least-Squares Smoothing Using Orthogonal Polynomials,"
%   Analytical Chemistry, vol. 63, pp. 630-635, 1991.

if nargin<5
    flag_coupling=0;
end

if nargin<4
    py=1;
end

if nargin<3
    px=1;
end

[x,y]=meshgrid(x,y);
[ly,lx]=size(x);

x=x(:);
y=y(:);

A=[];

if flag_coupling
    for kx=px:-1:0
        for ky=py:-1:0
            A=[A x.^kx.*y.^ky];
        end
    end
else
    for k=px:-1:1
        A=[A x.^k];
    end
    for k=py:-1:0
        A=[A y.^k];
    end
end

h=inv(A'*A)*A';
h=h(size(h,1),:);
h=reshape(h,ly,lx);


text:
function timespderiv(options, flags)
% TIMESPDERIV   Measure the performance of SPINTERP with derivatives
%    TIMESPDERIV(OPTIONS)  Script file to plot a graph showing the
%    ratio of computing BOTH interpolated values AND gradients to
%    computing interpolated values ONLY (at 100 points simultaneously).
%    By default, the Clenshaw-Curtis-grid with sparse indices data is
%    used. Alternatively, the Chebyshev grid and other options
%    may be selected using an OPTIONS structure created with SPSET.
%
%    TIMESPDERIV(OPTIONS, FLAGS)  FLAGS is an optional struct that
%    customizes the plots.
%    The following FLAGS are available:
%       continuousDerivatives [ 'on' | 'off' ]  turns on the augmented
%         continuous derivatives available for the Clenshaw-Curtis
%         sparse grid.
%       absoluteTime [ 'on' | 'off' ]  plots the absolute measured
%         time in seconds instead of the ratio.
%       nPoints [ integer {100} ]  can be used to change the default
%         number of interpolated points.
%
%  Examples:
%    Plot ratios for CC-Grid, discontinuous derivatives:
%       timespderiv;
%
%    Plot ratios for CC-Grid, augmented continuous derivatives:
%       flags = struct('continuousDerivatives', 'on');
%       timespderiv([], flags);
%
%    Plot ratios for CGL-Grid:
%       timespderiv(spset('GridType','Chebyshev'));
%
%    Plot absolute time measurements for CGL-Grid:
%       flags = struct('absoluteTime', 'on');
%       timespderiv(spset('GridType','Chebyshev'), flags);
%
%    See also SPSET, SPINTERP, TESTFUNCTIONS.

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.0
% Date   : July 8, 2006

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

if nargin < 1, options = []; end
if nargin < 2, flags = []; end

continuousDerivatives = 0;
if isfield(flags, 'continuousDerivatives')
	if strcmpi(flags.continuousDerivatives, 'on')
		continuousDerivatives = 1;
	end
end

absoluteTime = 0;
if isfield(flags, 'absoluteTime')
	if strcmpi(flags.absoluteTime, 'on')
		absoluteTime = 1;
	end
end

% Number of points to interpolate
nsamplepoints = 100;
if isfield(flags, 'nPoints')
	nsamplepoints = flags.nPoints;
end

sparseIndices = spget(options, 'SparseIndices', 'auto');
gridtype = spget(options, 'GridType', 'Clenshaw-Curtis');

% Set the problem dimensions and the maximum discretization levels
% for each dimension.
switch lower(gridtype)
 case 'clenshaw-curtis'
	d = [1,2,4,8,16];
  n = [17,13,10,6,4];
 case 'noboundary'
	error('MATLAB:timespderiv:badopt',['Computing derivatives not ' ...
	      'supported for grid type ''' gridtype '''.']);
 case 'maximum'
	error('MATLAB:timespderiv:badopt',['Computing derivatives not ' ...
	      'supported for grid type ''' gridtype '''.']);
 case 'chebyshev'
	d = [1,2,4,8,16];
  n = [17,13,10,6,4];
 otherwise
	error('MATLAB:timespderiv:badopt','Unknown grid type.');
end

ndims = length(n);
markers = ['s', '.', 'd', '+', 'o', 'x', '*', 'v', '^'];

%  Compute random constants w and c for Gerz' test functions.
w = rand(d(end),1);
c = rand(d(end),1);
sumc = sum(c,1);
sumw = sum(w,1);
c = 1.5.*c/sumc;
w = w/sumw;

time = zeros(ndims,n(1)+1);
npoints = zeros(ndims,n(1)+1);

for k = 1:d(end)
	v{k} = rand(nsamplepoints,1);
end

for m = 1:ndims
	disp(['Current dim: ' num2str(d(m))]);
	z = [];
	for l = 0:n(m)
		disp(['Current level n = ' num2str(l) '...']);
		options = spset('MinDepth',l,'MaxDepth',l,'Vectorized','on',...
										'GridType', gridtype, 'SparseIndices', ...
										sparseIndices, 'PrevResults', z);
		z = spvals('testfunctions',d(m),[],options,1,c(1:d(m)),w(1: ...
																									d(m)));
		t = z.surplusCompTime;
		disp(['Computing sparse grid points and evaluating function took ' ...
					num2str(z.fevalTime) ' [s].']);
		disp(['Computing hierarchical surpluses took ' num2str(t) ' [s].']);
		npoints(m,l+1) = z.nPoints;
		if continuousDerivatives
			z.continuousDerivatives = 'on';
		end
		t = 0; k = 0;
		tic;
		while t < 0.5
			y = spinterp(z,v{1:d(m)});
			t = toc;
			k = k + 1;
		end
		t = t / k;
		t2 = 0; k = 0;
		tic;
		while t2 < 0.5
			[y,dy] = spinterp(z,v{1:d(m)});
			t2 = toc;
			k = k + 1;
		end
		t2 = t2 / k;
		if absoluteTime
			time(m,l+1) = t2;
		else
			time(m,l+1) = t2/t;
		end
		disp([num2str(nsamplepoints) ...
					' evals took ' num2str(t) ' [s], with gradients ' ...
					num2str(t2) ' [s]. Ratio: ' num2str(t2/t)]);
		disp(' ');
	end
	disp('Done!');
	disp(' ');
end

% Plot results
clf;
h = loglog(npoints',time','LineWidth',1);
for k = 1:ndims
	set(h(k), 'Marker', markers(k));
	hold on;
end
axis tight;
hold on;

% Add legend and axis labels
hold off;
if absoluteTime
	if ~continuousDerivatives
		title(['Time to compute ' num2str(nsamplepoints) ...
					' values with gradients [' gridtype ' grid].']);
	else
		title(['Time to compute ' num2str(nsamplepoints) ...
					' vals with CONTINUOUS grads [' gridtype ' grid].']);
	end
	ylabel('Time [s]')
else
	if ~continuousDerivatives
		title(['Ratio of interpolating ' num2str(nsamplepoints) ...
					' values with/without gradients [' gridtype ' grid].']);
	else
		title(['Ratio of interpol. ' num2str(nsamplepoints) ...
					' vals with/without CONTINUOUS grads [' gridtype ' grid].']);
	end
	ylabel('Ratio [-]')
	set(gca,'YTick',[1,2,5,10,20,50]);
	set(gca,'YLim',[0.8,50]);
end
xlabel('number of nodes N')
s = {};
for k = 1:ndims
	s{k} = ['d = ' num2str(d(k))];
end

if strcmpi(gridtype, 'chebyshev')
	legend(s,2);
else
	legend(s,1);
end


text:
function [x] = epp1(v, rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Element-wise Soft Thesholding
%
% x= sign(v) max( |v|- rho, 0)
%
% Authors: Yin Li @ PAMI Lab SJTU
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: andrewssobral

repo: lrslibrary


text:
% Plot nodes radially out from a given center. Equidistant nodes
% have the same radius, but different angles. Works best as a quick
% visualization for trees, or very sparse graphs.
%
% Note 1: If a center node is not specified, the nodes are ordered by
%   sum of neighbor degrees, and the node with highest sum is plotted
%   in the center.
% Note 2: The graph has to be connected.
%
% Inputs: adjacency matrix (nxn), and center node (optional)
% Outputs: plot
%
% Other routines used: sortNodesBySumNeighborDegrees.m,
%                      adj2adjL.m, diameter.m, kminNeighbors.m
% GB: last updated, December 6 2012

function [] = radialPlot(adj,i0)

d=diameter(adj);
if d==Inf; printf('The graph is not connected.\n'); return; end

if nargin==1; I=sortNodesBySumNeighborDegrees(adj); i0=I(1); end;

L=adj2adjL(adj);

circles={};   % per circle
circle{i0}=0; % per node

for k=1:d
  kneigh = kminNeighbors(adj,i0,k);
  if not(isempty(kneigh));
    circles{k}=kneigh;
    for nei=1:length(kneigh); circle{kneigh(nei)}=k; end
  end;
end

for i=1:size(adj,1); descendants{i} = 0; end;

for k=length(circles)-1:-1:1

  for i=1:length(circles{k})
    node = circles{k}(i);
    neighs=L{node};

    for nei=1:length(L{node})
      nei = L{node}(nei);
      if sum(find(circles{k+1}==nei))>0; descendants{node}=descendants{node}+descendants{nei}+1; end
    end
  end

end

descendants{i0}=descendants{i0} + length(L{i0});
for i=1:length(L{i0}); descendants{i0}=descendants{i0} + descendants{L{i0}(i)}; end;

start_angle{i0}=0;
range{i0}=2*pi;

queue=[i0]; visited=[];

while not(isempty(queue))

  node=queue(1); queue=queue(2:length(queue)); visited=[visited node];

  % select neighbours that have not been visited
  neigh=[];
  for ll=1:length(L{node})
    if sum(find(visited==L{node}(ll)))==0 && sum(find(queue==L{node}(ll)))==0;
      neigh=[neigh L{node}(ll)]; end;
  end

  if isempty(neigh); continue; end;

  range{neigh(1)} = range{node}*(descendants{neigh(1)}+1)/descendants{node};
  start_angle{neigh(1)}=start_angle{node};
  if sum(find(queue==neigh(1)))==0; queue=[queue neigh(1)]; end;

  for j=2:length(neigh)
    range{neigh(j)} = range{node}*(descendants{neigh(j)}+1)/descendants{node};
    start_angle{neigh(j)} = start_angle{neigh(j-1)} + range{neigh(j-1)};
    if sum(find(queue==neigh(j)))==0; queue=[queue neigh(j)]; end;
  end

end

set(gcf,'Color',[1 1 1],'Colormap',hot);   % Modify color scheme here
map=colormap('hot');

% plot with the angle and radius (circle) information
for i=1:size(adj,1)
  x = circle{i}*cos(start_angle{i}+range{i}/2);
  y = circle{i}*sin(start_angle{i}+range{i}/2);
  plot(x,y,'.','Color',map(mod(5*circle{i},length(map))+1,:)); hold off; hold on;
end

xs = cos([0:1:30]*2*pi/30);
ys = sin([0:1:30]*2*pi/30);
for c=1:length(circles)
  plot(c*xs,c*ys,':','Color',[244/255,244/255,244/255]);
  hold off; hold on;
end

% add the links between nodes
for i=1:size(adj,1)
  for j=i+1:size(adj,1)
    if adj(i,j)>0
      xi=circle{i}*cos(start_angle{i}+range{i}/2);
      yi=circle{i}*sin(start_angle{i}+range{i}/2);

      xj=circle{j}*cos(start_angle{j}+range{j}/2);
      yj=circle{j}*sin(start_angle{j}+range{j}/2);

      line([xi xj],[yi yj],'Color',map(mod(5*circle{i},length(map))+1,:))
    end
  end
end

axis tight
axis off

author: aeolianine

repo: octave-networks-toolbox


text:
function [u, uu, erriter, num, timet] = test_CMF3D_ML
%%
%   Function test_CMF3D_ML
%
%   The matlab function to show how to use the functions:
%               CMF3D_ML_mex and CMF3D_ML_GPU
%
%   Before using the function CMF3D_ML_mex, you should compile its source file as follows:
%
%       >> mex CMF3D_ML_mex.c
%
%   Before using the function CMF3D_ML_GPU, you should compile its
%   GPU source file (for Windows):
%
%       >> nvmex -f nvmexopts.bat CMF3D_ML_GPU.cu -IC:\cuda\v4.0\include ...
%       >> -LC:\cuda\v4.0\lib\x64 -lcufft -lcudart
%
%   GPU compilation depends on your system and its configurations!
%
%   After compilation, you can define all the parameters (penalty, C_t, para) as follows:
%
%        - penalty: point to the edge-weight penalty parameters to
%                   total-variation function.
%
%          For the case without incorporating image-edge weights,
%          penalty is given by the constant everywhere. For the case
%          with image-edge weights, penalty is given by the pixelwise
%          weight function:
%
%          for example, penalty(x) = b/(1 + a*| grad f(x)|) where b,a > 0 .
%
%        - C_t(x,i=1...nlab): point to the capacities of
%                 sink flow fields pt(x,i=1...nlab);
%
%        - para: a sequence of parameters for the algorithm
%             para[0,1,2]: rows, cols, heights of the given 3D volume
%             para[3]: the number of labels or regions
%             para[4]: the maximum iteration number
%             para[5]: the error bound for convergence
%             para[6]: cc for the step-size of augmented Lagrangian method
%             para[7]: the step-size for the graident-projection step to the
%                    total-variation function. Its optimal range is [0.06, 0.12].
%
%   Outputs:
%
%       - u: the computed continuous labeling result u(x,i = 1...nlab) in [0,1],
%           where nlab is the number of labels/volumes.
%
%       - uu: the final labeled result uu(x), computed by the u(x,i).
%
%       - erriter: it returns the error evaluation of each iteration,
%           i.e. it shows the convergence rate. One can check the algorithm
%           performance.
%
%       - i: gives the total number of iterations, when the algorithm converges.
%
%       - timet: gives the total computation time.
%
%       Example:
%
%             >> [u, erriter, i, timet] = CMF3D_ML_GPU(single(penalty), single(Cs), single(Ct), single(para));
%
%             >> You can choose a 2D slice to demonstrate the final result.
%
%             >> figure, loglog(erriter,'DisplayName','erriterN');figure(gcf)
%
%
%   The original algorithm was proposed in the following papers:
%
%   [1] Yuan, J.; Bae, E.;  Tai, X.-C.
%       A Study on Continuous Max-Flow and Min-Cut Approaches
%       CVPR, 2010
%
%   [2] Yuan, J.; Bae, E.; Tai, X.-C.; Boycov, Y.
%       A Continuous Max-Flow Approach to Potts Model
%       ECCV, 2010
%
%   The mimetic finite-difference discretization method was proposed for
%   the total-variation function in the paper:
%
%   [1] Yuan, J.; Schn{\"o}rr, C.; Steidl, G.
%       Simultaneous Optical Flow Estimation and Decomposition
%       SIAM J.~Scientific Computing, 2007, vol. 29, page 2283-2304, number 6
%
%   This software can be used only for research purposes, you should cite ALL of
%   the aforementioned papers in any resulting publication.
%
%
%   Please email Jing Yuan (cn.yuanjing@gmail.com) for any questions,
%   suggestions and bug reports
%
%   The Software is provided "as is", without warranty of any kind.
%
%               Version 1.0
%   https://sites.google.com/site/wwwjingyuan/
%
%   Copyright 2011 Jing Yuan (cn.yuanjing@gmail.com)
%

nii = load_nii('T1w_fa18_restore_brain.nii');
ur = double(nii.img);
umax = max(max(max(ur)));
umin = min(min(min(ur)));
ur = (ur - umin)/(umax-umin);
ur = ur(21:120, 70:170, 69:170);

% define the label information

nlab = 4;
ulab(1) = 0;
ulab(2) = 0.16;
ulab(3) = 0.4;
ulab(4) = 0.53;

[rows, cols, heights] = size(ur);

% build up the data terms

Ct = zeros(rows,cols,heights,nlab);

for i=1:nlab
    Ct(:,:,:,i) = abs(ur - ulab(i))*10;
end

% set up the input parameter vector
cc = 0.25;
varParas = [rows; cols; heights; nlab; 200; 1e-4; cc; 0.11];
%                para 0,1,2 - rows, cols, heights
%                para 3 - the total number of labels
%                para 4 - the maximum number of iterations
%                para 5 - the error bound for convergence
%                para 6 - cc for the step-size of augmented Lagrangian method
%                para 7 - the step-size for the graident-projection of the spatial flow fields p(x,i)

penalty = ones(rows,cols,heights);

% ----------------------------------------------------------------------
%  Use the function CMF3D_ML_mex to run the algorithm on CPU
% ----------------------------------------------------------------------

% [u, erriter,num,timet] = CMF3D_ML_mex(single(penalty), single(Ct), single(varParas));

% ----------------------------------------------------------------------
%  Use the function CMF3D_ML_GPU to run the algorithm on GPU
% ----------------------------------------------------------------------

[u, erriter,num,timet] = CMF3D_ML_GPU(single(penalty), single(Ct), single(varParas));

% The output u must be multiplied by cc to scale it in [0,1],
% due to the implementation
u = cc*u;

[um,I] = max(u, [], 4);

uu = zeros(rows,cols,heights);

for k=1:rows
    for j=1:cols
        for i=1:heights
            uu(k,j,i) = ulab(I(k,j,i));
        end
    end
end


text:
function [ Tree,Cost ] =  UndirectedMaximumSpanningTree (CostMatrix)
% The function takes CostMatrix as input and returns the maximum spanning tree T
% Uses Kruskal's Algorithm
% Extract the edge weights from the cost matrix
% Sort the edges in a non decreasing order of weights
% This algorithm is revised by Lowell Guangdi at 2009/06/11.

n = size (CostMatrix,1); %Number of vertices
EdgeWeights = 0;         %Edges and corresponding weights
EdgeWeightsCounter = 0;
for i = 1:n
    for j = (i+1):n
        if ((CostMatrix(i,j))~=0)
            EdgeWeightsCounter = EdgeWeightsCounter + 1;
            EdgeWeights(EdgeWeightsCounter,1) = CostMatrix(i,j);
            EdgeWeights(EdgeWeightsCounter,2) = i;
            EdgeWeights(EdgeWeightsCounter,3) = j;
        end
    end
end

SortedEdgeWeights = 0;
SortedEdgeWeights = sortrows(EdgeWeights);
% First column of SortedEdgeWeights are the weights
% Second and third column are the vertices that the edges connect
m = size(SortedEdgeWeights,1); % number of edges

% We use the Disjoint sets data structures to detect cycle while adding new
% edges. Union by Rank with path compression is implemented here.

% Assign parent pointers to each vertex. Initially each vertex points to
% itself. Now we have a conceptual forest of n trees representing n disjoint
% sets
global ParentPointer ;
ParentPointer = 0;
ParentPointer(1:n) = 1:n;

% Assign a rank to each vertex (root of each tree). Initially all vertices
% have the rank zero.
TreeRank = 0;
TreeRank(1:n) = 0;

% Visit each edge in the sorted edges array
% If the two end vertices of the edge are in different sets (no cycle), add
% the edge to the set of edges in maximum spanning tree
MSTreeEdges = 0;
MSTreeEdgesCounter = 0; i = m;
while ((MSTreeEdgesCounter < (n-1)) && (i>=1))
%Find the roots of the trees that the selected edge's two vertices
%belong to. Also perform path compression.
    root1=0; root2=0; temproot=0;
    temproot = SortedEdgeWeights(i,2);
    root1 = FIND_PathCompression(temproot);

    temproot = SortedEdgeWeights(i,3);
    root2 = FIND_PathCompression(temproot);

    if (root1 ~= root2)
        MSTreeEdgesCounter = MSTreeEdgesCounter + 1;
        MSTreeEdges(MSTreeEdgesCounter,1:3) = SortedEdgeWeights(i,:);
        if (TreeRank(root1)>TreeRank(root2))
            ParentPointer(root2)=root1;
        else
            if (TreeRank(root1)==TreeRank(root2))
               TreeRank(root2)=TreeRank(root2) + 1;
            end
            ParentPointer(root1)=root2;
        end
    end
    i = i - 1;
end

MSTreeEdgesCounter = 0;
Tree = 0;
Tree(1:n,1:n)=0;
while (MSTreeEdgesCounter < (n-1))
    MSTreeEdgesCounter = MSTreeEdgesCounter + 1;
    Tree(MSTreeEdges(MSTreeEdgesCounter,2),MSTreeEdges(MSTreeEdgesCounter,3))=1;
    Tree(MSTreeEdges(MSTreeEdgesCounter,3),MSTreeEdges(MSTreeEdgesCounter,2))=1;
end
%T

Cost = 0;
for p = 1:n
    for q = p+1:n
       if Tree( p,q ) == 1
          Cost = Cost + CostMatrix( p,q );
       end
    end
end

end

function [parent] = FIND_PathCompression(temproot)

global ParentPointer;
ParentPointer(temproot);
if (ParentPointer(temproot)~=temproot)
    ParentPointer(temproot) = FIND_PathCompression(ParentPointer(temproot));
end
parent = ParentPointer(temproot);
end


text:
function ft_flush_event(filename, varargin)

% FT_FLUSH_EVENT removes all events from the event queue
%
% Use as
%   ft_flush_event(filename, ...)
%
% See also FT_FLUSH_HEADER, FT_FLUSH_DATA

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

% TODO implement filtering

% set the defaults
eventformat = keyval('eventformat', varargin); if isempty(eventformat), eventformat = ft_filetype(filename); end

switch eventformat
  case 'disp'
    % nothing to do

  case 'fcdc_buffer'
    [host, port] = filetype_check_uri(filename);
    buffer('flush_evt', [], host, port);

  case 'fcdc_mysql'
    % open the database
    [user, password, server, port] = filetype_check_uri(filename);
    if ~isempty(port)
      server = sprintf('%s:%d', server, port);
    end
    mysql('open', server, user, password);
    % remove all previous event information
    cmd = 'TRUNCATE TABLE fieldtrip.event';
    mysql(cmd);
    mysql('close');

  case 'matlab'
    if exist(filename, 'file')
      warning('deleting existing file ''%s''', filename);
      delete(filename);
    end

  otherwise
    error('unsupported data format');
end

author: PatternRecognition

repo: OpenBMI


text:
function denu=tide_oload(tut,odisp)

D2R=pi/180;
args=[1.40519E-4, 2.0,-2.0, 0.0, 0.00;  % M2
      1.45444E-4, 0.0, 0.0, 0.0, 0.00;  % S2
      1.37880E-4, 2.0,-3.0, 1.0, 0.00;  % N2
      1.45842E-4, 2.0, 0.0, 0.0, 0.00;  % K2
      0.72921E-4, 1.0, 0.0, 0.0, 0.25;  % K1
      0.67598E-4, 1.0,-2.0, 0.0,-0.25;  % O1
      0.72523E-4,-1.0, 0.0, 0.0,-0.25;  % P1
      0.64959E-4, 1.0,-3.0, 1.0,-0.25;  % Q1
      0.53234E-5, 0.0, 2.0, 0.0, 0.00;  % Mf
      0.26392E-5, 0.0, 1.0,-1.0, 0.00;  % Mm
      0.03982E-5, 2.0, 0.0, 0.0, 0.00];  % Ssa
ep1975=[1975,1,1,0,0,0];
dp=zeros(3,1);

%angular argument
ep=time2epoch(tut);
fday=ep(4)*3600.0+ep(5)*60.0+ep(6);
ep(4)=0;ep(5)=0;ep(6)=0.0;
days=timediff(epoch2time(ep),epoch2time(ep1975))/86400.0+1.0;
t=(27392.500528+1.000000035*days)/36525.0;
t2=t*t; t3=t2*t;
a(1)=fday;
a(2)=(279.69668+36000.768930485*t+3.03E-4*t2)*D2R;% H0
a(3)=(270.434358+481267.88314137*t-0.001133*t2+1.9E-6*t3)*D2R;% S0
a(4)=(334.329653+4069.0340329577*t-0.010325*t2-1.2E-5*t3)*D2R;% P0
a(5)=2.0*pi;

%dispalcements by 11 constituents
for i=1:11
    ang=0;
    for j=1:5
        ang=ang+a(j)*args(i,j);
    end
    for j=1:3
        dp(j)=dp(j)+odisp(i,j)*cos(ang-odisp(i,j+3)*D2R);
    end
end

denu(1)=-dp(2);
denu(2)=-dp(3);
denu(3)=dp(1);

return

repo: GINav


text:
function planC = removeItemsFromDICOM(planC)
% "removeItemsFromDICOM"
%   Reformats the DICOM planC data structure so that struct arrays
%   consisting of fields Item_1, Item_2, ... Item_n are concatenated
%   into a struct array 1...n under the same parent.
%   See help condenseFields for details on how this is done.
%
%   JRA 11/3/03
%
% Usage: planC = removeItemsFromDICOM(planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
planC = condenseFields(planC, 'Item_\d', 'N/A');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [center, Cn, PNR] = initComponents_residual_parallel(obj, K, save_avi, use_parallel, min_corr, min_pnr, seed_method)
%% initializing spatial/temporal components for the residual video
%% input:
%   K:  scalar, maximum number of neurons
%   frame_range: 1 X 2 vector indicating the starting and ending frames
%   save_avi: save the video of initialization procedure
%   use_parallel: boolean, do initialization in patch mode or not.
%       default(true); we recommend you to set it false only when you want to debug the code.

%% Output:
%   center: d*2 matrix, centers of all initialized neurons.
%   Cn:     correlation image
%   PNR:    peak to noise ratio
%% Author: Pengcheng Zhou, Columbia University, 2017
%% email: zhoupc1988@gmail.com

%% process parameters

try
    % map data
    mat_data = obj.P.mat_data;

    % folders and files for saving the results
    log_folder = obj.P.log_folder;
    log_file = obj.P.log_file;
    log_data_file = obj.P.log_data;
    log_data = matfile(log_data_file, 'Writable', true); %#ok<NASGU>

    % dimension of data
    dims = mat_data.dims;
    d1 = dims(1);
    d2 = dims(2);
    obj.options.d1 = d1;
    obj.options.d2 = d2;

    % parameters for patching information
    patch_pos = mat_data.patch_pos;
    block_pos = mat_data.block_pos;

    % number of patches
    [nr_patch, nc_patch] = size(patch_pos);
catch
    error('No data file selected');
end
fprintf('\n---------------PICK NEURONS FROM THE RESIDUAL----------------\n');

if exist('min_corr', 'var') && ~isempty(min_corr)
    obj.options.min_corr = min_corr;
end
if exist('min_pnr', 'var') && ~isempty(min_pnr)
    obj.options.min_pnr = min_pnr;
end
if exist('seed_method', 'var') && ~isempty(seed_method)
    obj.options.seed_method = seed_method;
end
% frames to be loaded for initialization
frame_range = obj.frame_range;
T = diff(frame_range) + 1;

% maximum neuron number in each patch
if (~exist('K', 'var')) || (isempty(K))
    % if K is not specified, use a very large number as default
    K = round((d1*d2));
end

% exporting initialization procedures as a video
if ~exist('save_avi', 'var')||isempty(save_avi)
    save_avi = false; %don't save initialization procedure
elseif save_avi
    use_parallel = false;
end

% use parallel or not
if ~exist('use_parallel', 'var')||isempty(use_parallel)
    use_parallel = true; %don't save initialization procedure
end

% parameter for avoiding using boundaries pixels as seed pixels
options = obj.options;
if ~isfield(options, 'bd') || isempty(options.bd')
    options.bd = options.gSiz;   % boundary pixesl to be ignored during the process of detecting seed pixels
end
bd = options.bd;

if strcmpi(obj.options.seed_method, 'manual')
    use_parallel = false;
end
% no centering of the raw video
% options.center_psf = false;

%% prepare for the variables for computing the background.
bg_model = obj.options.background_model;
bg_ssub = obj.options.bg_ssub;
W = obj.W;
b0 = obj.b0;
b = obj.b;
f = obj.f;

%% the extracted neurons' signals
A = cell(nr_patch, nc_patch);
C = cell(nr_patch, nc_patch);
sn = cell(nr_patch, nc_patch);
ind_neurons = cell(nr_patch, nc_patch);

AA = cell(nr_patch, nc_patch);   % save the ai^T*ai for each neuron

for mpatch=1:(nr_patch*nc_patch)
    if strcmpi(bg_model, 'ring')
        tmp_block = block_pos{mpatch};
    else
        tmp_block = patch_pos{mpatch};
    end
    % find the neurons that are within the block
    mask = zeros(d1, d2);
    mask(tmp_block(1):tmp_block(2), tmp_block(3):tmp_block(4)) = 1;
    ind = (reshape(mask(:), 1, [])* obj.A>0);
    A{mpatch}= obj.A(logical(mask), ind);
    AA{mpatch}= sum(A{mpatch}.^2, 1);
    sn{mpatch} = obj.P.sn(logical(mask));
    C{mpatch} = obj.C(ind, :);
    ind_neurons{mpatch} = find(ind);    % indices of the neurons within each patch
end

%% start initialization
% save the log infomation
k_options = obj.P.k_options +1;
eval(sprintf('log_data.options_%d=options; ', k_options));
obj.P.k_options = k_options;

flog = fopen(log_file, 'a');
fprintf(flog, '[%s]\b', get_minute());
fprintf(flog, 'Start picking neurons from the residual......\n\tThe collection of options are saved as intermediate_results.options_%d\n', k_options);

Ain = cell(nr_patch, nc_patch); % save spatial components of neurons in each patch
Cin = cell(nr_patch, nc_patch); % save temporal components of neurons in each patch, denoised trace
Sin = cell(nr_patch, nc_patch); % save temporal components of neurons in each patch, deconvolved trace
Cin_raw = cell(nr_patch, nc_patch); % save temporal components of neurons in each patch, raw trace
kernel_pars = cell(nr_patch, nc_patch); % save temporal components of neurons in each patch
center = cell(nr_patch, nc_patch);     % save centers of all initialized neurons
Cn = zeros(d1, d2);
PNR = zeros(d1, d2);
default_kernel = obj.kernel;

results = cell(nr_patch*nc_patch, 1);
if use_parallel
    parfor mpatch=1:(nr_patch*nc_patch)
        % get the indices corresponding to the selected patch
        tmp_patch = patch_pos{mpatch};
        if strcmpi(bg_model, 'ring')
            % when it's ring model, log pixels surrounding the patch
            tmp_block = block_pos{mpatch};
        else
            tmp_block = tmp_patch;
        end
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);

        % use ind_patch to indicate pixels within the patch
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % get the neural activity
        C_patch = C{mpatch};                % previous estimation of neural activity
        if isempty(C_patch)
            C_patch = 0;
            A_patch = 0;
        else
            A_patch = A{mpatch};
        end

        % boundaries pixels to be avoided for detecting seed pixels
        tmp_options = options;
        tmp_options.visible_off = true;
        tmp_options.bd = bd*([r==1, r==nr_patch, c==1, c==nc_patch]);

        % patch dimension
        tmp_options.d1 = diff(tmp_patch(1:2))+1;
        tmp_options.d2 = diff(tmp_patch(3:4))+1;

        % parameter for calcium indicators. This one may not be used if the
        % selected deconvolution algorithm doesn't need it
        tmp_options.kernel = default_kernel;

        % file names for saving avi file
        if save_avi
            tmp_save_avi = sprintf('%sinitialization_res_%d_%d_%d_%d.avi', log_folder, tmp_block(1), tmp_block(2), tmp_block(3), tmp_block(4));
        else
            tmp_save_avi = save_avi;
        end
        % get data
        if strcmpi(bg_model, 'ring')
            % including areas outside of the patch for recorving background
            % in the ring model
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        else
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, false);
        end
        [nr_block, nc_block, ~] = size(Ypatch);

        Ypatch = double(reshape(Ypatch, [], T)) - A_patch*C_patch;
        % get background
        if strcmpi(bg_model, 'ring')
            W_ring = W{mpatch};
            b0_ring = b0{mpatch};
            %             Ypatch = bsxfun(@minus, Ypatch(ind_patch,:)- W_ring*Ypatch, b0_ring-W_ring*mean(Ypatch, 2));
            if bg_ssub==1
                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch,:))- W_ring*Ypatch, b0_ring-W_ring*mean(Ypatch, 2));
            else
                % get the dimension of the downsampled data
                [d1s, d2s] = size(imresize(zeros(nr_block, nc_block), 1/bg_ssub));
                % downsample data and reconstruct B^f
                temp = reshape(bsxfun(@minus, Ypatch, mean(Ypatch, 2)), nr_block, nc_block, []);
                temp = imresize(temp, 1./bg_ssub);
                Bf = reshape(W_ring*reshape(temp, [], T), d1s, d2s, T);
                Bf = imresize(Bf, [nr_block, nc_block]);
                Bf = reshape(Bf, [], T);

                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch, :)) - Bf(ind_patch, :), b0_ring);
            end
        elseif strcmpi(bg_model, 'nmf')
            b_nmf = b{mpatch};
            f_nmf = f{mpatch};
            Ypatch = Ypatch- b_nmf*f_nmf;
        else
            b_svd = b{mpatch};
            f_svd = f{mpatch};
            b0_svd = b0{mpatch};
            Ypatch = Ypatch - bsxfun(@plus, b_svd*f_svd, b0_svd);
        end

        [tmp_results, tmp_center, tmp_Cn, tmp_PNR, ~] = greedyROI_endoscope(Ypatch, K, tmp_options, [], tmp_save_avi);

        % put everthing into one struct variable
        tmp_results.center = tmp_center;
        tmp_results.Cn = tmp_Cn;
        tmp_results.PNR = tmp_PNR;
        results{mpatch} = tmp_results;
        %     eval(sprintf('results_patch_%d=tmp_results;', mpatch));  %#ok<PFBFN>
        % display initialization progress
        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
    end
else
    for mpatch=1:(nr_patch*nc_patch)
        % get the indices corresponding to the selected patch
        tmp_patch = patch_pos{mpatch};
        if strcmpi(bg_model, 'ring')
            % when it's ring model, log pixels surrounding the patch
            tmp_block = block_pos{mpatch};
        else
            tmp_block = tmp_patch;
        end
        [r, c] = ind2sub([nr_patch, nc_patch], mpatch);

        % use ind_patch to indicate pixels within the patch
        ind_patch = false(diff(tmp_block(1:2))+1, diff(tmp_block(3:4))+1);
        ind_patch((tmp_patch(1):tmp_patch(2))-tmp_block(1)+1, (tmp_patch(3):tmp_patch(4))-tmp_block(3)+1) = true;

        % get the neural activity
        C_patch = C{mpatch};                % previous estimation of neural activity
        if isempty(C_patch)
            C_patch = 0;
            A_patch = 0;
        else
            A_patch = A{mpatch};
        end

        % boundaries pixels to be avoided for detecting seed pixels
        tmp_options = options;
        tmp_options.visible_off = true;
        tmp_options.bd = bd*([r==1, r==nr_patch, c==1, c==nc_patch]);

        % patch dimension
        tmp_options.d1 = diff(tmp_patch(1:2))+1;
        tmp_options.d2 = diff(tmp_patch(3:4))+1;

        % parameter for calcium indicators. This one may not be used if the
        % selected deconvolution algorithm doesn't need it
        tmp_options.kernel = default_kernel;

        % file names for saving avi file
        if save_avi
            tmp_save_avi = sprintf('%sinitialization_res_%d_%d_%d_%d.avi', log_folder, tmp_block(1), tmp_block(2), tmp_block(3), tmp_block(4));
        else
            tmp_save_avi = save_avi;
        end
        % get data
        if strcmpi(bg_model, 'ring')
            % including areas outside of the patch for recorving background
            % in the ring model
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
        else
            Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, false);
        end
        [nr_block, nc_block, ~] = size(Ypatch);

        Ypatch = double(reshape(Ypatch, [], T)) - A_patch*C_patch;
        % get background
        if strcmpi(bg_model, 'ring')
            W_ring = W{mpatch};
            b0_ring = b0{mpatch};
            %             Ypatch = bsxfun(@minus, Ypatch(ind_patch,:)- W_ring*Ypatch, b0_ring-W_ring*mean(Ypatch, 2));
            if bg_ssub==1
                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch,:))- W_ring*Ypatch, b0_ring-W_ring*mean(Ypatch, 2));
            else
                % get the dimension of the downsampled data
                [d1s, d2s] = size(imresize(zeros(nr_block, nc_block), 1/bg_ssub));
                % downsample data and reconstruct B^f
                temp = reshape(bsxfun(@minus, Ypatch, mean(Ypatch, 2)), nr_block, nc_block, []);
                temp = imresize(temp, 1./bg_ssub);
                Bf = reshape(W_ring*reshape(temp, [], T), d1s, d2s, T);
                Bf = imresize(Bf, [nr_block, nc_block]);
                Bf = reshape(Bf, [], T);

                Ypatch = bsxfun(@minus, double(Ypatch(ind_patch, :)) - Bf(ind_patch, :), b0_ring);
            end

        elseif strcmpi(bg_model, 'nmf')
            b_nmf = b{mpatch};
            f_nmf = f{mpatch};
            Ypatch = Ypatch- b_nmf*f_nmf;
        else
            b_svd = b{mpatch};
            f_svd = f{mpatch};
            b0_svd = b0{mpatch};
            Ypatch = Ypatch - bsxfun(@plus, b_svd*f_svd, b0_svd);
        end

        [tmp_results, tmp_center, tmp_Cn, tmp_PNR, ~] = greedyROI_endoscope(Ypatch, K, tmp_options, [], tmp_save_avi);

        % put everthing into one struct variable
        tmp_results.center = tmp_center;
        tmp_results.Cn = tmp_Cn;
        tmp_results.PNR = tmp_PNR;
        results{mpatch} = tmp_results;
        %     eval(sprintf('results_patch_%d=tmp_results;', mpatch));  %#ok<PFBFN>
        % display initialization progress
        fprintf('Patch (%2d, %2d) is done. %2d X %2d patches in total. \n', r, c, nr_patch, nc_patch);
    end
end

%% collect results
for mpatch=1:(nr_patch*nc_patch)
    % get the indices corresponding to the selected patch
    [mr, mc] = ind2sub([nr_patch, nc_patch], mpatch);
    tmp_patch = patch_pos{mr, mc};
    tmp_block = tmp_patch;
    r0 = tmp_block(1);
    r1 = tmp_block(2);
    c0 = tmp_block(3);
    c1 = tmp_block(4);
    %         ind_patch = true(r1-r0+1, c1-c0+1);
    %         ind_patch((tmp_patch(1):tmp_patch(2))-r0+1, (tmp_patch(3):tmp_patch(4))-c0+1) = false;
    %
    % unpack results
    tmp_results = results{mpatch};
    if isempty(tmp_results)
        continue;
    end

    tmp_Ain = tmp_results.Ain;
    tmp_ind = (sum(tmp_Ain, 1)>0);
    tmp_Ain = tmp_Ain(:, tmp_ind);
    tmp_Cin = tmp_results.Cin(tmp_ind,:);
    tmp_Cin_raw = tmp_results.Cin_raw(tmp_ind,:);
    tmp_center = tmp_results.center(tmp_ind,:) ;
    tmp_Cn = tmp_results.Cn;
    tmp_PNR = tmp_results.PNR;
    if options.deconv_flag
        tmp_Sin = tmp_results.Sin(tmp_ind,:);
        tmp_kernel_pars = tmp_results.kernel_pars(tmp_ind,:);
    end
    tmp_K = size(tmp_Ain, 2);   % number of neurons within the selected patch
    [tmp_d1, tmp_d2] = size(tmp_Cn);

    temp = zeros(d1, d2, tmp_K);  % spatial components of all neurons
    temp(r0:r1, c0:c1, :) = reshape(full(tmp_Ain), tmp_d1, tmp_d2, []);
    Ain{mr, mc} = reshape(temp, d1*d2, tmp_K);
    Cin{mr, mc} = tmp_Cin;      % temporal components of all neurons
    Cin_raw{mr, mc} = tmp_Cin_raw;
    if options.deconv_flag
        Sin{mr, mc} = tmp_Sin;
        kernel_pars{mr,mc} = tmp_kernel_pars;
    end
    center{mr, mc} = bsxfun(@plus, tmp_center, [r0-1, c0-1]); % centers

    Cn(r0:r1, c0:c1) = max(Cn(r0:r1, c0:c1), tmp_Cn);
    PNR(r0:r1, c0:c1) = max(PNR(r0:r1, c0:c1), tmp_PNR);
end
%% export the results
Ain = cell2mat(reshape(Ain, 1, []));
Cin = cell2mat(reshape(Cin, [], 1));
Cin_raw = cell2mat(reshape(Cin_raw, [], 1));
center = cell2mat(reshape(center, [], 1));
obj.A = [obj.A, sparse(Ain)];
obj.C = [obj.C; Cin];
obj.C_raw = [obj.C_raw; Cin_raw];
if options.deconv_flag
    Sin = cell2mat(reshape(Sin, [], 1));
    kernel_pars = cell2mat(reshape(kernel_pars, [], 1));
    obj.S = [obj.S; sparse(Sin)];
    obj.P.kernel_pars = [obj.P.kernel_pars; kernel_pars];
else
    obj.S = [obj.S; zeros(size(obj.C))];
end
K = size(Ain, 2);
k_ids = obj.P.k_ids;
obj.ids = [obj.ids, k_ids+(1:K)];
obj.tags = [obj.tags;zeros(K,1, 'like', uint16(0))];
obj.P.k_ids = K+k_ids;

%% save the results to log

fprintf(flog, '[%s]\b', get_minute());
fprintf(flog, 'Finished the initialization of neurons from the residual video.\n');
fprintf(flog, '\tIn total, %d neurons were detected. \n', size(Ain,2));

if obj.options.save_intermediate
    initialization_res.Ain = sparse(Ain);
    initialization_res.Cin = Cin;
    initialization_res.Cin_raw = Cin_raw;
    initialization_res.options = options;  %#ok<STRNU>

    tmp_str = get_date();
    tmp_str=strrep(tmp_str, '-', '_');
    eval(sprintf('log_data.initialization_res_%s = initialization_res;', tmp_str));

    fprintf(flog, '\tThe  results were saved as intermediate_results.initialization_res_%s\n\n', tmp_str);
end
fclose(flog);

author: zhoupc


text:
function formats = supported_wc_extensions()
    % Return a cell of arrays with the formats supported. Automatically updated
    formats ={};
    raw_readers_folder = fileparts(mfilename('fullpath'));
    files = dir(raw_readers_folder);
    for i = 1: length(files)
        fname = files(i).name;
        fend = regexp(fname,'_wc_reader.m');
        if ~isempty (fend)
             formats = [formats,{fname(1:fend-1)}];
        end
    end
end

author: csn-le


text:
function twRatio = computeTWRatio(throttle, ut, rVect, vVect, tankMasses, dryMass, stgStates, lvState, tankStates, bodyInfo, storageSoCs, powerStorageStates)
    altitude = norm(rVect) - bodyInfo.radius;
    presskPa = getPressureAtAltitude(bodyInfo, altitude);

    [~, totalThrust]= LaunchVehicleStateLogEntry.getTankMassFlowRatesDueToEngines(tankStates, tankMasses, stgStates, throttle, lvState, presskPa, ut, rVect, vVect, bodyInfo, [], storageSoCs, powerStorageStates, []);

    totalMass = (dryMass + sum(tankMasses))*1000; %kg
    totalThrust = totalThrust * 1000; % N

%     twRatio = computeSLThrustToWeight(bodyInfo, totalThrust, totalMass);
    twRatio = computeTrueThrustToWeight(bodyInfo, totalThrust, totalMass, altitude);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%IOPEN Morphological opening
%
% OUT = IOPEN(IM, SE, OPTIONS) is the image IM after morphological opening
% with the structuring element SE.  This is a morphological erosion followed
% by dilation.
%
% OUT = IOPEN(IM, SE, N, OPTIONS) as above but the structuring element
% SE is applied N times, that is N erosions followed by N dilations.
%
% Notes::
% - For binary image an opening operation can be used to eliminate small white
%   noise regions.
% - It is cheaper to apply a smaller structuring element multiple times than
%   one large one, the effective structuring element is the Minkowski sum
%   of the structuring element with itself N times.
% - Windowing options of IMORPH can be passed.  By default output image is
%   same size as input image.
%
% See also ICLOSE, IDILATE, IERODE, IMORPH.

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

function b = iopen(a, se, varargin)

    b = ierode(a, se, varargin{:});
    b = idilate(b, se, varargin{:});

author: petercorke

repo: machinevision-toolbox-matlab


text:
function M=minMatOverDim(C,minIdx)
%%MINMATOVERDIM Given a multidimensional matrix C, this function returns
%           matrix M that is obtained by minimizing C over the given
%           dimension and removing the given dimension. This is equivalent
%           to using a call to min(C,[],minIdx) in Matlab. However, here,
%           we do it without using any matrix/ vector commands so that it
%           can mirror how one might implement such a function in C.
%
%INPUTS: C A real n1Xn2Xn3..XnS S-dimensional hypermatrix.
%   minIdx The real index from 1 to S over which the minimization should be
%          performed.
%
%OUTPUTS: M An n1X...n(minIdx-1)X1Xn(minIdx+1)X...nS matrix holding the
%           minimum values over the specified dimension.
%
%The algorithm consists of two types of steps. First, a step in the linear
%indexation of C to go from one value of an index in spot minIdx is
%prod(nVals(1:(minIdx-1))), where nVals=size(C). Thus, we can go from one
%element to the next for the minimization. For a 3D matrix, the step to go
%from C(i1,1,i2) to C(i1+1,1,i2) increments by 1. However, to go from
%C(n1,1,i2) to C(1,1,i2+1) is a step of prod(nVals(1:minIdx)). The same
%type of thing applies to a matrix with more dimensions, since all
%dimensions before the one in question and after can be collapsed into 1
%dimension. Thus, this function puts the above rules together to go
%minimize the matrix.
%
%EXAMPLE:
%Here, we just show that the results are equivalent to using the min
%command with a resize.
% C=randn(13,18,11,6,9);
% minIdx=3;
% M=minMatOverDim(C,minIdx);
% MAlt=min(C,[],minIdx);
% all(M(:)==MAlt(:))
%The result is 1, indicating that the two values are equal.
%
%March 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

nVals=size(C);
S=numel(nVals);

if(minIdx<1)
    error('Invalid value of minIdx given.')
end

if(minIdx>S)
    M=C;
    return;
end

MDims=nVals;
MDims(minIdx)=1;
M=zeros(MDims);

totalNumElsM=prod(MDims);

%Note that prod([])=1
incrMinIdx=prod(nVals(1:(minIdx-1)));

incrBigStep=incrMinIdx*nVals(minIdx);

CStartIdx=0;
curMCol=0;
for curEl=0:(totalNumElsM-1)
    curIdx=CStartIdx;
    minVal=C(curIdx+1);

    curIdx=curIdx+incrMinIdx;
    for i=2:nVals(minIdx)
        curVal=C(curIdx+1);

        if(curVal<minVal)
            minVal=curVal;
        end

        curIdx=curIdx+incrMinIdx;
    end
    M(curEl+1)=minVal;

    %If a big step has to be taken.
    if(mod(curEl+1,incrMinIdx)==0)
        curMCol=curMCol+1;

        CStartIdx=incrBigStep*curMCol;
    else
        CStartIdx=CStartIdx+1;
    end
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


%  This m-file calculates the area of a circle,
%  and displays the result.
radius = 2.5;
area = pi * 2.5^2;
string = ['The area of the circle is ' num2str(area)];
disp(string);


text:
%RANDINIT Reset random number generator
%
% RANDINIT resets the defaul random number stream.  For example:
%
%         >> rand
%         ans =
%            0.8147
%         >> rand
%         ans =
%             0.9058
%         >> rand
%         ans =
%             0.1270
%         >> randinit
%         >> rand
%         ans =
%             0.8147

%
% See also RandStream.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function randinit(seed)

    stream = RandStream.getGlobalStream;
    stream.reset()

repo: spatialmath-matlab


text:
function all_ok_ = t_run_tests(test_names, verbose)
%T_RUN_TESTS  Run a series of tests.
%   ALL_OK = T_RUN_TESTS(TEST_NAMES, VERBOSE)
%
%   Runs a set of tests whose names are given in the cell array TEST_NAMES.
%   If the optional parameter VERBOSE is true, it prints the details of the
%   individual tests. Optionally returns an ALL_OK flag, equal to 1 if all
%   tests pass (and the number matches the expected number), 0 otherwise.
%
%   Example:
%       tests{end+1} = 't_loadcase';
%       tests{end+1} = 't_jacobian';
%       tests{end+1} = 't_hessian';
%       t_run_tests( tests, verbose );
%
%   See also T_BEGIN, T_END.

%   MP-Test
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Test.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mptest for more info.

if nargin < 2
    verbose = 0;
end

global t_num_of_tests;
global t_counter;
global t_ok_cnt;
global t_not_ok_cnt;
global t_skip_cnt;

%% figure out padding for printing
if ~verbose
    len = zeros(length(test_names), 1);
    for k = 1:length(test_names)
        len(k) = length(test_names{k});
    end
    maxlen = max(len);
end

%% initialize statistics
num_of_tests = 0;
counter = 0;
ok_cnt = 0;
not_ok_cnt = 0;
skip_cnt = 0;

t0 = tic;
for k = 1:length(test_names)
    if verbose
        fprintf('\n----------  %s  ----------\n', test_names{k});
    else
        pad = maxlen + 4 - length(test_names{k});
        fprintf('%s', test_names{k});
        for m = 1:pad, fprintf('.'); end
    end
    feval( test_names{k}, ~verbose );

    num_of_tests    = num_of_tests  + t_num_of_tests;
    counter         = counter       + t_counter;
    ok_cnt          = ok_cnt        + t_ok_cnt;
    not_ok_cnt      = not_ok_cnt    + t_not_ok_cnt;
    skip_cnt        = skip_cnt      + t_skip_cnt;
end

if verbose
    fprintf('\n\n----------  Summary  ----------\n');
end
if counter == num_of_tests && counter == ok_cnt + skip_cnt && not_ok_cnt == 0
    all_ok = 1;
    if skip_cnt
        fprintf('All tests successful (%d passed, %d skipped of %d)', ...
            ok_cnt, skip_cnt, num_of_tests);
    else
        fprintf('All tests successful (%d of %d)', ok_cnt, num_of_tests);
    end
else
    all_ok = 0;
    fprintf('Ran %d of %d tests: %d passed, %d failed', ...
        counter, num_of_tests, ok_cnt, not_ok_cnt);
    if skip_cnt
        fprintf(', %d skipped', skip_cnt);
    end
end
fprintf('\nElapsed time %.2f seconds.\n', toc(t0));

if nargout
    all_ok_ = all_ok;   %% copy to optional output arg
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower