text:
% performs argument checks/transformations similar to those found in
% gp.m from GPML, but Hessian friendly

% Copyright (c) 2014--2015 Roman Garnett.
% Modified 2016 Luigi Acerbi.

% function gpstruct = gpset(gpstruct, hyp, inf, mean, cov, lik, prior, bounds, x, y)
function gpstruct = gpset(gpstruct, varargin)

% Input arguments
args = {'hyp', 'inf', 'mean', 'cov', 'lik', 'prior', 'bounds'};

for idx = 1:numel(varargin)
    % Inputs are x and y
    if isnumeric(varargin{idx})
        if numel(varargin) ~= idx+1 || ~isnumeric(varargin{idx+1})
            error('gpstruct:trainingdata', ...
                'Training data X and Y need to be the last two inputs of GPSET.');
        end
        gpstruct.x = varargin{idx};
        gpstruct.y = varargin{idx+1};
        break;
    end

    if idx > numel(varargin)
        error('gpstruct:inputarguments', ...
            'Too many input arguments.');
    end

    % Assign input to struct if different than skip placeholder
    if ~strcmp(varargin{idx}, '~')
        gpstruct.(args{idx}) = varargin{idx};
    end
end

% Default gp values
defaultstruct = struct( ...
    'hyp', struct('cov',[],'lik',[],'mean',[]), ...
    'inf',{@exact_inference}, ...
    'cov', {@isotropic_sqdexp_covariance}, ...
    'lik',{@likGauss}, ...
    'mean', {@zero_mean}, ...
    'x', [], ...
    'y', [], ...
    's', [], ...
    'prior', [], ...
    'bounds', [] ...
    );

% Assign default values
for field = fields(defaultstruct)'
    if ~isfield(gpstruct,field{:}) || isempty(gpstruct.(field{:}))
        gpstruct.(field{:}) = defaultstruct.(field{:});
    end
end

% Allow string/function handle input; convert to cell arrays if necessary
for field = {'inf', 'cov', 'lik', 'mean'}
    if ischar(gpstruct.(field{:})) || isa(gpstruct.(field{:}), 'function_handle')
        gpstruct.(field{:}) = {gpstruct.(field{:})};
    end
end

% If HYP is a cell array, convert to struct array
if iscell(gpstruct.hyp) && isstruct(gpstruct.hyp{1})
    hyptemp = struct('cov','lik','mean');
    for i = 1:length(gpstruct.hyp)
        hyptemp(i) = gpstruct.hyp{i};
    end
    gpstruct.hyp = hyptemp;
end

% Ensure all hyperparameter fields exist
for field = {'cov', 'lik', 'mean'}
    for i = 1:numel(gpstruct.hyp)
        if (~isfield(gpstruct.hyp(i), field{:}))
            gpstruct.hyp(i).(field{:}) = [];
        end
    end
end

% Check dimension of hyperparameter fields if X is available
D = size(gpstruct.x, 2);
if D > 0
    shortfields = {'cov','lik','mean'};
    longfields = {'covariance','likelihood','mean'};
    for ifield = 1:numel(shortfields)
        field = shortfields{ifield};
        expression = feval(gpstruct.(field){:});
        for i = 1:numel(gpstruct.hyp)
            nparams = numel(gpstruct.hyp(i).(field));
            if (nparams ~= eval(expression))
                error('gpstruct:incorrect_specification', ...
                      'wrong number of %s hyperparameters (%i given, %s expected)', ...
                      longfields{ifield}, ...
                      nparams, ...
                      expression);
            end
        end
    end
end

% Check prior and bounds
for field = {'cov', 'lik', 'mean'}
    % Ensure all prior fields exist
    if (~isfield(gpstruct.prior, field{:}))
        gpstruct.prior.(field{:}) = [];
    end
    % Ensure all bounds fields exist
    if (~isfield(gpstruct.bounds, field{:}))
        gpstruct.bounds.(field{:}) = [];
    end
    % Assign [-Inf;Inf] bounds by default
    num_entries = numel(gpstruct.hyp(1).(field{:}));
    for i = 1:num_entries
        if numel(gpstruct.prior.(field{:})) < i
            gpstruct.prior.(field{:}){i} = [];
        end
        if numel(gpstruct.bounds.(field{:})) < i || ...
                isempty(gpstruct.bounds.(field{:}){i})
            gpstruct.bounds.(field{:}){i} = [-Inf,Inf];
        end
    end
end

% Add field with (empty) posterior
if ~isfield(gpstruct,'post'); gpstruct.post = []; end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 25
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% setup for figure 25.3

% generate a random signal of 3 seconds
srate = 1000;

randsig1 = randn(1,3*srate); % 3 seconds, with this sampling rate
randsig2 = randn(1,3*srate);

% now filter at 5 Hz
f       = 5; % frequency of wavelet in Hz
time    = -1:1/srate:1; % time for wavelet, from -1 to 1 second in steps of 1/sampling-rate
s       = 6/(2*pi*f); % width of Gaussian
wavelet = exp(2*pi*1i*f.*time) .* exp(-time.^2./(2*s^2));

% FFT parameters
n_wavelet            = length(wavelet);
n_data               = length(randsig1);
n_convolution        = n_wavelet+n_data-1;
half_of_wavelet_size = (length(wavelet)-1)/2;

% FFT of wavelet and EEG data
convolution_result_fft = ifft(fft(wavelet,n_convolution).*fft(randsig1,n_convolution),n_convolution)*sqrt(s)/10;
filtsig1  = real(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size));
anglesig1 = angle(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size));

convolution_result_fft = ifft(fft(wavelet,n_convolution).*fft(randsig2,n_convolution),n_convolution)*sqrt(s)/10;
filtsig2  = real(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size));
anglesig2 = angle(convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size));

figure
for i=1:2
    subplot(2,1,i)
    eval([ 'plot(randsig' num2str(i) ')' ]) % "eval" can be useful for increasing control over commands
    hold on
    eval([ 'plot(filtsig' num2str(i) ',''r'')' ])
end

%% Figure 25.3

% initialize output correlation matrix
correlations = zeros(5,round(1000/f));

for i=1:round(1000/f)

    % correlation of unfiltered random signal
    temp = corrcoef(randsig1(1:end-i),randsig1(i+1:end));
    correlations(1,i) = temp(1,2);

    % correlation of filtered signal
    temp = corrcoef(filtsig1(1:end-i),filtsig1(i+1:end));
    correlations(2,i) = temp(1,2);

    % phase clustering
    correlations(3,i) = abs(mean(exp(1i*( angle(anglesig1(1:end-i)-anglesig1(i+1:end))))));

    % difference of correlations of filtered signal
    temp = corrcoef(filtsig2(1:end-i),filtsig2(i+1:end));
    correlations(4,i) = temp(1,2) - correlations(2,i);

    % difference of phase clusterings
    correlations(5,i) = abs(mean(exp(1i*( angle(anglesig2(1:end-i)-anglesig2(i+1:end)))))) - correlations(3,i);
end

figure
plot(correlations')
xlabel('Lag (ms)')
ylabel('Connectivity strength')
legend({'unfiltered';'power corr';'ISPC';'corr diffs';'ISPC diffs'})

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function roidb = roidb_from_voc(imdb)
% roidb = roidb_from_voc(imdb)
%   Builds an regions of interest database from imdb image
%   database. Uses precomputed selective search boxes available
%   in the R-CNN data package.
%
%   Inspired by Andrea Vedaldi's MKL imdb and roidb code.

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Change by Ning Zhang

cache_file = ['./imdb/cache/roidb_' imdb.name];
try
  load(cache_file);
catch
  roidb.name = imdb.name;

  fprintf('Loading region proposals...');

  regions_file = sprintf('./data/selective_search_data/%s', roidb.name);
  % follow the selective search instructions and run selective search on bird images.
  try
    regions = load(regions_file);
  catch
    disp('Run selective search on cub2011 images first.');
    exit(-1);
  end
  fprintf('done\n');
  % if config file doesn't exist. run config = get_bird_data;
  load('caches/cub2011_config.mat');
  for i = 1:length(imdb.image_ids)
    tic_toc_print('roidb (%s): %d/%d\n', roidb.name, i, length(imdb.image_ids));
    voc_rec = get_pascal_format(config, imdb.image_ids{i});
    roidb.rois(i) = attach_proposals(voc_rec, regions.boxes{i}, imdb.class_to_id);
  end

  fprintf('Saving roidb to cache...');
  save(cache_file, 'roidb', '-v7.3');
  fprintf('done\n');
end

function voc_rec = get_pascal_format(config, img_path)
  % find image id for img_path
  full_path = [config.img_base img_path];
  train_idx = find(ismember(config.impathtrain, full_path));
  test_idx = find(ismember(config.impathtest, full_path));
  assert(isempth(train_idx) || isempty(test_idx));
  assert(~isempty(train_idx) || ~isempty(test_idx));
  if ~isempty(train_idx)
    voc_rec.object(1).bbox = config.train_box{1}(train_idx,:);
    voc_rec.object(1).class = 'bbox';
    count = 2;
    if config.train_box{2}(train_idx,1) ~= -1
      voc_rec.object(count).bbox = config.train_box{2}(train_idx,:);
      voc_rec.object(count).class = 'head';
      count = count + 1;
    end
    if config.train_box{3}(train_idx,1) ~= -1
      voc_rec.object(count).bbox = config.train_box{3}(train_idx,:);
      voc_rec.object(count).class = 'body';
    end
  end
  if ~isempty(test_idx)
    voc_rec.object(1).bbox = config.test_box{1}(test_idx,:);
    voc_rec.object(1).class = 'bbox';
    count = 2;
    if config.test_box{2}(test_idx,1) ~= -1
      voc_rec.object(count).bbox = config.test_box{2}(test_idx,:);
      voc_rec.object(count).class = 'head';
      count = count + 1;
    end
    if config.test_box{3}(test_idx,1) ~= -1
      voc_rec.object(count).bbox = config.test_box{3}(test_idx,:);
      voc_rec.object(count).class = 'body';
    end
  end
end

% ------------------------------------------------------------------------
function rec = attach_proposals(voc_rec, boxes, class_to_id)
% ------------------------------------------------------------------------

% change selective search order from [y1 x1 y2 x2] to [x1 y1 x2 y2]
boxes = boxes(:, [2 1 4 3]);

%           gt: [2108x1 double]
%      overlap: [2108x20 single]
%      dataset: 'voc_2007_trainval'
%        boxes: [2108x4 single]
%         feat: [2108x9216 single]
%        class: [2108x1 uint8]
if isfield(voc_rec, 'objects')
  gt_boxes = cat(1, voc_rec.objects(:).bbox);
  all_boxes = cat(1, gt_boxes, boxes);
  gt_classes = class_to_id.values({voc_rec.objects(:).class});
  gt_classes = cat(1, gt_classes{:});
  num_gt_boxes = size(gt_boxes, 1);
else
  gt_boxes = [];
  all_boxes = boxes;
  gt_classes = [];
  num_gt_boxes = 0;
end
num_boxes = size(boxes, 1);

rec.gt = cat(1, true(num_gt_boxes, 1), false(num_boxes, 1));
rec.overlap = zeros(num_gt_boxes+num_boxes, class_to_id.Count, 'single');
for i = 1:num_gt_boxes
  rec.overlap(:, gt_classes(i)) = ...
      max(rec.overlap(:, gt_classes(i)), boxoverlap(all_boxes, gt_boxes(i, :)));
end
rec.boxes = single(all_boxes);
rec.feat = [];
rec.class = uint8(cat(1, gt_classes, zeros(num_boxes, 1)));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function nfgCompareProjectomeToGold(phantomDir, projType, bWriteBundledProj, bDoTrueError, bDoArclength)
%Compare projectome with gold standard.
%
%   nfgCompareProjectomeToGold(phantomDir, projType, bWriteBundledProj)
%
% This loads the gold standard pathways and bundle info.  For each bundle,
% the corresponding pathways are found in the projectome solution and the
% volume of the gold standard and the projectome estimate bundles are
% calculated.
%
% Input the projectome type is a cell array and can be any or all of
%   ('stt', 'stt_bm', 'ctr_bm', 'trk', or 'gold').
%
% NOTES:

if ieNotDefined('bWriteBundledProj'); bWriteBundledProj=0; end
if isempty(projType) || ~iscellstr(projType)
    error('Must provide cell array of projectome types to compare!');
end

% Directories
fiberDir = nfgGetName('fiberDir',phantomDir);
% Input Files
volExFile = nfgGetName('volExFile',phantomDir);
goldPDBFile = nfgGetName('goldPDBFile',phantomDir);
goldInfoFile = nfgGetName('goldInfoFile',phantomDir);
projFile={};
for pp=1:length(projType)
    projFile{pp} = nfgGetName([projType{pp},'PDBFile'],phantomDir); %#ok<AGROW>
end
% Output
bundleOutFile = fullfile(fiberDir,'gold');

% Validate existence of directories
if ~isdir(fiberDir)
    error('Must provide a valid path to an NFG simulation phantom!');
end

% What types of comparisons should we do?
if ieNotDefined('bDoTrueError'); bDoTrueError = 0; end
if ieNotDefined('bDoArclength'); bDoArclength = 1; end
if ieNotDefined('bDoFiberCount'); bDoFiberCount = 0; end

% Get bundle ID and radius of each gold pathway from info file
gi = load(goldInfoFile);
g_bundleID = zeros(1,length(gi.strand_info));
g_radius = zeros(1,length(gi.strand_info));
for ii=1:length(gi.strand_info)
    g_bundleID(ii) = gi.strand_info(ii).bundleID+1;
    g_radius(ii) = gi.strand_info(ii).radius;
end

% XXX For testing
%g_bundleID(g_bundleID>30) = -1;

% Set the distance threshold to be equal to the diameter of the gold fiber
% bundle + the voxel resolution.  Ideally we would like the distance to be the radius of the gold
% plus the radius of the projectome, but we don't know the radii for all
% projectome estimates
vol = niftiRead(volExFile);
dThresh = max(vol.pixdim) + 2*max(g_radius);

% Compare Gold to the Projectome
disp('Importing gold projectome ...');
fgG = mtrImportFibers(goldPDBFile);
% Prepare output file
fg = dtiNewFiberGroup();
fg.fibers = fgG.fibers;
fgAllNonSel = dtiNewFiberGroup();
iAllNonSelGroup = [];
iGroupStat = zeros(1,length(fgG.fibers));
iBundleStat = g_bundleID;
% Radii to search through to compare projectomes
vR = 0.005:0.0005:0.07;
% For common summary info
fpCount = zeros(1,length(projFile));
tpCount = zeros(1,length(projFile));
% For arclength calculations
alE = zeros(length(projFile),max(g_bundleID));
fvol = zeros(length(projFile),max(g_bundleID));
fvolG = zeros(1,max(g_bundleID));
% For trueerror calculations
teE = zeros(length(projFile),max(g_bundleID));
teW = zeros(1,max(g_bundleID));
% For normalized fiber count comparison
%nfc = zeros(length(projFile),max(g_bundleID));
%nfcG = nfgGetNormalizeFiberCount(g_bundleID,g_bundleID);
numBundlesGold = sum(unique(g_bundleID)>0);
numFibersGold = sum(g_bundleID>0);
vID = {};
for pp=1:length(projFile)
    % Get the projectome to test
    disp(' ');disp(['Importing test projectome ' projFile{pp} ' ...']);
    fgProj = mtrImportFibers(projFile{pp});
    % Match the projectome fibers to bundles in the gold set
    [fibersSel, fibersNonSel, iBundleSel] = nfgMatchProjectome(fgG,g_bundleID,fgProj,dThresh);
    % False positives
    fpCount(pp) = length(fibersNonSel) / (length(fibersNonSel) + length(fibersSel)) * 100;
    % True positives
    tpCount(pp) = sum(unique(iBundleSel)>0) / numBundlesGold * 100;

    % Calculate the volume estimate of the test fibers compared with gold
    if bDoTrueError
        [teE(pp,:), teW] = nfgCompareAllVolumeTrueError(fgG,g_bundleID,g_radius,fibersSel,iBundleSel,phantomDir);
    end
    if bDoArclength
        %[alE(pp,:)] = nfgCompareAllVolumeArclength(fgG,g_bundleID,g_radius,vR,fibersSel,fibersNonSel,iBundleSel);
        bApplyHagmannW = 0;
        if strcmp(projType{pp},'trk')
            [fvol(pp,:), fvolG, fvolNon] = nfgCompareAllVolumeArclength(fgG,g_bundleID,g_radius,vR,fibersSel,fibersNonSel,iBundleSel,bApplyHagmannW);
        else
            [fvol(pp,:), fvolG, fvolNon] = nfgCompareAllVolumeArclength(fgG,g_bundleID,g_radius,vR,fibersSel,fibersNonSel,iBundleSel,0);
        end
    end
    if bDoFiberCount
        vID{pp} = iBundleSel; %#ok<AGROW>
        %nfc(pp,:) = nfgGetNormalizeFiberCount(iBundleSel,g_bundleID);
    end

    % Preparing for output pdb
    nSel = length(fibersSel);
    nNonSel = length(fibersNonSel);
    fg.fibers(end+1:end+nSel) = fibersSel;
    fg.fibers(end+1:end+nNonSel) = fibersNonSel;
    iGroupStat(end+1:end+nSel) = pp;
    iGroupStat(end+1:end+nNonSel) = pp;
    %fgAllNonSel.fibers(end+1:end+nNonSel) = fibersNonSel;
    %iAllNonSelGroup(end+1:end_nSel) = pp;
    iBundleStat(end+1:end+nSel) = iBundleSel;
    iBundleStat(end+1:end+nNonSel) = -1;
    [foo, projName] = fileparts(projFile{pp});
    bundleOutFile = [bundleOutFile '-' projName]; %#ok<AGROW>
end

[nfc, nfcG] = nfgGetNormalizeFiberCount(vID,g_bundleID);

% Plotting and Reporting
summary_info = {};

if bDoTrueError
    % First plot the unweighted trueerror volume error
    figure;
    bar(teE','group');
    xlabel('Bundle ID');
    ylabel('Within Bundle Volume Error (% compared to bundle volume)');
    title('Unnormalized trueError Estimation');
    legend(projType);

    % Second plot trueerror volume error normalized by total volume
    figure;
    teNE = teE.*repmat(teW/sum(teW),size(teE,1),1);
    bar(teNE','group');
    xlabel('Bundle ID');
    ylabel('Within Bundle Volume Error (% compared to total volume)');
    title('Normalized trueError Estimation');
    legend(projType);
end

if bDoFiberCount
    % Do plots by sorting NFC Gold
    [foo, iSort] = sort(nfcG,'descend');
    % First plot the normalized fiber counts
    figure;
    bar([nfcG(iSort)*numFibersGold; nfc(:,iSort)*numFibersGold]','group');
    xlabel('Bundle ID');
    ylabel('Number of fibers compared to smallest');
    title('Normalized Fiber Count');
    legend(['Gold', projType(:)']);

    % Second plot the normalized fiber count error
    figure;
    nfcE = abs(repmat(nfcG*numFibersGold,size(nfc,1),1)-nfc*numFibersGold);
    bar(nfcE(:,iSort)','group');
    xlabel('Bundle ID');
    ylabel('Difference betweeen NFC and NFC Gold');
    title('Normalized Fiber Count Error');
    legend(projType);
    strE = 'NFC Error divided by number of bundles: ';
    for pp=1:length(projType)
        strE = [strE num2str(sum(nfcE(pp,:))/numBundlesGold) ' (' projType{pp} ') ']; %#ok<AGROW>
    end
    summary_info{end+1} = strE;
end

if bDoArclength
    % Do plots by sorting NFC Gold
    [foo, iSort] = sort(fvolG,'descend');
    % First plot the normalized fiber counts
    figure;
    bar([fvolG(iSort)*100; fvol(:,iSort)*100]','group');
    xlabel('Bundle ID');
    ylabel('Bundle Volume Percent of Total');
    title('Normalized Fiber Volume');
    legend(['Gold', projType(:)']);

    % Second plot the normalized fiber count error
%     figure;
     fvolE = abs(repmat(fvolG,size(fvol,1),1)-fvol)*100;
%     bar(nfcE(:,iSort)','group');
%     xlabel('Bundle ID');
%     ylabel('Difference betweeen NFC and NFC Gold');
%     title('Normalized Fiber Count Error');
%     legend(projType);
    strE = 'Max/Sum Volume Error Percentage: ';
    for pp=1:length(projType)
        strE = [strE num2str(max(fvolE(pp,:))) '/' num2str(sum(fvolE(pp,:))) ' (' projType{pp} ') ']; %#ok<AGROW>
    end
    summary_info{end+1} = strE;
end

% False positive report
strFPs = 'False Positives: ';
for pp=1:length(projFile)
    strFPs = [strFPs num2str(fpCount(pp)) '% (' projType{pp} ') ']; %#ok<AGROW>
end
summary_info{end+1} = strFPs;
% True positive report
strTPs = 'True Positives: ';
for pp=1:length(projFile)
    strTPs = [strTPs num2str(tpCount(pp)) '% (' projType{pp} ') ']; %#ok<AGROW>
end
summary_info{end+1} = strTPs;

% Write out full report
disp(' '); disp('Summary Info');
for ii=1:length(summary_info)
    disp(summary_info{ii});
end

% Write out pathway file that combines bundled gold and bundled projectome
if bWriteBundledProj
    fg = dtiClearQuenchStats(fg);
    fg = dtiCreateQuenchStats(fg,'Bundle ID',iBundleStat);
    fg = dtiCreateQuenchStats(fg,'Group',iGroupStat);
    % Need to put the extension on the file
    bundleOutFile = [bundleOutFile '.pdb'];
    mtrExportFibers(fg,bundleOutFile);
end

return;

% Color wheel
% vC = {'b','g','r','c','m','y','k'};

%
%     % Calculate normalized true error volume error and store mean
%     teNE = teE.*teW/sum(teW);
%     pE1W(pp) = mean(teNE(teE<100 & teE>0));
%     % Store unnormalized volume error as well
%     pE1(pp) = mean(teE(teE<100 & teE>0));
%
%     % Plotting and Reporting Error
%     percentE = sum(abs(E),2)/sum(volG)*100;
%     %figure(100);
%     %plot(vR,percentE,vC{pp}); hold on;
%     figure(100);
%     teNE = teE.*teW/sum(teW);
%     plot(teNE,vC{pp}); hold on;
%     pE1W(pp) = mean(teNE(teE<100 & teE>0));
%     figure(101);
%     plot(teE,vC{pp}); hold on;
%     pE1(pp) = mean(teE(teE<100 & teE>0));
%     [minE, minI] = min(percentE);
%     disp(['Minimum within bundle error of ' num2str(minE) '% found with radius ' num2str(vR(minI))]);
%     percentENonBundle = volPNonBundle * vR(minI)^2 / sum(volG) * 100;
%
%     %disp([projName ' non-bundle error is ' num2str(percentENonBundle) '%']);
%     pE2(pp) = 100*length(fibersNonSel)/(length(fibersNonSel)+length(fibersSel));
%     %disp([projName ' non-bundle error is ' num2str(100*length(fibersNonSel)/(length(fibersNonSel)+length(fibersSel))) '%.']);
%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function IncEditDataType(topH, inc)

% function IncEditDataType(topH, inc);
%
% Update the dataType structure with the present information, then
% increment the scan index by the specified amount [inc] from its
% present value. Load the scan content fields that correspond to
% the new index. Input topH is the handle id of the top-level
% edit figure.
%
% DBR 4/99

UpdateEditDataType(topH);

% Increment the present contents of the scan-index field:
uiData = get(topH, 'UserData');
iScan = inc + round(str2num(get(uiData.hScan, 'string')));
nScans = length(uiData.dataType.scanParams);
if iScan > nScans
    iScan = nScans;
end
if iScan < 1
    iScan = 1;
end
set(uiData.hScan, 'string', int2str(iScan));
uiData.iScan = iScan;
set(topH, 'UserData', uiData);

% Load the fields based on the new scan index and
% the present contents of the dataType structure.
nFields = length(uiData.scanData);
for iField=1:nFields
    fieldName = uiData.scanData(iField).field;
    val = uiData.dataType.scanParams(iScan).(fieldName);
    set(uiData.scanData(iField).handle, 'string', num2str(val));
end

nFields = length(uiData.blockData);
for iField=1:nFields
    fieldName = uiData.blockData(iField).field;
    val = uiData.dataType.blockedAnalysisParams(iScan).(fieldName);
    set(uiData.blockData(iField).handle, 'string', num2str(val));
end

nFields = length(uiData.eventData);
for iField=1:nFields
    fieldName = uiData.eventData(iField).field;
    val = uiData.dataType.eventAnalysisParams(iScan).(fieldName);
    set(uiData.eventData(iField).handle, 'string', num2str(val));
end

set(topH, 'UserData', uiData);

return


text:
function h = supportFunction(polygon, varargin)
%SUPPORTFUNCTION Compute support function of a polygon.
%
%   H = supportFunction(POLYGON, N)
%   uses N points for suport function approximation
%
%   H = supportFunction(POLYGON)
%   assume 24 points for approximation
%
%   H = supportFunction(POLYGON, V)
%   where V is a vector, uses vector V of angles to compute support
%   function.
%
%   See also
%   polygons2d, convexification

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-12-20
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

N = 24;
u = 0:2*pi/N:2*pi*(1-1/N);

if length(varargin)==1
    var = varargin{1};
    if length(var)==1
        N = var;
        u = 0:2*pi/N:2*pi*(1-1/N);
    else
        u = var;
    end
end

% ensure u vertical vector
if size(u, 1)==1
    u=u';
end

h = zeros(size(u));

for i=1:length(u)

    v = repmat([cos(u(i)) sin(u(i))], [size(polygon, 1), 1]);

    h(i) = max(dot(polygon, v, 2));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [newstr,same]=samenum(same,num) %#codegen
coder.extrinsic('strrep','strcat')
newstr=[];

for n=1:1:(length(same)-num+1)
    if n>length(same)-num+1
        break
    end
num3=find(same(n:1:n+num-1)==same(n));
if length(num3)==num
   same= strrep(same,same(n:1:n+num-1),same(n));
   newstr=strcat(newstr,same(n));
   sound(sin(1:2:10))

else

end
end

end

author: Sable

repo: mcbench-benchmarks


text:
% coordinates of the three rings (blue, red and green)
[X1,Y1,Z1]=torus(5-0.3,30,0.3,2); % the top (blue)
[X2,Z2,Y2]=torus(5+0.3,30,0.3,2); % inner gimbal (green, semi-transparent)
[X3,Z3,Y3]=torus(5+0.9,30,0.3,1); % outer gimbal (red, semi-transparent)

% handles to surf plots of the three rings
H1=surf(X1,Y1,Z1,'EdgeColor','none','FaceColor','blue'); hold on;
H2=surf(X2,Y2,Z2,'EdgeColor','none','FaceColor','green','FaceAlpha',0.3);
H3=surf(X3,Y3,Z3,'EdgeColor','none','FaceColor','red','FaceAlpha',0.3);

% creates the bars inside the inner (blue) ring. The wide cylinders
% represent weights.
[Xbar1,Ybar1,Zbar1,Xbar2,Ybar2,Zbar2]=createInnerBars(0.2,0.5,0.2,0.5);
bar1=surf(Xbar1,Ybar1,Zbar1,'EdgeColor','none','FaceColor','blue');
bar2=surf(Xbar2,Ybar2,Zbar2,'EdgeColor','none','FaceColor','blue');

% the red motionless supporting bar
[ Xupper,Yupper,Zupper,Xlower,Ylower,Zlower ] = createOuterBars( 0.3 );
surf(Xlower,Ylower,Zlower,...
    'EdgeColor','none','FaceColor','red','FaceAlpha',0.3);

hold off;
 camlight left;
 lighting gouraud; % phong is more demanding but gives nicer results
 bb=7;
 axis([-bb bb -bb bb -bb bb]);
axis square;
xlabel('x'); ylabel('y'); zlabel('z');


text:
function varargout = spm_brainwarp(varargin)
% Part of nonlinear spatial normalisation - a compiled routine
%_______________________________________________________________________
% [Alpha,Beta,Var] = spm_brainwarp(VG,VF,Affine,basX,basY,basZ,...
%                                   dbasX,dbasY,dbasZ,T,fwhm,VW,VW2)
% VG    - Template volume(s) (see spm_vol)
% VF    - Object volume
% Affine    - The affine transformation which maps between the object
%     and template.
% basX  - Basis vectors in X. # rows must eq. VG(1)
% basY  - Basis vectors in Y. # rows must eq. VG(2)
% basZ  - Basis vectors in Z. # rows must eq. VG(3)
% dbasX - Derivatives of basis vectors in X. # rows must eq. VG(1)
% dbasY - Derivatives of basis vectors in Y. # rows must eq. VG(2)
% dbasZ - Derivatives of basis vectors in Z. # rows must eq. VG(3)
% T - The current parameter estimates.
% fwhm  - The approximate smoothness of the images.
% VW    - an optional weighting volume for determining which voxels
%         should be weighted more heavily in the fitting process.
%         This volume should have the same dimensions and position
%         as the volumes in VG.
% VW2   - another optional weighting volume for determining which voxels
%         should be weighted more heavily in the fitting process.
%         This volume should have the same dimensions and position
%         as the volumes in VF.
% Without the weighting volumes, all voxels are assigned weights that
% are uniformly one.
%
% Alpha - A*A - where A is the design matrix
% Beta  - A*b - where f is the object image
% Var   - the approximate chi^2 (corrected for number of resels).
%_______________________________________________________________________
%
% The voxels of g1, g2.. are sampled according to the smoothness of the
% image (fwhm). The corresponding voxels of f are determined according
% to the current parameter estimates and the affine transform.  See
% "spm_write_sn.m" for more details about how this is done.
%
%
%-----------------------------------------------------------------------
%
% The design matrix A is generated internally from:
%
% diag(w)*[diag(df/dx)*B diag(df/dy)*B diag(df/dz)*B ...
%      diag(g1)*[1 x y z] ...
%      diag(g2)*[1 x y z] ...]
%
% where df/dx, df/dy & df/dz are column vectors containing the gradient
%   of image f with respect to displacements in x, y & z
%   (in the space of g).
%
%   B is generated from kron(basZ,kron(basY,BasX)). Each column of
%   B is a basis image.
%
%   g1, g2.. are template images.
%
%   x, y & z are simply the spatial coordinates of the voxels of f.
%
%   s1, s2.. are the current estimates for the required scaling
%   factors. These are derived from T(3*prod(VG(1:3))+1),
%   T(3*prod(VG(1:3))+2)...
%
%       w is an optional vector of weights, where w = 1/(1/w1 + 1/w2)
%       where w1 and w2 are derived from the optional weighting images.
%
% The vector b contains [diag(w)*(f - diag(g1)*s1 - diag(g1)*x*s2 - ...)].
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

%-This is merely the help file for the compiled routine
error('spm_brainwarp.c not compiled - see Makefile')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
load hall1-200;
matlabpool open;
X = normalize(XO);
[P Q L] = onlineRPMF(X, 2, 1, 1, 1e-2, ones(size(X)));
% M = normalize(I0);
% [P2 Q2 L2] = onlineRPMF(M, 2, 1, 1, 1e-2, ones(size(M)));
%show(X, L, abs(X - L), [144 176]);
show(X, L, (X - L), [144 176]);
%show(M, L2, (M - L2), [vidHeight vidWidth]);
matlabpool close;

author: andrewssobral

repo: lrslibrary


text:
classdef ProjectorsTestsSuite < handle

    methods

        function obj = ProjectorsTestsSuite()
            warning('off', 'MATLAB:structOnObject')
            results = runtests("ProjectorsTests","ProcedureName","testProjectorsInFEM", 'Debug', true);
            table(results)
            warning('on', 'MATLAB:structOnObject')
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [blockInfo] = PL2ReadNextDataBlock(blockInfo)
% PL2ReadNextDataBlock(blockInfo): read next data block, return block info
%
% example:
%
% blockInfo = PL2ReadFirstDataBlock(filename)
% while blockInfo.Type > 0
%    % process data from the block
%    % ...
%    % read the next data block
%    blockInfo = internalPL2ReadNextDataBlock(blockInfo);
% end
%
%
% INPUT:
%   blockInfo - a structure returned from PL2ReadFirstDataBlock or PL2ReadNextDataBlock
%
% OUTPUT:
%   blockInfo - information about the data block
%
%   blockInfo.Type - type of data block
%                    0 - not a data block, did not find data in the file
%                    1 - spike data type
%                    2 - analog data type
%                    3 - digital event data type
%                    4 - start/stop event data type
%   blockInfo.Source - block source
%   blockInfo.Channel - block channel inside source
%   if blockInfo.Type is 1 (spike), SpikeData structure is filled
%       blockInfo.SpikeData.WaveLength - number of samples in each waveform
%       blockInfo.SpikeData.NumSpikes - number of spikes
%       blockInfo.SpikeData.Timestamps  - spike timestamps in seconds
%       blockInfo.SpikeData.Units - spike units
%       blockInfo.SpikeData.Waveforms - spike waveforms in milliVolts
%   if blockInfo.Type is 2 (analog), AnalogData structure is filled
%       blockInfo.AnalogData.NumSamples - number of analog samples
%       blockInfo.AnalogData.Timestamp - timestamps of the first sample (in seconds)
%       blockInfo.AnalogData.Values - sample values in milliVolts
%   if blockInfo.Type is 3 (event), EventData structure is filled
%       blockInfo.EventData.NumEvents - number of events
%       blockInfo.EventData.Timestamps - event timestamps in seconds
%       blockInfo.EventData.Values - event values
%   if blockInfo.Type is 4 (start/stop event), StartStopEventData structure is filled
%       blockInfo.StartStopEventData.NumEvents - number of events
%       blockInfo.StartStopEventData.Timestamps - event timestamps in seconds
%       blockInfo.StartStopEventData.Values - event values

blockInfo = internalPL2ReadNExtDataBlock(blockInfo);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%   This subroutine selcts the Ni closest earthquakes
%   around a userdefined point
%

report_this_filefun(mfilename('fullpath'));

new = a;
ni  = input('Please input  number of events ni:')
%ni = 100
axes(h6)
[xa0,ya0]  = ginput(1);

l = sqrt(((a.Longitude-xa0)*cos(pi/180*ya0)*111).^2 + ((a.Latitude-ya0)*111).^2) ;
[s,is] = sort(l);
new = a(is(:,1),:) ;
plos1 = plot(new(1:ni,1),new(1:ni,2),'xw','EraseMode','back');
%plos1 = plot(new(1:ni,1),new(1:ni,2),'xw')

figure_w_normalized_uicontrolunits(2)
clf
h3 = gcf;
newt = new(1:ni,:);
[st,ist] = sort(newt);
newt2 = newt(ist(:,3),:);

newt2(:,9) = newt2.Date + newt2.Date.Month/12 + newt2.Date.Day/365;
[st,ist] = sort(newt2);
newt3 = newt2(ist(:,9),:);

%figure_w_normalized_uicontrolunits(3)
%clf
%plode = plot(newt3(:,9),-newt3(:,8),'o')
%axis([81 92.6 -20 0])
%grid
%xlabel('Time')
%ylabel('Depth in km')

timeplot

author: CelsoReyes


text:
function displayInterRes( imgIn )
%DISPLAYINTERRES display intermediate results
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% set all parameters here inside

persistent paraIR;

if(ischar(imgIn))
    clear paraIR;
    return;
end

if(~isfield(paraIR,'dim'))
    paraIR.dim = ndims(imgIn);
    if(paraIR.dim == 3)
        % y-x-cha
        paraIR.subs = {':',':',':'};
        paraIR.permRule = [];

    elseif(paraIR.dim == 4)
        % 3D: y-z-x-cha

        paraIR.subs = {':', ':', ceil(size(imgIn,2)/2), ':'};
        paraIR.permRule = [1 3 2 4]; % => y-x-z-cha

%         % OR
%         % 2Dt: t-y-x-cha
%         paraIR.subs = {':', ':', ceil(size(imgIn,1)/2), ':'};
%         paraIR.permRule = [2 3 1 4]; % => y-x-t-cha

    elseif(paraIR.dim == 5)
        % t-y-z-x-cha
        paraIR.subs = {':', ':', ceil(size(imgIn,3)/2), ceil(size(imgIn,1)/2), ':'};
        paraIR.permRule = [2 4 3 1 5]; % => y-x-z-t-cha
    end

    paraIR.range = [0 1];
    paraIR.iter = 1;
    paraIR.normF = [];
    paraIR.screensize = get(0,'ScreenSize');

    paraIR.plotNorm = false;
end

if(~isempty(paraIR.permRule))
    imgIn = permute(imgIn,paraIR.permRule);
end
imgIn = imgIn(paraIR.subs{:});
imgIn = ((imgIn - min(imgIn(:))) * (paraIR.range(2)-paraIR.range(1)))./(max(imgIn(:)) - min(imgIn(:)));

% show channel combined abs image
imgIn = sqrt(sum(abs(imgIn).^2, ndims(imgIn)));
paraIR.normF(end+1) = norm(imgIn(:));

paraIR.hfig = figure(999);
set(paraIR.hfig,'position',[paraIR.screensize(1)+5 paraIR.screensize(2)+50 paraIR.screensize(3)-10 paraIR.screensize(4)-130]);
if(paraIR.plotNorm)
    subplot(2,1,1);
end
imagesc( imgIn, paraIR.range );
title(sprintf('\\bfiteration %g \nFrobenius norm: %g', paraIR.iter, norm(imgIn(:))))
axis image
colormap(gray(256))
% drawnow
if(paraIR.plotNorm)
    subplot(2,1,2);
    plot(paraIR.normF);
    xlabel('Iteration');
    ylabel('Frobenius norm');
    title('Frobenius norm behaviour');
end
drawnow
pause(0.5);

paraIR.iter = paraIR.iter + 1;

end

author: thomaskuestner


text:
function Hd = gen_iir_coeffs
%GEN_IIR_COEFFS Returns a discrete-time filter object.

%
% M-File generated by MATLAB(R) 7.7 and the Signal Processing Toolbox 6.10.
%
% Generated on: 16-Dec-2008 10:06:34
%

% Elliptic Lowpass filter designed using FDESIGN.LOWPASS.

% All frequency values are in Hz.
Fs = 1024;  % Sampling Frequency

Fpass = 100;     % Passband Frequency
Fstop = 125;     % Stopband Frequency
Apass = 1;       % Passband Ripple (dB)
Astop = 80;      % Stopband Attenuation (dB)
match = 'both';  % Band to match exactly

% Construct an FDESIGN object and call its ELLIP method.
h  = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);
set(Hd, 'Arithmetic', 'double');

% [EOF]


text:
function [D0,phisq0,st] = lme_fit_FSinit(X,Zcols,y,ni,e)
% [phisq0,D0,st] = lme_fit_FSinit(X,Zcols,y,ni,e)
%
% Starting values for linear mixed-effects estimation.This function is
% intended to be used to provide starting values for the lme_fit_NR
% function by mean of some initial iterations of the Fisher scoring algorithm.
%
% Input
% X: Ordered design Matrix (according to time for each subject).
% Zcols: Vector with the indices of the colums of X that will be considered
% as random effects.
% y: Ordered data vector (according to X).
% ni: Vector whose entries are the number of repeated measures for each
% subject (ordered according to X).
% e: Convergence epsilon. Default 10^-3;
%
% Output
% phisq0: Estimated within-subject variability.
% D0: Estimated random effects covariance matrix.
% st: Termination state (1 for convergence and 0 otherwise).
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:57 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:57 $
%    $Revision: 1.2 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%
if nargin < 4
    error('Too few inputs');
elseif nargin < 5
    e = 10^-3;
end;
nit = 50;
st = 1;
m = length(ni);
p = size(X,2);
q = length(Zcols);
ind = [false(q*q,1);true];
for k=1:q
    for j=1:k
        ind((k-1)*q+j) = true;
    end;
end;
n = sum(ni);
Z = X(:,Zcols);
W = zeros(n,max(ni));
SIGMA = W;

%Starting values
[D,phisq] = lme_fit_init(X,Zcols,y,ni);
L = chol(D);
phi = sqrt(phisq);
theta = [vec(L);phi];

%% Iterations
tf = true;
it = 0;
display('Starting initial Fisher scoring iterations');
while tf
    it = it+1;
    %Computation of W = SIGMA^-1 and H.
    posi = 1; H = 0; Term = 0;
    scInvD = D\eye(q)*phisq;
    for i=1:m
        posf = posi+ni(i)-1;
        Zi = Z(posi:posf,:);
        Wi = (eye(ni(i))-Zi/(Zi'*Zi+scInvD)*Zi')/phisq;
        W(posi:posf,1:ni(i)) = Wi;
        SIGMA(posi:posf,1:ni(i)) = Zi*D*Zi'+ eye(ni(i))*phisq;
        Xi = X(posi:posf,:);
        Ti = Xi'*Wi;
        H = H + Ti*Xi;
        Term = Term + Ti*y(posi:posf);
        posi = posf+1;
    end;
    invH = H\eye(p);
   %Estimation
    Bhat = invH*Term;
    r = y-X*Bhat;
    posi = 1; lreml = 0;
    for i=1:m
        posf = posi+ni(i)-1;
        Wi = W(posi:posf,1:ni(i));
        ri = r(posi:posf);
        lreml = lreml + log(det(Wi))-ri'*Wi*ri;
        posi = posf+1;
    end;
    gr = lme_Gradient(X,Zcols,W,invH,L,phi,r,ni);
    EI = lme_EI(X,Zcols,W,invH,SIGMA,L,phi,ni);
    theta(ind) = theta(ind) + EI\gr;
    L = reshape(theta(1:end-1),q,q);
    D = L'*L;
    phi = theta(end);
    phisq = phi*phi;

    %Restricted log-likelihood
    lreml = 0.5*(lreml - log(det(H)));
    display(['Likelihood at FS iteration ' num2str(it) ' : ' num2str(lreml)]);
    eps = norm(gr);
    display(['Gradient norm: ' num2str(eps)]);

    %Termination
    if (it==nit) || (eps<e)
        tf = false;
        phisq0 = phisq;
        D0 = D;
        if it == nit
            st = 0;
            display(['Initial FS does not converge after ' num2str(nit)...
                                                        ' iterations!!!']);
        end;
    end;
end


text:
classdef FourthOrder3DTensor < AbstractTensor ...
                                        & FourthOrderDescriptor ...
                                        & TensorRepresentation ...
                                        & Elasticity3dDescriptor

    methods (Access = public)

        function obj = FourthOrder3DTensor()
        end

        function createRandomTensor(obj)
            obj.createRandomTensor@AbstractTensor();
        end

    end

    methods (Access = protected)

        function loadTensorSize(obj)
            obj.tensorSize = [3,3,3,3];
        end
    end

end


text:
function candoit = canGetHessian(problem)
% Checks whether the Hessian can be computed for a problem structure.
%
% function candoit = canGetHessian(problem)
%
% Returns true if the Hessian of the cost function can be computed given
% the problem description, false otherwise.
%
% See also: canGetCost canGetDirectionalDerivative canGetGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    candoit = isfield(problem, 'hess') || ...
              (isfield(problem, 'ehess') && canGetEuclideanGradient(problem));

    % Display an extra warning message to the user in anticipation of
    % common mistakes.
    if ~candoit && ...
           (isfield(problem, 'ehess') && ~canGetEuclideanGradient(problem))
        warning('manopt:canGetHessian', ...
               ['If the Hessian is supplied as a Euclidean Hessian (ehess),\n' ...
                'then the Euclidean gradient must also be supplied (egrad).']);
    end

end

author: NicolasBoumal

repo: manopt


text:
classdef DescriptorMatcher < handle
    %DESCRIPTORMATCHER  Common interface for matching keypoint descriptors
    %
    % Class for matching keypoint descriptors.
    %
    % Matchers of keypoint descriptors in OpenCV have wrappers with a common
    % interface that enables you to easily switch between different algorithms
    % solving the same problem. This section is devoted to matching
    % descriptors that are represented as vectors in a multidimensional space.
    % All objects that implement vector descriptor matchers inherit the
    % DescriptorMatcher interface.
    %
    % It has two groups of match methods: for matching descriptors of an image
    % with another image or with an image set.
    %
    % ## Example
    %
    %     X = rand(100,10);
    %     Y = rand(100,10);
    %     matcher = cv.DescriptorMatcher('BruteForce');
    %     matcher.add(X);
    %     matcher.train();  % Optional for BruteForce matcher
    %     matches = matcher.match(Y);
    %
    % See also: cv.DescriptorExtractor, cv.FeatureDetector, cv.drawMatches,
    %  matchFeatures
    %

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

    methods
        function this = DescriptorMatcher(matcherType, varargin)
            %DESCRIPTORMATCHER  Creates a descriptor matcher by name
            %
            %     matcher = cv.DescriptorMatcher(type)
            %     matcher = cv.DescriptorMatcher(type, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __type__ In the first variant, it creates a descriptor matcher
            %   of a given type with the default parameters (using default
            %   constructor). The following types are recognized:
            %
            %   * __BruteForce__ (default) L2 distance
            %   * __BruteForce-SL2__ L2SQR distance
            %   * __BruteForce-L1__ L1 distance
            %   * __BruteForce-Hamming__, __BruteForce-HammingLUT__
            %   * __BruteForce-Hamming(2)__
            %   * __FlannBased__ Flann-based indexing
            %
            %   In the second variant, it creates a matcher of the given type
            %   using the specified parameters. The following descriptor
            %   matcher types are supported:
            %
            %   * __BFMatcher__ Brute-force descriptor matcher. For each
            %     descriptor in the first set, this matcher finds the closest
            %     descriptor in the second set by trying each one. This
            %     descriptor matcher supports masking permissible matches of
            %     descriptor sets.
            %   * __FlannBasedMatcher__ Flann-based descriptor matcher. This
            %     matcher trains `flann::Index_` on a train descriptor
            %     collection and calls its nearest search methods to find the
            %     best matches. So, this matcher may be faster when matching a
            %     large train collection than the brute force matcher.
            %     `FlannBasedMatcher` does not support masking permissible
            %     matches of descriptor sets because `flann::Index` does not
            %     support this.
            %
            % ## Options
            % The Brute-force matcher constructor (`BFMatcher`) accepts the
            % following options:
            %
            % * __NormType__ See cv.DescriptorExtractor.defaultNorm, default
            %   'L2'. One of:
            %   * __L1__, __L2__ L1 and L2 norms are preferable choices for
            %     cv.SIFT and cv.SURF descriptors.
            %   * __Hamming__ should be used with cv.ORB, cv.BRISK and
            %     cv.BriefDescriptorExtractor.
            %   * __Hamming2__ should be used with cv.ORB when `WTA_K` equals
            %     3 or 4 (see cv.ORB.WTA_K description).
            % * __CrossCheck__ If it is false, this is will be default
            %   `BFMatcher` behaviour when it finds the `k` nearest neighbors
            %   for each query descriptor. If `CrossCheck==true`, then the
            %   cv.DescriptorMatcher.knnMatch method with `k=1` will only
            %   return pairs `(i,j)` such that for i-th query descriptor the
            %   j-th descriptor in the matcher's collection is the nearest and
            %   vice versa, i.e. the `BFMatcher` will only return consistent
            %   pairs. Such technique usually  produces best results with
            %   minimal number of outliers when there are enough matches. This
            %   is alternative to the ratio test, used by *D. Lowe* in SIFT
            %   paper. default false
            %
            % The Flann-based matcher constructor (`FlannBasedMatcher`) takes
            % the following optional arguments:
            %
            % * __Index__ Type of indexer, default 'KDTree'. One of the below.
            %   Each index type takes optional arguments (see IndexParams
            %   options below). You can specify the indexer by a cell array
            %   that starts from the type name followed by option arguments:
            %   `{'Type', 'OptionName',optionValue, ...}`.
            %   * __Linear__ Brute-force matching, linear search
            %   * __KDTree__ Randomized kd-trees, parallel search
            %   * __KMeans__ Hierarchical k-means tree
            %   * __HierarchicalClustering__ Hierarchical index
            %   * __Composite__ Combination of KDTree and KMeans
            %   * __LSH__ Multi-probe LSH
            %   * __Autotuned__ Automatic tuning to one of the above
            %     (`Linear`, `KDTree`, `KMeans`)
            %   * __Saved__ Load saved index from a file
            %
            % * __Search__ Option in matching operation. Takes a cell array of
            %   option pairs:
            %   * __Checks__ The number of times the tree(s) in the index
            %     should be recursively traversed. A higher value for this
            %     parameter would give better search precision, but also take
            %     more time. If automatic configuration was used when the
            %     index was created, the number of checks required to achieve
            %     the specified precision was also computed, in which case
            %     this parameter is ignored. -1 for unlimited. default 32
            %   * __EPS__ search for eps-approximate neighbours. default 0
            %   * __Sorted__ only for radius search, require neighbours sorted
            %     by distance. default true
            %
            % ## IndexParams Options for `FlannBasedMatcher`
            %
            % The following are the options for FLANN indexers
            % (Fast Library for Approximate Nearest Neighbors):
            %
            % ### `Linear`
            % Linear index takes no options.
            %
            % ### `Saved`
            % Saved index takes only one argument specifying the filename.
            %
            % ### `KDTree` and `Composite`
            % * __Trees__ The number of parallel kd-trees to use. Good values
            %   are in the range [1..16]. default 4
            %
            % ### `KMeans` and `Composite`
            % * __Branching__ The branching factor to use for the hierarchical
            %   k-means tree. default 32
            % * __Iterations__ The maximum number of iterations to use in the
            %   k-means clustering stage when building the k-means tree. A
            %   value of -1 used here means that the k-means clustering should
            %   be iterated until convergence. default 11
            % * __CentersInit__ The algorithm to use for selecting the initial
            %   centers when performing a k-means clustering step. The
            %   possible values are (default is 'Random'):
            %   * __Random__ picks the initial cluster centers randomly
            %   * __Gonzales__ picks the initial centers using Gonzales
            %     algorithm
            %   * __KMeansPP__ picks the initial centers using the algorithm
            %     suggested in [ArthurKmeansPP2007]
            %   * __Groupwise__ chooses the initial centers in a way inspired
            %     by Gonzales (by Pierre-Emmanuel Viel).
            % * __CBIndex__ This parameter (cluster boundary index) influences
            %   the way exploration is performed in the hierarchical kmeans
            %   tree. When `CBIndex` is zero the next kmeans domain to be
            %   explored is chosen to be the one with the closest center. A
            %   value greater then zero also takes into account the size of
            %   the domain. default 0.2
            %
            % ### `HierarchicalClustering`
            % * __Branching__ same as above.
            % * __CentersInit__ same as above.
            % * __Trees__ same as above.
            % * __LeafSize__ maximum leaf size. default 100
            %
            % ### `LSH`
            % * __TableNumber__ The number of hash tables to use (usually
            %   between 10 and 30). default 20
            % * __KeySize__ The length of the key in the hash tables (usually
            %   between 10 and 20). default 15
            % * __MultiProbeLevel__ Number of levels to use in multi-probe
            %   (0 is regular LSH, 2 is recommended). default 0
            %
            % ### `Autotuned`
            % * __TargetPrecision__ Is a number between 0 and 1 specifying the
            %   percentage of the approximate nearest-neighbor searches that
            %   return the exact nearest-neighbor. Using a higher value for
            %   this parameter gives more accurate results, but the search
            %   takes longer. The optimum value usually depends on the
            %   application. default 0.8
            % * __BuildWeight__ Specifies the importance of the index build
            %   time raported to the nearest-neighbor search time. In some
            %   applications it is acceptable for the index build step to take
            %   a long time if the subsequent searches in the index can be
            %   performed very fast. In other applications it is required that
            %   the index be build as fast as possible even if that leads to
            %   slightly longer search times. default 0.01
            % * __MemoryWeight__ Is used to specify the tradeoff between time
            %   (index build time and search time) and memory used by the
            %   index. A value less than 1 gives more importance to the time
            %   spent and a value greater than 1 gives more importance to the
            %   memory usage. default 0
            % * __SampleFraction__ Is a number between 0 and 1 indicating what
            %   fraction of the dataset to use in the automatic parameter
            %   configuration algorithm. Running the algorithm on the full
            %   dataset gives the most accurate results, but for very large
            %   datasets can take longer than desired. In such case using just
            %   a fraction of the data helps speeding up this algorithm while
            %   still giving good approximations of the optimum parameters.
            %   default 0.1
            %
            % ## Example
            % For example, `KDTree` with tree size = 4 is specified by:
            %
            %     matcher = cv.DescriptorMatcher('FlannBasedMatcher', ...
            %         'Index',  {'KDTree', 'Trees', 4}, ...
            %         'Search', {'Sorted', true})
            %
            % Here is an example for loading a saved index:
            %
            %     matcher = cv.DescriptorMatcher('FlannBasedMatcher', ...
            %         'Index', {'Saved', '/path/to/saved/index.xml'})
            %
            % ## References
            % [ArthurKmeansPP2007]:
            % > Arthur and S. Vassilvitskii
            % > *"k-means++: the advantages of careful seeding"*,
            % > Proceedings of the eighteenth annual ACM-SIAM symposium
            % > on Discrete algorithms, 2007
            %
            % [LSH]:
            % > *Multi-Probe LSH: Efficient Indexing for High-Dimensional Similarity Search*
            % > by Qin Lv, William Josephson, Zhe Wang, Moses Charikar, Kai Li.,
            % > Proceedings of the 33rd International Conference on
            % > Very Large Data Bases (VLDB). Vienna, Austria. September 2007
            %
            % See also: cv.DescriptorMatcher
            %
            if nargin < 1, matcherType = 'BruteForce'; end
            this.Type = matcherType;
            this.id = DescriptorMatcher_(0, 'new', matcherType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     matcher.delete()
            %
            % See also: cv.DescriptorMatcher
            %
            if isempty(this.id), return; end
            DescriptorMatcher_(this.id, 'delete');
        end

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

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the train descriptor collection
            %
            %     matcher.clear()
            %
            % See also: cv.DescriptorMatcher.empty
            %
            DescriptorMatcher_(this.id, 'clear');
        end

        function status = empty(this)
            %EMPTY  Returns true if there are no train descriptors in the collection
            %
            %     status = matcher.empty()
            %
            % ## Output
            % * __status__ boolean status
            %
            % See also: cv.DescriptorMatcher.clear
            %
            status = DescriptorMatcher_(this.id, 'empty');
        end

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

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

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

    %% DescriptorMatcher
    methods
        function status = isMaskSupported(this)
            %ISMASKSUPPORTED  Returns true if the descriptor matcher supports masking permissible matches
            %
            %     status = matcher.isMaskSupported()
            %
            % ## Output
            % * __status__ boolean status.
            %
            % Brute-force matchers support masking, while Flann-based matchers
            % do no support masking.
            %
            % See also: cv.DescriptorMatcher.match
            %
            status = DescriptorMatcher_(this.id, 'isMaskSupported');
        end

        function descriptors = getTrainDescriptors(this)
            %GETTRAINDESCRIPTORS  Returns the train descriptor collection
            %
            %     descriptors = matcher.getTrainDescriptors()
            %
            % ## Output
            % * __descriptors__ Set of train descriptors. A cell array of
            %   matrices.
            %
            % See also: cv.DescriptorMatcher.add
            %
            descriptors = DescriptorMatcher_(this.id, 'getTrainDescriptors');
        end

        function add(this, descriptors)
            %ADD  Adds descriptors to train a descriptor collection
            %
            %     matcher.add(descriptors)
            %
            % If the collection is not empty, the new descriptors are added
            % to existing train descriptors.
            %
            % ## Input
            % * __descriptors__ Descriptors to add. Each `descriptors{i}` is a
            %   set of descriptors from the same train image. Can be either a
            %   matrix or a cell array of matrices (matrices of type `uint8`
            %   or `single`)
            %
            % See also: cv.DescriptorMatcher.getTrainDescriptors
            %
            DescriptorMatcher_(this.id, 'add', descriptors);
        end

        function train(this)
            %TRAIN  Trains a descriptor matcher
            %
            %     matcher.train()
            %
            % Trains a descriptor matcher (for example, the flann index). In all
            % methods to match, the method `train` is run every time before
            % matching. Some descriptor matchers (for example,
            % `BruteForceMatcher`) have an empty implementation of this method.
            % Other matchers really train their inner structures (for example,
            % `FlannBasedMatcher` trains `flann::Index`).
            %
            DescriptorMatcher_(this.id, 'train');
        end

        function matches = match(this, queryDescriptors, varargin)
            %MATCH  Finds the best match for each descriptor from a query set
            %
            %     matches = matcher.match(queryDescriptors, trainDescriptors)
            %     matches = matcher.match(queryDescriptors)
            %     [...] = matcher.match(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ Query set of descriptors.
            % * __trainDescriptors__ Train set of descriptors. This set is not
            %   added to the train descriptors collection stored in the class
            %   object.
            %
            % ## Output
            % * __matches__ Matches. If a query descriptor is masked out in
            %   `Mask`, no match is added for this descriptor. So,
            %   `matches` size may be smaller than the query descriptors
            %   count. A 1-by-N structure array with the following fields:
            %   * __queryIdx__ query descriptor index (zero-based index)
            %   * __trainIdx__ train descriptor index (zero-based index)
            %   * __imgIdx__ train image index (zero-based index)
            %   * __distance__ distance between descriptors (scalar)
            %
            % ## Options
            % * __Mask__ default empty
            %   * In the first form, mask specifying permissible matches
            %   between an input query and train matrices of descriptors.
            %   Matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors,1)]`.
            %   * In the second form, set of masks. Each `masks{i}` specifies
            %   permissible matches between the input query descriptors and
            %   stored train descriptors from the i-th image
            %   `trainDescCollection{i}`. Cell array of length
            %   `length(trainDescriptors)`, each a matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors{i},1)]`.
            %
            % In the first variant of this method, the train descriptors are
            % passed as an input argument. In the second variant of the
            % method, train descriptors collection that was set by
            % cv.DescriptorMatcher.add is used. Optional mask (or masks) can
            % be passed to specify which query and training descriptors can be
            % matched. Namely, `queryDescriptors(i,:)` can be matched with
            % `trainDescriptors(j,:)` only if `mask(i,j)` is non-zero.
            %
            % See also: cv.DescriptorMatcher.knnMatch,
            %  cv.DescriptorMatcher.radiusMatch
            %
            matches = DescriptorMatcher_(this.id, 'match', queryDescriptors, varargin{:});
        end

        function matches = knnMatch(this, queryDescriptors, varargin)
            %KNNMATCH  Finds the k best matches for each descriptor from a query set
            %
            %     matches = matcher.knnMatch(queryDescriptors, trainDescriptors, k)
            %     matches = matcher.knnMatch(queryDescriptors, k)
            %     [...] = matcher.knnMatch(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ Query set of descriptors.
            % * __trainDescriptors__ Train set of descriptors. This set is not
            %   added to the train descriptors collection stored in the class
            %   object.
            % * __k__ Count of best matches found per each query descriptor or
            %   less if a query descriptor has less than `k` possible matches
            %   in total.
            %
            % ## Output
            % * __matches__ Matches. Each `matches{i}` is `k` or less matches
            %   for the same query descriptor. A cell array of length
            %   `size(queryDescriptors,1)`, each cell is a 1-by-(k or less)
            %   structure array that has the following fields:
            %   * __queryIdx__ query descriptor index (zero-based index)
            %   * __trainIdx__ train descriptor index (zero-based index)
            %   * __imgIdx__ train image index (zero-based index)
            %   * __distance__ distance between descriptors (scalar)
            %
            % ## Options
            % * __Mask__ default empty
            %   * In the first form, mask specifying permissible matches
            %   between an input query and train matrices of descriptors.
            %   Matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors,1)]`.
            %   * In the second form, set of masks. Each `masks{i}` specifies
            %   permissible matches between the input query descriptors and
            %   stored train descriptors from the i-th image
            %   `trainDescCollection{i}`. Cell array of length
            %   `length(trainDescriptors)`, each a matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors{i},1)]`.
            % * __CompactResult__ Parameter used when the mask (or masks) is
            %   not empty. If `CompactResult` is false, the `matches` vector
            %   has the same size as `queryDescriptors` rows. If
            %   `CompactResult` is true, the matches vector does not contain
            %   matches for fully masked-out query descriptors. default false
            %
            % This extended variant of cv.DescriptorMatcher.match method finds
            % several best matches for each query descriptor. The matches are
            % returned in the distance increasing order. See
            % cv.DescriptorMatcher.match for the details about query and train
            % descriptors.
            %
            % See also: cv.DescriptorMatcher.match,
            %  cv.DescriptorMatcher.radiusMatch, cv.batchDistance
            %
            matches = DescriptorMatcher_(this.id, 'knnMatch', ...
                queryDescriptors, varargin{:});
        end

        function matches = radiusMatch(this, queryDescriptors, varargin)
            %RADIUSMATCH  For each query descriptor, finds the training descriptors not farther than the specified distance
            %
            %     matches = matcher.radiusMatch(queryDescriptors, trainDescriptors, maxDistance)
            %     matches = matcher.radiusMatch(queryDescriptors, maxDistance)
            %     [...] = matcher.radiusMatch(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ Query set of descriptors.
            % * __trainDescriptors__ Train set of descriptors. This set is not
            %   added to the train descriptors collection stored in the class
            %   object.
            % * __maxDistance__ Threshold for the distance between matched
            %   descriptors. Distance here means metric distance (e.g. Hamming
            %   distance), not the distance between coordinates (which is
            %   measured in Pixels)!
            %
            % ## Output
            % * __matches__ Found matches. A cell array of length
            %   `size(queryDescriptors,1)`, each cell is a structure array
            %   that has the following fields:
            %   * __queryIdx__ query descriptor index (zero-based index)
            %   * __trainIdx__ train descriptor index (zero-based index)
            %   * __imgIdx__ train image index (zero-based index)
            %   * __distance__ distance between descriptors (scalar)
            %
            % ## Options
            % * __Mask__ default empty
            %   * In the first form, mask specifying permissible matches
            %   between an input query and train matrices of descriptors.
            %   Matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors,1)]`.
            %   * In the second form, set of masks. Each `masks{i}` specifies
            %   permissible matches between the input query descriptors and
            %   stored train descriptors from the i-th image
            %   `trainDescCollection{i}`. Cell array of length
            %   `length(trainDescriptors)`, each a matrix of size
            %   `[size(queryDescriptors,1),size(trainDescriptors{i},1)]`.
            % * __CompactResult__ Parameter used when the mask (or masks) is
            %   not empty. If `CompactResult` is false, the `matches` vector
            %   has the same size as `queryDescriptors` rows. If
            %   `CompactResult` is true, the matches vector does not contain
            %   matches for fully masked-out query descriptors. default false
            %
            % For each query descriptor, the methods find such training
            % descriptors that the distance between the query descriptor and
            % the training descriptor is equal or smaller than `maxDistance`.
            % Found matches are returned in the distance increasing order.
            %
            % See also: cv.DescriptorMatcher.match,
            %  cv.DescriptorMatcher.knnMatch, cv.batchDistance
            %
            matches = DescriptorMatcher_(this.id, 'radiusMatch',...
                queryDescriptors, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function Q1 = mpc(T_meas,T_sp)

s = 'http://byu.apmonitor.com';
c = 'mpc';

% input measurement
apm_meas(s,c,'TC',T_meas);

% input setpoint with deadband +/- DT
DT = 0.1;
apm_option(s,c,'TC.sphi',T_sp+DT);
apm_option(s,c,'TC.splo',T_sp-DT);

% solve MPC
output = apm(s,c,'solve');

% test for successful solution
if (apm_tag(s,c,'apm.appstatus')==1)
    % retrieve the first Q value
    Q1 = apm_tag(s,c,'Q1.Newval');
else
    % display output for debugging
    print(output)
    % not successful, set voltage to zero
    Q1 = 0;
end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function display(t)
%DISPLAY Command window display of a tenmat.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

disp(t,inputname(1));

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function value = error_f ( x )

%*****************************************************************************80
%
%% ERROR_F computes the error function.
%
%  Discussion:
%
%    This function was renamed "ERROR_F" from "ERF", to avoid a conflict
%    with the name of a corresponding routine often, but not always,
%    supplied as part of the math support library.
%
%    The definition of the error function is:
%
%      ERF(X) = ( 2 / SQRT ( PI ) ) * Integral ( 0 <= T <= X ) EXP ( -T**2 ) dT
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    FORTRAN77 original version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1989,
%    ISBN: 0-13-627258-4,
%    LC: TA345.K34.
%
%  Parameters:
%
%    Input, real X, the argument of the error function.
%
%    Output, real VALUE, the value of the error function at X.
%
  persistent nterf;
  persistent sqeps;
  persistent xbig;

  erfcs = [ ...
    -0.049046121234691808, -0.14226120510371364, ...
     0.010035582187599796, -0.000576876469976748, ...
     0.000027419931252196, -0.000001104317550734, ...
     0.000000038488755420, -0.000000001180858253, ...
     0.000000000032334215, -0.000000000000799101, ...
     0.000000000000017990, -0.000000000000000371, ...
     0.000000000000000007 ];

  sqrtpi = 1.7724538509055160;
%
%  Initialize the Chebyshev series.
%
  if ( size ( nterf ) == 0 )
    nterf = inits ( erfcs, 13, 0.1 * eps );
    xbig = sqrt ( - log ( sqrtpi * eps ) );
    sqeps = sqrt ( 2.0 * eps );
  end

  y = abs ( x );

  if ( y <= sqeps )
    value = 2.0 * x / sqrtpi;
  elseif ( y <= 1.0 )
    value = x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
  elseif ( y <= xbig )
    value = r8_sign ( x ) * ( 1.0 - error_fc ( y ) );
  else
    value = r8_sign ( x );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef SliderCheckboxChangedData < event.EventData & dynamicprops
    % Event data for widget property value changes
    %
    % Syntax:
    %           obj = wt.eventdata.SliderCheckboxChangedData(name,index,prop,state,value)
    %

    % Copyright 2020-2021 The MathWorks, Inc.

    %% Properties
    properties (SetAccess = protected)
        Name (1,1) string
        Index (1,1) double
        Property (1,1) string
        State (1,1) logical
        Value (1,1) double
    end %properties

    %% Constructor / destructor
    methods
        function obj = SliderCheckboxChangedData(name,index,prop,state,value)

            % Set the changed properties
            obj.Name = name;
            obj.Index = index;
            obj.Property = prop;
            obj.State = state;
            obj.Value = value;

        end %constructor
    end %methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [Lb]=logicRemoveInterior(L)

% function [Lb]=logicRemoveInterior(L)
% ------------------------------------------------------------------------
% This function removes the interior entries for the input logic L.
% Interior entries are those fully surrounded by neighbours (e.g. in 3D all
% entries which are attached to a top, bottom, left, right, front, and a
% back neighbour).
% Vectors and n-dimensional arrays are supported.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log:
% 2018/06/14 Created as alternative to bwmorph3 which requires a special
% toolbox and is from R2018a.
%------------------------------------------------------------------------

if ~isempty(L)
    if isvector(L) %Handle vector
        h=ones(1,3);
        if ~isrow(L)
            h=h'; %Transpose for column
        end
    else
        nd=ndims(L);
        h=zeros(3*ones(1,nd)); %Initialize h
        ijkMid=2*ones(1,nd); %Middle coordinate
        I=eye(nd,nd); %Identify matrix to offset indices
        ijkFilter=[ijkMid; ijkMid(ones(nd,1),:)+I; ijkMid(ones(nd,1),:)-I]; %Indices for filer
        ind=sub2indn(size(h),ijkFilter);%Linear indices for filer
        h(ind)=1; %Set ones for ND cross shape
    end
    L=L>0;
    LC=convn(double(L),h,'same');
    Lb=L & ~(LC==sum(h(:)));
else
    Lb=[];
end

%%
% _*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 dfupdateppdists(dffig)
%DFUPDATEPPDISTS Update distribution list for probability plots

%   $Revision: 1.1.6.5 $  $Date: 2004/01/24 09:36:02 $
%   Copyright 2003-2004 The MathWorks, Inc.

if nargin<1 || isempty(dffig)
   dffig = dfgetset('dffig');
end

% Get handle to control containing the distribution list
hsel = getappdata(dffig,'selectioncontrols');

% Determine plot type, don't continue unless it's a probability plot
h = hsel(3);   % handle of display type control
choice = get(hsel(3),'Value');
ftypes = getappdata(hsel(3),'codenames');
ftype = ftypes{choice};
if ~isequal(ftype, 'probplot');
   return
end

% Look at all plotted data sets and check for negative data
dsdb = getdsdb;
dslist = find(dsdb);
dslist(dslist==dsdb) = [];
nds = length(dslist);
xmin = 1;
for j=1:nds
   ds = dslist(j);
   if ds.plot
      xlimits = xlim(ds);
      if (xlimits(1) < xmin)
         xmin = xlimits(1);
         if (xmin<0)
            break
         end
      end
   end
end

% Get the current selection, and try to re-select it later
h = hsel(5);   % handle of distribution control
dlist = get(h, 'String');
choice = get(h, 'Value');
if choice<=length(dlist)
   cursel = dlist{choice};
else
   cursel = [];
end

% Get the distribution entries only
dlist = getappdata(h, 'allfullnames');
codenames = getappdata(h, 'allcodenames');
ndists = length(dlist);

% With 0 or negative data showing, omit some distributions
if (xmin<=0)
   ok = true(ndists,1);
   allspec = getappdata(h,'alldistspec');

   % ------- Could vectorize the following
   for j=1:ndists
      spec = allspec(j);
      ok(j) = checkdist(spec,xmin);
   end
   % -------

   dlist = dlist(ok);
   codenames = codenames(ok);
   ndists = length(dlist);
end
setappdata(h,'okcodenames',codenames);

% Create a list of fit names
fitdb = getfitdb;
flist = find(fitdb);
flist(flist==fitdb) = [];
nfits = length(flist);

% Weed out fits that cannot be used for probability plotting
for j=nfits:-1:1
   fj = flist(j);
   if ~fj.isgood || isequal(fj.fittype,'smooth')
      flist(j) = [];
   else
      distspec = fj.distspec;
      if isempty(distspec) || ~distspec.iscontinuous ...
                           || ~checkdist(distspec,xmin)
         flist(j) = [];
      end
   end
end
nfits = length(flist);

% Create combined list
newdlist = cell(ndists + nfits, 1);
newdlist(1:ndists) = dlist(:);
savelist = cell(nfits,1);
for j=1:nfits
   fj = flist(j);
   newdlist{ndists+j} = sprintf('Estimated %s (%s)',fj.distspec.name,fj.name);
   savelist{j} = fj.name;
end
setappdata(h, 'fitnames', savelist);

% Re-select the previous selection
set(h, 'String', newdlist);
choice = 1;
if ~isempty(cursel)
   choices = strmatch(cursel,newdlist,'exact');
   if numel(choices) == 1
      choice = choices;
   end
end
set(h, 'Value', choice);

% -----------------------------
function ok = checkdist(spec,xmin)
%Check distribution against minimum data value

lobnd = spec.support(1);
strict = ~spec.closedbound(1);
if strict && lobnd>=xmin
   ok = false;
elseif ~strict && lobnd>xmin
   ok = false;
else
   ok = true;
end

author: zouchuhang

repo: LayoutNet


text:
function op = prox_l1linf( q )

%PROX_L1LINF    L1-LInf block norm: sum of L-inf norms of rows.
%    OP = PROX_L1LINF( q ) implements the nonsmooth function
%        OP(X) = q * sum_{i=1:m} norm(X(i,:),Inf)
%    where X is a m x n matrix.  If n = 1, this is equivalent
%    to PROX_L1
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    then it must be positive and real.
%    If Q is a vector, it must be m x 1, and in this case,
%    the weighted norm OP(X) = sum_{i} Q(i)*norm(X(i,:),Inf)
%    is calculated.

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || any(q <= 0),
	error( 'Argument must be positive.' );
end
op = tfocs_prox( @(x)f(x,q), @(x,t)prox_f(x,t,q), 'vector' );
end % end of main function

function v = f(x,q)
    if numel(q) ~= 1 && size(q,1) ~= size(x,1)
        error('Weight must be a scalar or a column vector');
    end
    v = sum( q.* max(abs(x),[],2) );
end

function x = prox_f(x,t,q)
  if nargin < 2,
      error( 'Not enough arguments.' );
  end

  [n,d] = size(x);
  dim   = 2;

  % Option 1: explicitly call prox_linf on the rows:
  % slow, but the chief benefit is that it is low memory
  % This would probably be faster if the matrix was transposed before and after
%   for k= 1:n
%       if isscalar(q), qk = q;
%       else, qk = q(k);
%       end
%       x(k,:) = prox_linf_q( qk, x(k,:).', t ).';
%   end
%   return;

% Option 2: vectorize the call.  By far, more efficient than option 1

  %s     = sort( abs(x), dim, 'descend' );
  %cs    = cumsum(s,dim);
  % Since Matlab stores matrices in column-major order, this method is more cache friendly:
  s     = sort( abs(x)', q, 'descend' );
  cs    = cumsum(s,1)';
  s     = s';

  s     = [s(:,2:end), zeros(n,1)];

  ndx1 = zeros(n,1);
  ndx2 = zeros(n,1);

  if isscalar(q),
      tq = repmat( t*q, n, d );
  else
      tq = repmat( t*q, 1, d );
  end
  %Z = cs - s*diag(1:d);
  % The above may require a lot of memory, so use spdiag or this:
  Z = cs - bsxfun(@times,s,1:d);

  Z = ( Z >= tq );
  Z = Z.';
  % now Z is d x n (typically d is large, n is small)

  % Not sure how to vectorize the find.
  % One option is to use the [i,j] = find(...) form,
  % but that's also extra work, since we can't just find the "first".
  if n > 5
      % avoid the for-loop
      ndx1  = (d+1)-sum(Z)'; % this is the first row with Z > 0
      ndx2  = ndx1;
      ndx2( ndx2 == d+1 ) = Inf;
      ndx1( ndx1 == d+1)  = d; % arbitrary, but do this so we don't have a special case later
  else
      % this might be slightly less memory, so keep the code
      for k = 1:n
          % This is why we transposed Z: due to column-major order,
          %     find( columnVector ) is faster than find( rowVector )
          ndxk = find( Z(:,k), 1 );
          if ~isempty(ndxk)
              ndx1(k) = ndxk;
              ndx2(k) = ndxk;
          else
              ndx1(k) = 1; % value doesn't matter
              ndx2(k) = Inf;
          end
      end
  end
  indx_cs = sub2ind( [n,d], (1:n)', ndx1 );
  tau = (cs(indx_cs) - tq(:,1))./ndx2;
  tau = repmat( tau, 1, d );
  tau_noZeros = tau;
  tau_noZeros( ~x ) = 1;
  x   = x .* (  tau ./ max( abs(x), tau_noZeros ) );
  % Another, but not really better, way is to not do the rempat stuff and do:
  % x     = sign(x).*bsxfun( @min, tau, abs(x) );

end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% Copyright (C) 2013-2016 Nathanael Perraudin, Johan Paratte, David I Shuman.
% This file is part of GSPbox version {VERSION}
%
% 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/>.

% If you use this toolbox please kindly cite
%     N. Perraudin, J. Paratte, D. Shuman, V. Kalofolias, P. Vandergheynst,
%     and D. K. Hammond. GSPBOX: A toolbox for signal processing on graphs.
%     ArXiv e-prints, Aug. 2014.
% http://arxiv.org/abs/1408.5781

repo: gspbox


text:
%TRINTERP2 Interpolate SE(2) homogeneous transformations
%
% TRINTERP2(T0, T1, S) is a homogeneous transform (3x3) interpolated
% between T0 when S=0 and T1 when S=1.  T0 and T1 are both homogeneous
% transforms (4x4).  If S (Nx1) then T (3x3xN) is a sequence of
% homogeneous transforms corresponding to the interpolation values in S.
%
% TRINTERP2(T1, S) as above but interpolated between the identity matrix
% when S=0 to T1 when S=1.
%
% TRINTERP2(T0, T1, M) as above but M is a positive integer and return a
% sequence (4x4xM) of homogeneous transforms linearly interpolating between
% T0 and T1 in M steps.
%
% TRINTERP2(T1, M) as above but return a sequence (4x4xM) of
% homogeneous interpolating between identity matrix and T1 in M steps.
%
% Notes::
% - T0 or T1 can also be an SO(2) rotation matrix (2x2).
% - Rotation angle is linearly interpolated.
% - To obtain smooth continuous motion S should also be smooth and continuous,
%   such as computed by tpoly or lspb.
%
% See also TRINTERP, SE3.interp, UnitQuaternion, TPOLY, LSPB.

%## 2d homogeneous trajectory

% 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 T = trinterp2(A, B, C)

    switch nargin
        case 2
            % trinterp(T, s)
            T1 = A; s = B;

            th0 = 0;
            th1 = atan2(T1(2,1), T1(1,1));
            if ~isrot2(T1)
                p0 = [0 0]';
                p1 = transl2(T1);
            end
        case 3
            % trinterp(T1, T2, s)
            T0 = A; T1 = B; s = C;
            assert(all(size(A) == size(B)), 'SMTB:trinterp2:badarg', '2 matrices must be same size');
            th0 = atan2(T0(2,1), T0(1,1));
            th1 = atan2(T1(2,1), T1(1,1));
            if ~isrot2(T0)
                p0 = transl2(T0);
                p1 =transl2(T1);
            end
        otherwise
            error('SMTB:trinterp2:badarg', 'must be 2 or 3 arguments');
    end

    if length(s) == 1 && s > 1 && (s == floor(s))
        % integer value
        s = [0:(s-1)] / (s-1);
    elseif any(s<0 | s>1)
        error('SMTB:trinterp2:badarg', 'values of S outside interval [0,1]');
    end

    if isrot2(T1)

        % SO(2) case
        for i=1:length(s)
            th = th0*(1-s(i)) + s(i)*th1;

            T(:,:,i) = rot2(th);
        end
    else
        % SE(2) case
        for i=1:length(s)
            th = th0*(1-s(i)) + s(i)*th1;
            pr = p0*(1-s(i)) + s(i)*p1;

            T(:,:,i) = rt2tr(rot2(th), pr);
        end
    end

end

author: petercorke

repo: spatialmath-matlab


text:
function spm_dem_cue_movie(DEM,q)
% creates a movie of cued pointing
% FORMAT spm_dem_cue_movie(DEM,q)
%
% DEM - DEM structure from reaching simulations
% q   - flag switching from true to perceived reaching
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dem_set_movie.m 4231 2011-03-07 21:00:02Z karl $

% Dimensions
%--------------------------------------------------------------------------
N    = size(DEM.pU.v{1},2);
n    = size(DEM.pP.P{1},2);

% evaluate true location (targets)
%----------------------------------------------------------------------
for i = 1:N
    L(:,:,i) = DEM.pP.P{1};
end

if nargin > 1

    % evaluate perceived positions (motor plant)
    %----------------------------------------------------------------------
    x    = tan(DEM.qU.v{1}(1:2,:));                  % finger location
    c    = DEM.qU.v{1}(4 + (1:n),:);                 % target contrast

else

    % evaluate true positions (motor plant)
    %----------------------------------------------------------------------
    x    = tan(DEM.pU.v{1}(1:2,:));                  % finger location
    c    = DEM.pU.v{1}(4 + (1:n),:);                 % target contrast

end
c    = c - min(c(:)) + 1/32;
c    = c/max(c(:));
fin  = imread('pointfinger.jpg');

% movie
%--------------------------------------------------------------------------
s     = 2;
for i = 1:N

    cla
    axis image ij
    hold on

    % finger
    %----------------------------------------------------------------------
    imagesc(([-1 0] + .68)*s + x(1,i),([-1 0] + .96)*s + x(2,i),fin);
    hold on

    % trajectory
    %----------------------------------------------------------------------
    plot(x(1,1:i),x(2,1:i),'k:')

    % targets
    %----------------------------------------------------------------------
    for j = 1:n
        plot(L(1,j,i),L(2,j,i),'.','MarkerSize',64,'color',[c(j,i) (1 - c(j,i)) 0])
    end
    axis([-1 1 -1 1]*2)
    hold off
    drawnow

    % save
    %----------------------------------------------------------------------
    M(i) = getframe(gca);

end

% set ButtonDownFcn
%--------------------------------------------------------------------------
h = findobj(gca,'type','image');
set(h(1),'Userdata',{M,16})
set(h(1),'ButtonDownFcn','spm_DEM_ButtonDownFcn')

author: spm


text:
function [lineextimg c] = orient_from_lines(lines, vp, imgwidth, imgheight)

%%
ls = sample_line(lines);
linesamples = cat(1,ls(:).sample);
linesampleclass = cat(1,ls(:).lineclass);

%%
% lineextimg = cell(3,3);
% for i=1:9, lineextimg{i} = zeros(imgheight,imgwidth); end
lineextimg = cell(3,3,2);
for i=1:18, lineextimg{i} = zeros(imgheight,imgwidth); end

%%
% poly = extend_line(line, vp{1}, stoppinglines_sample, imgwidth, imgheight);
for i = 1:length(lines)
    lc = lines(i).lineclass;
    if lc~=0
    for extdir = setdiff(1:3, lc)
        targetdir = setdiff(1:3, [lc extdir]);

%         poly = extend_line_old(lines(i), vp{extdir}, linesamples(linesampleclass==targetdir,:), imgwidth, imgheight);
%         lineextimg{lc,extdir} = lineextimg{lc,extdir} + poly2mask(poly(:,1), poly(:,2), imgheight, imgwidth);

        poly = extend_line(lines(i), vp{extdir}, 1, linesamples(linesampleclass==targetdir,:), imgwidth, imgheight);
        lineextimg{lc,extdir,1} = lineextimg{lc,extdir,1} + poly2mask(poly(:,1), poly(:,2), imgheight, imgwidth);

        poly = extend_line(lines(i), vp{extdir}, -1, linesamples(linesampleclass==targetdir,:), imgwidth, imgheight);
        lineextimg{lc,extdir,2} = lineextimg{lc,extdir,2} + poly2mask(poly(:,1), poly(:,2), imgheight, imgwidth);
    end
    end
end

%%
function poly = extend_line(line, vp, towards_or_away, stoppinglines_sample, imgwidth, imgheight)
% towards_or_away: 1 or -1

p1 = line.point1;
p2 = line.point2;

curp1 = p1; curp2 = p2;
move_amount = 128;
while move_amount>=1
    [newp1 newp2 atvp] = move_line_towards_vp(curp1, curp2, vp, towards_or_away * move_amount);

    failed = 0;
    if atvp==1
%         move_amount = 0; % exit now.
        failed = 1;

    elseif (newp1(1)>imgwidth || newp1(1)<1 || newp1(2)>imgheight || newp1(2)<1) && ...
       (newp2(1)>imgwidth || newp2(1)<1 || newp2(2)>imgheight || newp2(2)<1)
        failed = 1;

    else
        isstop = inpolygon(stoppinglines_sample(:,1), stoppinglines_sample(:,2), ...
            [p1(1) p2(1) newp2(1) newp1(1) p1(1)], [p1(2) p2(2) newp2(2) newp1(2) p1(2)]);

        if any(isstop)
            failed = 1;
        end
    end

    if failed
        move_amount = move_amount/2;
    else
        curp1 = newp1;
        curp2 = newp2;
    end
end
% poly = [curp1(:)'; curp2(:)'];
poly = [p1(:)'; p2(:)'; curp2(:)'; curp1(:)'];

% curp1 = p1; curp2 = p2;
% move_amount = 32;
% while move_amount>=1
%     [newp1 newp2 atvp] = move_line_towards_vp(curp1, curp2, vp, -move_amount);
%
%     failed = 0;
%     if atvp==1
%         move_amount = 0; % exit now.
%
%     elseif (newp1(1)>imgwidth || newp1(1)<1 || newp1(2)>imgheight || newp1(2)<1) && ...
%        (newp2(1)>imgwidth || newp2(1)<1 || newp2(2)>imgheight || newp2(2)<1)
%         failed = 1;
%
%     else
%         isstop = inpolygon(stoppinglines_sample(:,1), stoppinglines_sample(:,2), ...
%             [p1(1) p2(1) newp2(1) newp1(1) p1(1)], [p1(2) p2(2) newp2(2) newp1(2) p1(2)]);
%
%         if any(isstop)
%             failed = 1;
%         end
%     end
%
%     if failed
%         move_amount = move_amount/2;
%     else
%         curp1 = newp1;
%         curp2 = newp2;
%     end
% end
% poly = [poly; curp2(:)'; curp1(:)'];
% poly = [poly; poly(1,:)];

%%

%%
function [newp1 newp2 atvp] = move_line_towards_vp(linep1, linep2, vp, amount)

% d = dist_line_to_point(linep1, linep2, vp);
% r = amount / d;
n1 = norm(vp-linep1);
n2 = norm(vp-linep2);
dir1 = (vp - linep1) / n1;
dir2 = (vp - linep2) / n2;
ratio21 = n2 / n1;

% if n1>amount && n2<amount
%     fprintf('check');
% end

if n1 < amount
    newp1 = linep1;
    newp2 = linep2;
    atvp = 1;
else
    newp1 = linep1 + dir1 * amount;
    newp2 = linep2 + dir2 * amount * ratio21;
    atvp = 0;
end

%%
function [newp1 newp2 atvp] = move_line_towards_vp_old(linep1, linep2, vp, amount)

d = dist_line_to_point(linep1, linep2, vp);
r = amount / d;
if r > 1
    newp1 = vp;
    newp2 = vp;
    atvp = 1;
else
    newp1 = linep1 + (vp-linep1)*r;
    newp2 = linep2 + (vp-linep2)*r;
    atvp = 0;
end


text:
% Summarise EEG events
% axs Jun 2018
function [ev_count] = event_code_count(EEG)

evt = struct2table(EEG.event);
evtype = evt.type;

n_ev = numel(evtype);

ev2 = zeros(n_ev,1);

str_flag = 0;
for i = 1:n_ev

    ev_here = evt.type{i};

    ev_num = str2double(ev_here);

    if isnan(ev_num)
        if strcmp(ev_here,'boundary') == 1
            ev_num = -99;
        else
            ev_num = 0;
            str_flag = 1;
        end
    end

    ev2(i) = ev_num;
end

% Now parse and display
n_ev_types = unique(ev2);

ev_count = zeros(numel(n_ev_types),2);
ev_count(1:end,1) = n_ev_types;

for i = 1:numel(n_ev_types)
    indx_here = find(ev2 == n_ev_types(i));
    count_here = numel(indx_here);
    ev_count(i,2) = count_here;
end

disp('Your event code counts are:')
disp(ev_count)

if str_flag == 1
    disp('Event codes of unknown strings were taken as 0s');
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% Function MYGRID_SAND	Read bathymetry data from Sandwell Database
%      [image_data,vlat,vlon] = mygrid_sand(region,iopt)
%
% program to get bathymetry from topo_8.2.img  (Smith and Sandwell bathymetry)
%  (values are even numbered if interpolated, odd-numbered if from a ship sounding)
% WARNING: change DatabasesDir to the correct one for your machine
%						Catherine de Groot-Hedlin
% latitudes must be between -72.006 and 72.006;
%	input:
%		region =[south north west east];
%               iopt = 1 if bathymetry is wanted (default)
%		       2 if ship tracks are wanted
%	output:
%		image_data
%                (for iopt = 1) - matrix of sandwell bathymetry/topography
%                (for iopt = 2) - matrix of ones and zeros, where 1 represents
%                    a ship location, 0 represents a depth based on interpolation
%		vlat - vector of latitudes associated with image_data
%      		vlon - vector of longitudes

%
function  [image_data,vlat,vlon] = mygrid_sand(region,iopt)

    DatabasesDir = '/export/home/grunion/cdh/airforce/sandwell.d';

    % determine the requested region
    blat = region(1);
    tlat = region(2);
    wlon = region(3);
    elon = region(4);

    % Setup the parameters for reading Sandwell data
    db_res         = 2/60;		% 2 minute resolution
    db_loc         = [-72.006 72.006 0.0 360-db_res];
    db_size        = [6336 10800];
    nbytes_per_lat = db_size(2)*2;	% 2-byte integers
    image_data     = [];

    % Determine if the database needs to be read twice (overlapping prime meridian)
    if ((wlon<0)&(elon>=0))
        wlon      = [wlon           0];
        elon      = [360-db_res  elon];
    end

    % Calculate number of "records" down to start (latitude) (0 to db_size(1)-1)
    % (mercator projection)
    rad=pi/180;arg1=log(tan(rad*(45+db_loc(1)/2)));
    arg2=log(tan(rad*(45+blat/2)));
    iblat = fix(db_size(1) +1 - (arg2-arg1)/(db_res*rad))

    arg2=log(tan(rad*(45+tlat/2)));
    itlat = fix(db_size(1) +1 - (arg2-arg1)/(db_res*rad))

    if (iblat < 0 ) | (itlat > db_size(1)-1)
        errordlg([' Requested latitude is out of file coverage ']);
    end

    % Go ahead and read the database
    for i = 1:length(wlon);

        % Open the data file
        fid = fopen([DatabasesDir '/topo_8.2.img'], 'r');
        if (fid < 0)
            errordlg(['Could not open database: ' DatabasesDir '/topo_8.2.img'],'Error');
        end

        % Make sure the longitude data goes from 0 to 360
        if wlon(i) < 0
            wlon(i) = 360 + wlon(i);
        end

        if elon(i) < 0
            elon(i) = 360 + elon(i);
        end

        % Calculate the longitude indices into the matrix (0 to db_size(1)-1)
        iwlon(i) = fix((wlon(i)-db_loc(3))/db_res)
        ielon(i) = fix((elon(i)-db_loc(3))/db_res)
        if (iwlon(i) < 0 ) | (ielon(i) > db_size(2)-1)
            errordlg([' Requested longitude is out of file coverage ']);
        end

        % allocate memory for the data
        data = zeros(iblat-itlat+1,ielon(i)-iwlon(i)+1);

        % Skip into the appropriate spot in the file, and read in the data
        disp('Reading in bathymetry data');
        for ilat = itlat:iblat
            offset = ilat*nbytes_per_lat + iwlon(i)*2;
            status = fseek(fid, offset, 'bof');
            data(iblat-ilat+1,:)=fread(fid,[1,ielon(i)-iwlon(i)+1],'integer*2');
        end

        % close the file
        fclose(fid);

        % put the two files together if necessary
        if (i>1)
            image_data = [image_data data];
        else
            image_data = data;
        end
    end

    % Determine the coordinates of the image_data
    vlat=zeros(1,iblat-itlat+1);
    arg2 = log(tan(rad*(45+db_loc(1)/2.)));
    for ilat=itlat+1:iblat+1;
        arg1 = rad*db_res*(db_size(1)-ilat+0.5);
        term=exp(arg1+arg2);
        vlat(iblat-ilat+2)=2*atan(term)/rad -90;
    end
    vlon=db_res*((iwlon+1:ielon+1)-0.5);
    % now choose between bathymetry and ship track
    if iopt==2
        image_data=mod(image_data,2);
    end
    % to plot it up
    if iopt ==2
        imagesc(vlon,vlat,image_data),axis('xy'),colormap(1-gray)
        title('ship track soundings')
    else
        imagesc(vlon,vlat,image_data),axis('xy'),colormap(jet),colorbar('vert')
        title('Smith and Sandwell bathymetry')
    end
    xlabel('longitude'),ylabel('latitude')