text:
% TEST_CNNIMAGERETRIEVAL  Code to evaluate (not train) the methods presented in the papers:
% F. Radenovic, G. Tolias, O. Chum, Fine-tuning CNN Image Retrieval with No Human Annotation, TPAMI 2018
% F. Radenovic, G. Tolias, O. Chum, CNN Image Retrieval Learns from BoW: Unsupervised Fine-Tuning with Hard Examples, ECCV 2016
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

clear;

%---------------------------------------------------------------------
% Set data folder and testing parameters
%---------------------------------------------------------------------

% Set data folder, change if you have downloaded the data somewhere else
data_root = fullfile(get_root_cnnimageretrieval(), 'data');
% Check, and, if necessary, download train data (for whiten)
download_train(data_root);
% Check, and, if necessary, download test data (Oxf5k and Par6k), and fine-tuned networks
download_test(data_root);

% Set test options
test_datasets = {'oxford5k', 'paris6k', 'roxford5k', 'rparis6k'};  % list of datasets to evaluate on
test_imdim = 1024;  % choose test image dimensionality
use_ms = 1; % use multi-scale representation, otherwise use single-scale
use_rvec = 0;  % use regional representation (R-MAC, R-GeM), otherwise use global (MAC, GeM)
use_gpu = [1];  % use GPUs (array of GPUIDs), if empty use CPU

% Choose ECCV16 fine-tuned CNN network
% network_file = fullfile(data_root, 'networks', 'retrieval-SfM-30k', 'retrievalSfM30k-siamac-alex.mat');
% network_file = fullfile(data_root, 'networks', 'retrieval-SfM-30k', 'retrievalSfM30k-siamac-vgg.mat');

% Choose TPAMI18 fine-tuned CNN network
% network_file = fullfile(data_root, 'networks', 'retrieval-SfM-30k', 'retrievalSfM30k-gem-alex.mat');
% network_file = fullfile(data_root, 'networks', 'retrieval-SfM-120k', 'retrievalSfM120k-gem-vgg.mat');
network_file = fullfile(data_root, 'networks', 'retrieval-SfM-120k', 'retrievalSfM120k-gem-resnet101.mat');

% After running the training script train_cnnimageretrieval.m you can evaluate fine-tuned network
% network_file = fullfile(data_root, 'networks', 'exp', 'resnet101_gem_test', 'net-epoch-30');

%---------------------------------------------------------------------
% Set dependent variables
%---------------------------------------------------------------------

% Choose training data for whitening and set up data folder
% train_whiten_file = fullfile(data_root, 'train', 'dbs', 'retrieval-SfM-30k-whiten.mat'); % less images, faster
train_whiten_file = fullfile(data_root, 'train', 'dbs', 'retrieval-SfM-120k-whiten.mat'); % more images, better results but slower

% Set folder where original training images are stored, whitening learned on them
ims_whiten_dir = fullfile(data_root, 'train', 'ims');

% Prepare function for desc extraction
if ~use_rvec
	if ~use_ms
		descfun = @(x, y) cnn_vecms (x, y, 1);
	else
		descfun = @(x, y) cnn_vecms (x, y, [1, 1/sqrt(2), 1/2]);
	end
else
	if ~use_ms
		descfun = @(x, y) cnn_vecrms (x, y, 3, 1);
	else
		descfun = @(x, y) cnn_vecrms (x, y, 3, [1, 1/sqrt(2), 1/2]);
	end
end

%---------------------------------------------------------------------
% Testing
%---------------------------------------------------------------------
[~, network_name, ~] = fileparts(network_file);
fprintf('>> %s: Evaluating CNN image retrieval...\n', network_name);

% Load pre-trained CNN network
load(network_file);
net = dagnn.DagNN.loadobj(net);

% prepare GPUs if necessary
numGpus = numel(use_gpu);
if numGpus, fprintf('>> Preparing GPU(s)...\n'); end
if numGpus > 1
	% check parallel pool integrity as it could have timed out
	pool = gcp('nocreate');
	if ~isempty(pool) && pool.NumWorkers ~= numGpus
		delete(pool);
	end
	pool = gcp('nocreate');
	if isempty(pool)
		parpool('local', numGpus);
	end
end
if numGpus >= 1
	if numGpus == 1
		gpuinfo = gpuDevice(use_gpu);
		net.move('gpu');
		fprintf('>>>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
	else
		spmd
			gpuinfo = gpuDevice(use_gpu(labindex));
			fprintf('>>>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
		end
	end
end

% Load training data filenames and pairs for whitening
train_whiten = load(train_whiten_file);
if isfield(train_whiten, 'train') && isfield(train_whiten, 'val')
	cids  = [train_whiten.train.cids train_whiten.val.cids];
	qidxs = [train_whiten.train.qidxs train_whiten.val.qidxs+numel(train_whiten.train.cids)]; % query indexes
	pidxs = [train_whiten.train.pidxs train_whiten.val.pidxs+numel(train_whiten.train.cids)]; % positive indexes
else
	cids  = train_whiten.cids;
	qidxs = train_whiten.qidxs; % query indexes
	pidxs = train_whiten.pidxs; % positive indexes
end

% learn whitening
fprintf('>> whitening: Extracting CNN descriptors for training images...\n');
vecs_whiten = cell(1, numel(cids));
if numGpus <= 1
	progressbar(0);
	for i=1:numel(cids)
		vecs_whiten{i} = descfun(imresizemaxd(imread(cid2filename(cids{i}, ims_whiten_dir)), test_imdim, 0), net);
		progressbar(i/numel(cids));
	end
else
	time = tic;
	parfor i=1:numel(cids)
		if strcmp(net.device, 'cpu'), net.move('gpu'); end
		vecs_whiten{i} = descfun(imresizemaxd(imread(cid2filename(cids{i}, ims_whiten_dir)), test_imdim, 0), net);
	end
	fprintf('>>>> done in %s\n', htime(toc(time)));
end
vecs_whiten = cell2mat(vecs_whiten);
fprintf('>> whitening: Learning...\n');
Lw = whitenlearn(vecs_whiten, qidxs, pidxs);

% extract and evaluate
for d = 1:numel(test_datasets)
	fprintf('>> %s: Processing test dataset...\n', test_datasets{d});
	cfg = configdataset (test_datasets{d}, fullfile(data_root, 'test/')); % config file for the dataset

	fprintf('>> %s: Extracting CNN descriptors for db images...\n', test_datasets{d});
	vecs = cell(1, cfg.n);
	if numGpus <= 1
		progressbar(0);
		for i = 1:cfg.n
			vecs{i} = descfun(imresizemaxd(imread(cfg.im_fname(cfg, i)), test_imdim, 0), net);
			progressbar(i/cfg.n);
		end
	else
		time = tic;
		parfor i = 1:cfg.n
			if strcmp(net.device, 'cpu'), net.move('gpu'); end
			vecs{i} = descfun(imresizemaxd(imread(cfg.im_fname(cfg, i)), test_imdim, 0), net);
		end
		fprintf('>>>> done in %s\n', htime(toc(time)));
	end
	vecs = cell2mat(vecs);

	fprintf('>> %s: Extracting CNN descriptors for query images...\n', test_datasets{d});
	qvecs = cell(1, cfg.nq);
	if numGpus <= 1
		progressbar(0);
		for i = 1:cfg.nq
			qvecs{i} = descfun(crop_qim(imread(cfg.qim_fname(cfg, i)), cfg.gnd(i).bbx, test_imdim), net);
			progressbar(i/cfg.nq);
		end
	else
		time = tic;
		parfor i = 1:cfg.nq
			if strcmp(net.device, 'cpu'), net.move('gpu'); end
			qvecs{i} = descfun(crop_qim(imread(cfg.qim_fname(cfg, i)), cfg.gnd(i).bbx, test_imdim), net);
		end
		fprintf('>>>> done in %s\n', htime(toc(time)));
	end
	qvecs = cell2mat(qvecs);

	vecsLw = whitenapply(vecs, Lw.m, Lw.P); % apply whitening on database descriptors
	qvecsLw = whitenapply(qvecs, Lw.m, Lw.P); % apply whitening on query descriptors

	fprintf('>> %s: Retrieval...\n', test_datasets{d});
	if strcmp(test_datasets{d}, 'oxford5k') || strcmp(test_datasets{d}, 'paris6k')
		% % raw descriptors
		% sim = vecs'*qvecs;
		% [sim, ranks] = sort(sim, 'descend');
		% map = compute_map (ranks, cfg.gnd);
		% fprintf('>> %s: mAP = %.4f, without whiten\n', test_datasets{d}, map);
		% with learned whitening
		sim = vecsLw'*qvecsLw;
		[sim, ranks] = sort(sim, 'descend');
		map = compute_map (ranks, cfg.gnd);
		fprintf('>> %s: mAP = %.4f\n', test_datasets{d}, map);
	elseif strcmp(test_datasets{d}, 'roxford5k') || strcmp(test_datasets{d}, 'rparis6k')
		sim = vecsLw'*qvecsLw;
		[sim, ranks] = sort(sim, 'descend');
		% evaluate ranks
		ks = [1, 5, 10];
		% search for easy (E setup)
		for i = 1:numel(cfg.gnd), gnd(i).ok = [cfg.gnd(i).easy]; gnd(i).junk = [cfg.gnd(i).junk, cfg.gnd(i).hard]; end
		[mapE, apsE, mprE, prsE] = compute_map (ranks, gnd, ks);
		% search for easy & hard (M setup)
		for i = 1:numel(cfg.gnd), gnd(i).ok = [cfg.gnd(i).easy, cfg.gnd(i).hard]; gnd(i).junk = cfg.gnd(i).junk; end
		[mapM, apsM, mprM, prsM] = compute_map (ranks, gnd, ks);
		% search for hard (H setup)
		for i = 1:numel(cfg.gnd), gnd(i).ok = [cfg.gnd(i).hard]; gnd(i).junk = [cfg.gnd(i).junk, cfg.gnd(i).easy]; end
		[mapH, apsH, mprH, prsH] = compute_map (ranks, gnd, ks);
		fprintf('>> %s: mAP E: %.2f, M: %.2f, H: %.2f\n', test_datasets{d}, 100*mapE, 100*mapM, 100*mapH);
		fprintf('>> %s: mP@k[%d %d %d] E: [%.2f %.2f %.2f], M: [%.2f %.2f %.2f], H: [%.2f %.2f %.2f]\n', test_datasets{d}, ks(1), ks(2), ks(3), 100*mprE, 100*mprM, 100*mprH);
	end
end

author: filipradenovic

repo: cnnimageretrieval


text:
function [precision,recall] = evaluate(NAME,boundingBoxes, threshold)
tmpName = extractBefore(NAME,'.jpg');
PATH_REAL = "Data/" + tmpName + ".txt";
fileID = fopen(PATH_REAL{1},'r');
formatSpec = ' %d ';
BoundingData = fscanf(fileID,formatSpec);

[realObjectCount,~] = size(BoundingData);
realObjectCount = (realObjectCount - 1)/4;
correctObjects  = 0;
[~,detectedCount] = size(boundingBoxes);

for i=1:detectedCount
    ourBox = boundingBoxes{1,i};
    topX = ourBox(1,1);
    topY = ourBox(1,2);
    lowX = ourBox(1,3);
    lowY = ourBox(1,4);
    for b=1:realObjectCount
        rtopX = BoundingData(1 + (4*(b-1)),1);
        rtopY = BoundingData(2 + (4*(b-1)),1);
        rlowX = BoundingData(3 + (4*(b-1)),1);
        rlowY = BoundingData(4 + (4*(b-1)),1);
        width = 0;
        if lowX > topX
            width =  lowX-topX;
        else
             width =  topX-lowX;
        end
        height = 0;
        if lowY > topY
            height =  lowY-topY ;
        else
             height =  topY-lowY;
        end

        realWidth = 0;
        if rlowX > rtopX
            realWidth =  rlowX-rtopX;
        else
             realWidth =  rtopX-rlowX;
        end
        realHeight = 0;
        if rlowY > rtopY
            realHeight =   rlowY-rtopY ;
        else
             realHeight =  rtopY-rlowY;
        end

        ourBox  = [ topX,   topY, width ,height  ];
        realBox = [ rtopX, rtopY,realWidth ,realHeight];
        overlapRatio = bboxOverlapRatio(ourBox,realBox);
        if overlapRatio > threshold
            correctObjects = correctObjects + 1;
            break;
        end
    end
end

recall = correctObjects /realObjectCount ;
precision = correctObjects /detectedCount ;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/evaluate.m


text:
function [DI]=TriScatteredInterp_ND(DT,D,XI,InterpMethod)

% function [DI]=TriScatteredInterp_ND(DT,D,XI,InterpMethod)
% ------------------------------------------------------------------------
%
%
%
% ------------------------------------------------------------------------

%%

if ~isa(DT,'DelaunayTri') %if DT is not a delaunay tesselation
    DT=delaunayTriangulation(DT); %assuming DT are coordinates, replace by Delaunay tesselation
end

DI=nan(size(XI,1),size(D,2)); %Allocate DI
for q=1:size(D,2)% loop over dimensions
    switch InterpMethod
        case 'nat_near' %natural in chull, neirest outside chull
            [DI(:,q),~]=TriScatteredInterp_nat_near(DT,D(:,q),XI);
        otherwise %TriScatterInterp can handle other methods
            F = scatteredInterpolant(DT,D(:,q),InterpMethod); %Construct interpolator
            DI(:,q)=F(XI); %Interpolate
    end
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 M = spm_mesh_reduce(M,t)
% Reduce the number of triangles in a mesh
% FORMAT M = spm_mesh_reduce(M,f)
% M        - a patch structure
% t        - desired number of triangles
%
% M        - reduced patch structure
%__________________________________________________________________________
%
% References:
%
% M. Garland and P. Heckbert. Surface Simplification Using Quadric Error
% Metrics. In Proceedings of SIGGRAPH 97.
% http://mgarland.org/files/papers/quadrics.pdf
%
% M. Garland and P. Heckbert. Simplifying Surfaces with Color and Texture
% using Quadric Error Metrics. In Proceedings of IEEE Visualization 98.
% http://mgarland.org/files/papers/quadric2.pdf
%
% Wrapper around an implementation by Sven Forstmann, MIT licence:
% https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_reduce.m 7421 2018-09-20 10:58:01Z guillaume $

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

M = reducepatch(M,t);

author: spm


text:
%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% 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.
%============================================================================

%source: http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=05291740
% A Triaxial Accelerometer Calibration Method Using a Mathematical Model
% by S.P. Won, F. Golnaraghi

function [B, G] = IMUcalibIteration(S_x, S_y, S_z, B, G)
    % gravitation at Helsinki
    g0 = 9.8189;

    % current estimate of accelerations
    A_x = (S_x - B(1)) / G(1);
    A_y = (S_y - B(2)) / G(2);
    A_z = (S_z - B(3)) / G(3);

    % squares
    A_x2 = (A_x.*A_x);
    A_y2 = (A_y.*A_y);
    A_z2 = (A_z.*A_z);

    % error estimate
    E = A_x2 + A_y2 + A_z2 - (g0*g0);

    ACCEL = [A_x2 A_y2 A_z2 A_x A_y A_z];

    CAL = ACCEL \ E;

    % gain change
    G_change2 = 1 ./ (1 - CAL(1:3));
    G_change = sqrt(abs(G_change2));

    % bias change
    B_change = CAL(4:6) .* G .* G_change2 .* 0.5;

    % update estimates
    G = G .* G_change;
    B = B + B_change;
end

author: hhyyti

repo: dcm-imu

save_path: github-repos/MATLAB/hhyyti-dcm-imu


text:
function kern = linearKernParamInit(kern)

% LINEARKERNPARAMINIT Linear kernel parameter initialisation.

% KERN

% KERN

kern.variance = 1;
kern.nParams = 1;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function prepare = spm_cfg_eeg_prepare
% Configuration file for the prepare tool
%__________________________________________________________________________
% Copyright (C) 2012-2016 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_prepare.m 7169 2017-09-19 10:42:27Z vladimir $

rev = '$Rev: 7169 $';

D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the M/EEG mat file.'};

defaulttype = cfg_const;
defaulttype.tag = 'defaulttype';
defaulttype.name = 'Set channel types to default';
defaulttype.val  = {1};
defaulttype.help = {'Reset all channel types to default.'};

bidschantype = cfg_files;
bidschantype.tag = 'bidschantype';
bidschantype.name = 'Set channel types from BIDS tsv';
bidschantype.filter = '.*_channels.tsv$';
bidschantype.num = [1 1];
bidschantype.help = {'Select BIDS tsv file and set channel types based on it.'};

newtype = cfg_menu;
newtype.tag = 'newtype';
newtype.name = 'New channel type';
newtype.labels = {'EEG', 'EOG', 'ECG', 'EMG', 'LFP', 'PHYS', 'ILAM', 'Other'};
newtype.val = {'Other'};
newtype.values = {'EEG', 'EOG', 'ECG', 'EMG', 'LFP', 'PHYS', 'ILAM', 'Other'};
newtype.help = {'Select the new channel type to set.'};

settype = cfg_branch;
settype.tag = 'settype';
settype.name = 'Set channel type';
settype.val = {spm_cfg_eeg_channel_selector, newtype};
settype.help = {'Select the new channel type to set.'};

rawmeg = cfg_files;
rawmeg.tag = 'rawmeg';
rawmeg.name = 'Select MEG dataset';
rawmeg.filter = 'any';
rawmeg.num = [1 1];
rawmeg.help = {'Select the MEG dataset to copy sensors from.'};

loadmegsens = cfg_branch;
loadmegsens.tag = 'loadmegsens';
loadmegsens.name = 'Load MEG sensors';
loadmegsens.val = {rawmeg};
loadmegsens.help = {'Reload MEG sensor representation from raw dataset.'};

megfid = cfg_files;
megfid.tag = 'megfid';
megfid.name = 'Select MEG headshape file';
megfid.filter = 'any';
megfid.num = [1 1];
megfid.help = {'Select the file to read MEG headshape from.'};

fidname = cfg_entry;
fidname.tag = 'fidname';
fidname.name = 'MEG fiducial label';
fidname.strtype = 's';
fidname.help = {'Label of a fiducial point as specified in the MEG dataset.'};

hsname = cfg_entry;
hsname.tag = 'hsname';
hsname.name = 'Surface point label';
hsname.strtype = 's';
hsname.help = {'Label of a fiducial point as specified in the surface points file.'};

matching = cfg_branch;
matching.tag = 'matching';
matching.name = 'Matching pair';
matching.help = {'Specify a matching pair of labeled point in the dataset and surface points.'};
matching.val  = {fidname, hsname};

fiducials = cfg_repeat;
fiducials.tag = 'fiducials';
fiducials.name = 'Fiducials';
fiducials.help = {'Specify at least 3 matching pairs of fiducials to coregister the surface points to the sensors.'};
fiducials.num  = [3 Inf];
fiducials.values  = {matching};
fiducials.val = {matching matching matching};

headshape = cfg_branch;
headshape.tag = 'headshape';
headshape.name = 'Load MEG fiducials\headshape';
headshape.val = {megfid, fiducials};
headshape.help = {'Load MEG fiducials or headshape from a file.'};

nasfid = cfg_entry;
nasfid.tag = 'nasfid';
nasfid.name = 'Nasion';
nasfid.strtype = 's';
nasfid.val = {'nas'};
nasfid.help = {'Enter the nasion fiducial label.'};

lpafid = cfg_entry;
lpafid.tag = 'lpafid';
lpafid.name = 'Left';
lpafid.strtype = 's';
lpafid.val = {'lpa'};
lpafid.help = {'Enter the left fiducial label.'};

rpafid = cfg_entry;
rpafid.tag = 'rpafid';
rpafid.name = 'Right';
rpafid.strtype = 's';
rpafid.val = {'rpa'};
rpafid.help = {'Enter the right fiducial label.'};

multimodal = cfg_branch;
multimodal.tag = 'multimodal';
multimodal.name = 'Specify MEG fiducial labels';
multimodal.val  = {nasfid, lpafid, rpafid};
multimodal.help = {'ONLY for multimodal datasets the labels of MEG fiducials should be specified.'};

defaulteegsens = cfg_branch;
defaulteegsens.tag = 'defaulteegsens';
defaulteegsens.name = 'Assign default EEG sensors';
defaulteegsens.val  = {multimodal};
defaulteegsens.help = {'Set EEG sensor locations to SPM template.'};

eegsens = cfg_files;
eegsens.tag = 'eegsens';
eegsens.name = 'Select EEG sensors file';
eegsens.filter = 'any';
eegsens.num = [1 1];
eegsens.help = {'Select the file with EEG electrode coordinates (e.g. Polhemus or SFP).'};

nomatch      = cfg_const;
nomatch.tag  = 'nomatch';
nomatch.name = 'Not necessary';
nomatch.val  = {1};
nomatch.help = {''};

megmatch = cfg_choice;
megmatch.name = 'Match fiducials to MEG';
megmatch.tag = 'megmatch';
megmatch.values = {nomatch, fiducials};
megmatch.val = {nomatch};
megmatch.help = {'Match EEG fiducials to MEG fiducials (only for multimodal datasets).'};

loadeegsens = cfg_branch;
loadeegsens.tag = 'loadeegsens';
loadeegsens.name = 'Load EEG sensors';
loadeegsens.val = {eegsens, megmatch};
loadeegsens.help = {'Load EEG electrode locations.'};

refsens = cfg_branch;
refsens.tag = 'refsens';
refsens.name = 'Select reference sensors';
refsens.val = {spm_cfg_eeg_channel_selector('notype')};
refsens.help = {'Select the sensors to which the EEG recording was referenced',...
    '(select ''All'' for average reference).'};

montage = cfg_files;
montage.tag = 'montage';
montage.name = 'Select montage file';
montage.filter = 'mat';
montage.num = [1 1];
montage.help = {'Select montage file that specifies the referencing.'};

seteegref = cfg_choice;
seteegref.tag = 'seteegref';
seteegref.name = 'Define EEG referencing';
seteegref.values = {refsens, montage};
seteegref.val = {refsens};
seteegref.help = {'Define the way EEG channels were derived from sensors.'};

project3dEEG = cfg_const;
project3dEEG.tag = 'project3dEEG';
project3dEEG.name = 'Project EEG sensors to 2D';
project3dEEG.val  = {1};
project3dEEG.help = {'Project EEG sensor locations to 2D.'};

project3dMEG = cfg_const;
project3dMEG.tag = 'project3dMEG';
project3dMEG.name = 'Project MEG sensors to 2D';
project3dMEG.val  = {1};
project3dMEG.help = {'Project MEG sensor locations to 2D.'};

loadtemplate = cfg_files;
loadtemplate.tag = 'loadtemplate';
loadtemplate.name = 'Load channel template file';
loadtemplate.filter = 'mat';
loadtemplate.num = [1 1];
loadtemplate.help = {'Specify 2D channel locations by loading a template file.'};

status = cfg_menu;
status.tag = 'status';
status.name = 'Status to set';
status.labels = {'GOOD', 'BAD'};
status.val = {1};
status.values = {0, 1};
status.help = {'Select the new channel type to set.'};

setbadchan = cfg_branch;
setbadchan.tag = 'setbadchan';
setbadchan.name = 'Set/unset bad channels';
setbadchan.val = {spm_cfg_eeg_channel_selector, status};
setbadchan.help = {'Set or clear bad flag for channels.'};

bidschanstatus = cfg_files;
bidschanstatus.tag = 'bidschanstatus';
bidschanstatus.name = 'Set bad channels from BIDS tsv';
bidschanstatus.filter = '.*_channels.tsv$';
bidschanstatus.num = [1 1];
bidschanstatus.help = {'Select BIDS tsv file and set channel status based on it.'};

fname = cfg_entry;
fname.tag = 'fname';
fname.name = 'Output file name';
fname.strtype = 's';
fname.val = {'avref_montage.mat'};
fname.help = {''};

avref = cfg_branch;
avref.tag = 'avref';
avref.name = 'Create average reference montage';
avref.help = {'Create an average reference montage for the specific dataset, ',...
    'taking into account bad channels.'};
avref.val  = {fname};

condlabel = cfg_entry;
condlabel.tag = 'label';
condlabel.name = 'Condition label';
condlabel.strtype = 's';
condlabel.help = {''};

condlist = cfg_repeat;
condlist.tag = 'condlist';
condlist.name = 'Specify conditions list';
condlist.num = [1 Inf];
condlist.values = {condlabel};
condlist.help = {'Specify the conditions order manually as a list.'};

condfile = cfg_files;
condfile.tag = 'condfile';
condfile.name = 'Read from conditions file';
condfile.filter = 'mat';
condfile.num = [1 1];
condfile.help = {'Select a mat-file with cell array of condition labels.'};

sortconditions = cfg_choice;
sortconditions.name = 'Sort conditions';
sortconditions.tag = 'sortconditions';
sortconditions.values = {condlist, condfile};
sortconditions.val = {condlist};
sortconditions.help = {'Specify conditions order.'};

eventfile = cfg_files;
eventfile.tag = 'eventfile';
eventfile.name = 'BIDS tsv file';
eventfile.filter = '.*_events.tsv$';
eventfile.num = [1 1];
eventfile.help = {'File to load events from'};

replace = cfg_menu;
replace.tag = 'replace';
replace.name = 'Replace existing events';
replace.labels = {'Replace', 'Add'};
replace.val = {1};
replace.values = {1, 0};
replace.help = {'Replace exisitng events or add to them'};

bidsevents = cfg_branch;
bidsevents.tag = 'bidsevents';
bidsevents.name = 'Load events from BIDS tsv file';
bidsevents.help = {'Load events from BIDS tsv file.'};
bidsevents.val  = {eventfile, replace};

task = cfg_repeat;
task.tag = 'task';
task.name = 'Select task(s)';
task.num = [1 Inf];
task.values = {defaulttype, bidschantype, settype, loadmegsens, headshape,...
    defaulteegsens, loadeegsens, seteegref, project3dEEG, project3dMEG,...
    loadtemplate, setbadchan, bidschanstatus, avref, sortconditions, bidsevents};
task.help = {'Select task(s).'};

prepare = cfg_exbranch;
prepare.tag = 'prepare';
prepare.name = 'Prepare';
prepare.val = {D, task};
prepare.help = {'Prepare EEG/MEG data.'};
prepare.prog = @eeg_prepare;
prepare.vout = @vout_eeg_prepare;
prepare.modality = {'EEG'};

function out = eeg_prepare(job)

D = spm_eeg_load(job.D{1});
for i = 1:numel(job.task)
    S = [];
    S.D = D;
    switch  char(fieldnames(job.task{i}))
        case 'defaulttype'
            S.task = 'defaulttype';
            D = spm_eeg_prep(S);
        case 'bidschantype'
            S.task = 'bidschantype';
            S.filename = char(job.task{i}.bidschantype);
            D = spm_eeg_prep(S);
        case 'settype'
            S.task = 'settype';
            S.type = job.task{i}.settype.newtype;
            S.ind  = D.selectchannels(spm_cfg_eeg_channel_selector(job.task{i}.settype.channels));
            D = spm_eeg_prep(S);
        case 'loadmegsens'
            S.task = 'loadmegsens';
            S.source = char(job.task{i}.loadmegsens.rawmeg{1});
            D = spm_eeg_prep(S);
        case 'headshape'
            S.task = 'headshape';
            S.source = 'convert';
            S.source = char(job.task{i}.loadmegsens.megfid{1});
            S.regfid = {};
            for j = 1:numel(job.task{i}.headshape.matching)
                S.regfid{j, 1} = job.task{i}.headshape.matching(j).hsname;
                S.regfid{j, 2} = job.task{i}.headshape.matching(j).fidname;
            end
            D = spm_eeg_prep(S);
        case 'defaulteegsens'
            S.task = 'defaulteegsens';
            S.regfid ={
                job.task{i}.defaulteegsens.multimodal.nasfid 'spmnas'
                job.task{i}.defaulteegsens.multimodal.lpafid 'spmlpa'
                job.task{i}.defaulteegsens.multimodal.lpafid 'spmrpa'
                };
            D = spm_eeg_prep(S);
        case 'loadeegsens'
            S.task = 'loadeegsens';
            S.source = 'locfile';

            S.sensfile = char(job.task{i}.loadeegsens.eegsens);
            if isfield(job.task{i}.loadeegsens.megmatch, 'fiducials')
                for j = 1:numel(job.task{i}.loadeegsens.megmatch.matching)
                    S.regfid{j, 1} = job.task{i}.loadeegsens.megmatch.matching(j).hsname;
                    S.regfid{j, 2} = job.task{i}.loadeegsens.megmatch.matching(j).fidname;
                end
            end
            D = spm_eeg_prep(S);
        case 'seteegref'
            S.task = 'sens2chan';
            if isfield(job.task{i}.seteegref, 'refsens')
                S.refelec = spm_cfg_eeg_channel_selector(job.task{i}.seteegref.refsens.channels);
            else
                S.montage = char(job.task{i}.seteegref.montage);
            end
            D = spm_eeg_prep(S);
        case 'project3dEEG'
            S.task = 'project3D';
            S.modality = 'EEG';
            D = spm_eeg_prep(S);
        case 'project3dMEG'
            S.task = 'project3D';
            S.modality = 'MEG';
            D = spm_eeg_prep(S);
        case 'loadtemplate'
            S.task = 'loadtemplate';
            S.P = char(job.task{i}.loadtemplate);
            D = spm_eeg_prep(S);
        case 'setbadchan'
            S.task = 'setbadchan';
            S.channels = spm_cfg_eeg_channel_selector(job.task{i}.setbadchan.channels);
            S.status   =  job.task{i}.setbadchan.status;
            D = spm_eeg_prep(S);
        case 'bidschanstatus'
            S.task = 'bidschanstatus';
            S.filename = char(job.task{i}.bidschanstatus);
            D = spm_eeg_prep(S);
        case 'avref'
            eegchan  = D.indchantype('EEG');
            goodind  = D.indchantype('EEG', 'GOOD');

            goodind = find(ismember(eegchan, goodind));

            tra              =  eye(length(eegchan));
            tra(: ,goodind)  =  tra(:, goodind) - 1/length(goodind);

            montage          = [];
            montage.labelorg = D.chanlabels(eegchan);
            montage.labelnew = D.chanlabels(eegchan);
            montage.tra      = tra;

            [p, f]           = fileparts(job.task{i}.avref.fname);
            if isempty(p)
                p = D.path;
            end

            out.avrefname    = {fullfile(p, [f '.mat'])};

            save(char(out.avrefname), 'montage', spm_get_defaults('mat.format'));
        case 'sortconditions'
            S.task = 'sortconditions';
            if isfield(job.task{i}.sortconditions, 'label')
                S.condlist = job.task{i}.sortconditions.label;
            else
                S.condlist = char(job.task{i}.sortconditions.condfile);
            end
            D = spm_eeg_prep(S);
        case 'bidsevents'
            S.task = 'loadbidsevents';
            S.filename = char(job.task{i}.bidsevents.eventfile);
            S.replace  = job.task{i}.bidsevents.replace;

            D = spm_eeg_prep(S);
    end
end

save(D);

out.D = D;
out.Dfname = {fullfile(D)};

function dep = vout_eeg_prepare(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Prepared M/EEG Data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Prepared Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

if ismember('avref', cellfun(@fieldnames, job.task))
    dep(3) = cfg_dep;
    dep(3).sname = 'Average reference montage';
    dep(3).src_output = substruct('.','avrefname');
    % this can be entered into any file selector
    dep(3).tgt_spec   = cfg_findspec({{'filter','mat'}});
end


text:
function [S,t,f,R,Serr]=mtspecgramtrigpb(data,E,win,movingwin,params,fscorr)
% Multi-taper event triggered time-frequency spectrum - binned point process
%
% Usage:
%
% [S,t,f,R,Serr]=mtspecgramtrigpb(data,E,win,movingwin,params,fscorr)
% Input:
%       data        (single channel data vector) -- required
%       E           (event times) - required
%       win         (in the form [winl,winr] i.e window around each event
%                                                 required
%       movingwin         (in the form [window,winstep] i.e length of moving
%                                                 window and step size) -
%                                                 required
%                                                 Note that units for the windows have
%                                                 to be consistent with
%                                                 units of Fs
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%                   Note that T has to be equal to movingwin(1).
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over events when 1, don't average when 0) - optional. Default 0
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
% Output:
%       S       (triggered time-frequency spectrum in form time x frequency x events if segave=0;
%               or in the form time x frequency segave=1)
%       t       (times)
%       f       (frequencies)
%       R       (spike rate)
%       Serr    (error bars) - only for err(1)>=1

if nargin < 4; error('Need data, events and parameters for the windows'); end;
if nargin < 5; params=[]; end;

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear tapers pad fpass err trialave
if nargin < 6 || isempty(fscorr); fscorr=0; end;
data=change_row_to_column(data);
data=createdatamatpb(data,E,Fs,win);
if nargout==5;[S,t,f,R,Serr]=mtspecgrampb(data,movingwin,params,fscorr);
else [S,t,f,R]=mtspecgrampb(data,movingwin,params,fscorr);end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function plotCI(data, m1, s1, ci1, m2, s2, ci2)

% plotCI(data, m1, s1, mci1, m2, s2, mci2)
%   Plot profiles, average profiles, and confidence intervals
%
% Input:
%   data - historical energy usage data
%   m1 - mean daily profile
%   s1 - standard deviation of daily profile
%   ci1 - 95% confidence interval of daily profile
%   m2 - mean month profile
%   s2 - standard deviation of month profile
%   ci2 - 95% confidence interval of month profile

% Copyright 2006-2009 The MathWorks, Inc.

%% Get number of days and number of hours
  numDays  = size(data, 1);
  numHours = size(data, 2);

%% Create Figure
  figure('Units', 'Pixels', 'Position', [100, 100, 800, 400]);

  % Create subplot for visualizing daily profile
  subplot(1,2,1);
  plot(data', 'Color', [.8 .8 .8]);
  patch([1:numHours,numHours:-1:1],[ci1(1,:), ci1(2,end:-1:1)], [.8 .8 .8])
  line(1:numHours, m1, 'Linewidth', 2, 'Color', 'r');
  grid on; xlim([1 numHours]);
  xlabel('hours');
  ylabel('system load (MW)');
  title('Daily Profile (mean & CI)');

  % Create subplot for visualizing month profile
  subplot(1,2,2);
  plot(data, 'Color', [.8 .8 .8]);
  patch([1:numDays,numDays:-1:1],[ci2(1,:), ci2(2,end:-1:1)], [.8 .8 .8])
  line(1:numDays, m2, 'Linewidth', 2, 'Color', 'r');
  grid on;xlim([1 numDays]);
  xlabel('days');
  ylabel('system load (MW)');
  title('Month Profile (mean & CI)');

author: Sable

repo: mcbench-benchmarks


text:
%DemoMV.m
% OVERVIEW, Main script, runs the eval_mvm_twa script which calls the ComputeMVM and ComputeTWA in the Tools/MVM and Tools/TWA subfolders respectively.
%
% The expected output for this example script is given below.
% The results for MVM analysis are stored in the MVMResult structure and TWA analysis are stored in the TWAResult structure.
% The fields of each structure with expected values are listed below.
%
%   Structure.Field                 Expected value
%   MVMResult.energyinband_array    6.60e-07
%   MVMResult.sqi_array             0.9975
%   MVMResult.heart_rate_est_arr    80.11
%
%   TWAResult.HR            [80.32,80.21,79.79,80.65,80.65,80.43,80,80,80.21,80.32,80,79.47]
%   TWAResult.VAlt          [38.84,40.10,37.81,37.38,37.70,38.14,39.32,39.74,38.31,40.39,39.07,42.20]
%   TWAResult.VAlt_Sig      [38.84,40.10,37.81,37.38,37.70,38.14,39.32,39.74,38.31,40.39,39.07,42.20]
%   TWAResult.Noise_Median  [4.28,3.92,4.40,3.75,3.68,4.14,3.00,3.21,4.15,3.95,3.78,7.20]
%   TWAResult.Noise_95      [11.19,10.75,9.61,12.67,11.69,9.72,9.15,9.84,10.29,11.03,9.89,19.40]
%   TWAResult.VAltPt        [132,116,115,125,116,122,130,117,91,111,101,159]
%   TWAResult.successful    1
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

% Add path to the readheader function
addpath(genpath('../../PhysioNet-Cardiovascular-Signal-Toolbox-master/'))
% Load data
sig = load('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twam');
siginfo = readheader('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twam.hea');
fs = siginfo.freq;
ecg = (sig.val - siginfo.adczero)./siginfo.gain; % Adjust signal according to gain and dc offset
clear sig siginfo
if (size(ecg,2) > size(ecg,1))
   ecg = ecg';
end

% Load annotations
ann.QRSon = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','qrson')'; ann.Q = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','q')'; ann.R = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','r')'; ann.S = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','s')'; ann.QRSoff = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','qrsoff')'; ann.Toff = read_ann('../Tools/ECG_Analysis_Tools/MV/Demos/testdata/twa/sample_ecg_twa','toff')';
% Compute MVM and TWA in sample ECG.
%ann.Q = [];
[MVMResult,TWAResult] = eval_mvm_twa(ecg,ann,fs);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function m = imMoment(img, p, q)
%IMMOMENT  Compute simple moment(s) of an image
%
%   M_PQ = imMoment(IMG, P, Q)
%   P is order for x coord, Q is order for y coord.
%
%   Example
%   % generate image
%   img = zeros([5 6]);
%   img(2:3, 2:5) = 1;
%   % compute total mass of image
%   m00  = imMoment(img, 0, 0);
%   % compute centroid
%   cx  = imMoment(img, 1, 0)/m00;
%   cy  = imMoment(img, 0, 1)/m00;
%
%
%   See also
%     imEquivalentEllipse, imPrincipalAxes
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2008-10-08,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% image dimension
dim = size(img);
Dx = dim(2);
Dy = dim(1);

% compute x and y for each pixel
Ix = repmat((1:Dx), [Dy 1]);
Iy = repmat((1:Dy)', [1 Dx]);

% compute moment
m = zeros(size(p));
for i=1:length(p(:))
    m(i) = sum(Ix(:).^p(i) .* Iy(:).^q(i) .* img(:));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%kkapnmexport 'Flexible PNM export'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kapnmexport.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkapnmexport(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kapnmexport - Flexible PNM export
%
%  DESCRIPTION
% This kroutine exports the input object (specified with the [-i] parameter) in PNM format (presently only binary PPM is supported). The output file (specified with the [-o] parameter) will serve as the final name for the output object if there is only one frame in the input object, if there are more than one frames (i.e. if d and/or t > 1) the output file [-o] will be a basename, and a separate file will be created for each frame.
% If there is mask in the image, masked pixels will be substituted by a color (specified by the parameter [-maskcolor] for exporting.
%
%
%
%  EXAMPLES
% All examples for the Annotate toolbox are listed on the $ANNOTATE/manual/index.html on-line manual.
%
%  "SEE ALSO"
% kman page for ka_color_decodecolor or on-line manual for information on the color specification schemes.
%
%  RESTRICTIONS
% If there are more than one frames in the input file, the output file name will be used as a basename, and it will not be available for subsequent operators in a cantata workspace.
%
%  REFERENCES
% All references for the Annotate toolbox are listed on the $ANNOTATE/manual/index.html on-line manual.
%
%  COPYRIGHT
% Annotate Toolbox v1.0 Copyright (C) 1997, Rafael Santos.
%

function varargout = kkapnmexport(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkapnmexport(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kapnmexport"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [x] = qdToState(qd)
% Converts state vector for simulation to qd struct used in hardware.
% x is 1 x 13 vector of state variables [pos vel quat omega]
% qd is a struct including the fields pos, vel, euler, and omega

x = zeros(1,13); %initialize dimensions

x(1:3) = qd.pos;
x(4:6) = qd.vel;

Rot = RPYtoRot_ZXY(qd.euler(1), qd.euler(2), qd.euler(3));
quat = RotToQuat(Rot);

x(7:10) = quat;
x(11:13) = qd.omega;

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
classdef My_copy_layer < dagnn.ElementWise
  properties
  end

  methods
    function outputs = forward(obj, inputs, params)
        outputs = inputs ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
        derInputs = derOutputs ;
        derParams = {} ;
    end

    % ---------------------------------------------------------------------
    function obj = DropOut(varargin)
      obj.load(varargin{:}) ;
    end

    function obj = reset(obj)
      reset@dagnn.ElementWise(obj) ;
    end

  end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [E, E1, E2, Et, Et1, Et2] = testClassSeperability(N,MEAN1,STD1,MEAN2,STD2, PLOT)

%
% This function generates two random Gaussian sequences and computes
% both the theoretical error and the histogram-based method, using
% respectively the functions computeHistError() and theoreticalError().
%
%
% ARGUMENTS:
% N: number of samples per class
% MEAN1: average value of the samples (1st class)
% STD1: standard deviation of the samples (1st class)
% MEAN2: average value of the samples (2nd class)
% STD2: standard deviation of the samples (2nd class)
% PLOT: this is 1 if results are to be plotted
%
% --------------------------------------------
% Theodoros Giannakopoulos
% Dep. of Informatics and Telecommunications
% University of Athens, Greece
% http://www.di.uoa.gr/~tyiannak
% --------------------------------------------
%

if (N<100)
    fprintf('N has to be at lest 1000!\n');
    return;
end

% GENERATE DATASETS (normal distribution)
x1 = randn(N,1) * sqrt(STD1) + MEAN1;
x2 = randn(N,1) * sqrt(STD2) + MEAN2;

if (N<1000)
    nBins = 10;
else
    if (N<10000)
        nBins = 50;
    else
        nBins = 100;
    end
end

% CALCULATE HISTOGRAMS FOR BOTH CLASSES:
step = (max(MEAN1,MEAN2) - min(MEAN1,MEAN2)) / nBins;
X = min(MEAN1-3*STD1,MEAN2-3*STD2): step : max(MEAN1+3*STD1,MEAN2+3*STD2);
[H1, X1] = hist(x1, X);
[H2, X2] = hist(x2, X);

% COMPUTE THE ERROR PROBABILITIES FOR BOTH CLASSES using the historam method:
%
% N O T E : this function computes the error classification probability for
% the training data using the histograms as a pdf estimation method.
% Therefore, it can be used for any distribution of input data.
%

[E1, E2] = computeHistError(x1, x2);
E1 = 100* E1;
E2 = 100* E2;

% OVERALL ERROR PROBABILITY:
E = (E1 + E2) / 2;

if (PLOT==1)
    figure;
    subplot(2,1,1);
    str = sprintf('Estimated Errors: E1 = %.2f%%, E2 = %.2f%%, E = %.2f%%.', E1, E2, E);
    plot(X1, H1);
    hold on;
    plot(X2, H2,'r');
    legend('Class A','Class B');
    title(str);
    axis([min(MEAN1-3*STD1,MEAN2-3*STD2) max(MEAN1+3*STD1,MEAN2+3*STD2) 0 max([H1 H2])]);
    subplot(2,1,2);
    [Et, Et1, Et2, x] = theoreticalError(MEAN1, STD1, MEAN2, STD2);
    str = sprintf('Theoretical Errors: E1 = %.2f%%, E2 = %.2f%%, E = %.2f%%.', Et1, Et2, Et);
    title(str);
end


text:
x=@(alpha,beta) 4*cos(alpha);
y=@(alpha,beta) (5+4*sin(alpha)).*cos(beta);
z=@(alpha,beta) (5+4*sin(alpha)).*sin(beta);
ezsurf(x,y,z)

author: Eurus-Holmes


text:
function m = t_idx( i1, i2, N )
if i1 > i2
    tmp = i1;
    i1 = i2;
    i2 = tmp;
end
m = (i1-1)*N - i1*(i1-1)/2 - i1 + i2;
end

% return i1*N - i1*(i1+1)/2 - i1 + i2 - 1

repo: pTVreg


text:
function RHSdiv = divergenceTermCylindrical1D(F)
% This function calculates the divergence of a field using its face
% average value and the vector u, which is a face vector
%
% SYNOPSIS:
%       RHSdiv = divergenceTermCylindrical1D(F)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = F.domain.dims(1);
G = 1:Nr+2;
DX = F.domain.cellsize.x(2:end-1);
rp = F.domain.cellcenters.x;
rf = F.domain.facecenters.x;

% define the vector of cell index
row_index = reshape(G(2:Nr+1),Nr,1); % main diagonal (only internal cells)

% calculate the flux vector
% note: size(Fx) = [1:m+1]
Fx = F.xvalue;

% reassign the east, west, north, and south flux vectors for the
% code readability
Fe = Fx(2:Nr+1);		Fw = Fx(1:Nr);
re = rf(2:Nr+1);     rw = rf(1:Nr);

% compute the divergence
div_x = (re.*Fe - rw.*Fw)./(rp.*DX);

% define the RHS Vector
RHSdiv = zeros(Nr+2,1);

% assign the values of the RHS vector
RHSdiv(row_index) = reshape(div_x,Nr,1);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
classdef PTKAirways < PTKPlugin
    % PTKAirways. Plugin for segmenting the pulmonary airways from CT data
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAirways calls the PTKTopOfTrachea plugin to find the trachea
    %     location, and then runs the library routine
    %     PTKAirwayRegionGrowingWithExplosionControl to obtain the
    %     airway segmentation. The results are stored in a heirarchical tree
    %     structure.
    %
    %     The output image generated by GenerateImageFromResults creates a
    %     colour-coded segmentation image with true airway points shown as blue
    %     and explosion points shown in red.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Airways'
        ToolTip = 'Shows a segmentation of the airways illustrating deleted points'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Version = 2
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            trachea_results = dataset.GetResult('PTKTopOfTrachea');

            if dataset.IsGasMRI
                threshold = dataset.GetResult('PTKThresholdGasMRIAirways');
                coronal_mode = false;
            elseif strcmp(dataset.GetImageInfo.Modality, 'MR')
                lung_threshold = dataset.GetResult('PTKMRILungThreshold');
                threshold = lung_threshold.LungMask;
                threshold_raw = threshold.RawImage;
                se = ones(1, 3, 10);
                threshold_raw_c = imopen(threshold_raw, se);
                threshold.ChangeRawImage(threshold_raw_c);
                coronal_mode = true;

            else
                threshold = dataset.GetResult('PTKThresholdLungInterior');
                coronal_mode = false;
            end

            % We use results from the trachea finding to remove holes in the
            % trachea, which can otherwise cause early branching of the airway
            % algorithm
            threshold.SetIndexedVoxelsToThis(trachea_results.trachea_voxels, true);

            start_point = trachea_results.top_of_trachea;

            if PTKSoftwareInfo.FastMode
                maximum_number_of_generations = 10;
            else
                maximum_number_of_generations = 15;
            end

            explosion_multiplier = 4;

            debug_mode = PTKSoftwareInfo.GraphicalDebugMode;
            results = PTKAirwayRegionGrowingWithExplosionControl(threshold, start_point, maximum_number_of_generations, explosion_multiplier, coronal_mode, reporting, debug_mode);
        end

        function results = GenerateImageFromResults(airway_results, image_templates, reporting)
            template_image = image_templates.GetTemplateImage(PTKContext.LungROI);
            results = PTKGetImageFromAirwayResults(airway_results.AirwayTree, template_image, false, reporting);
            results_raw = results.RawImage;
            explosion_points = template_image.GlobalToLocalIndices(airway_results.ExplosionPoints);
            results_raw(explosion_points) = 3;
            results.ChangeRawImage(results_raw);
       end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function swave=smoothwavelet(wave,dt,period,dj,scale)
% Smoothing as in the appendix of Torrence and Webster "Inter decadal changes in the ENSO-Monsoon System" 1998
%
% used in wavelet coherence calculations
%
%
% Only applicable for the Morlet wavelet.
%
% (C) Aslak Grinsted 2002-2014
% http://www.glaciology.net/wavelet-coherence

% -------------------------------------------------------------------------
%The MIT License (MIT)
%
%Copyright (c) 2014 Aslak Grinsted
%
%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.
%---------------------------------------------------------------------------

% TODO: take mother argument
%

n=size(wave,2);

%swave=zeros(size(wave));
twave=zeros(size(wave));

% %filter in time:....
% for i=1:size(wave,1)
%     sc=period(i)/dt; % time/cycle / time/sample = samples/cycle
%     t=(-round(sc*3):round(sc*3))*dt;
%     f=exp(-t.^2/(2*scale(i)^2));
%     f=f/sum(f); %filter must have unit weight
%
%     smooth=conv(wave(i,:),f); %slowest line of code in the wtcsig calculation. should be done like in wavelet with fft and ifft.
%     cutlen=(length(t)-1)*.5;
%     twave(i,:)=smooth((cutlen+1):(end-cutlen)); %remove paddings
% end
%
%filter in time:....
%
% qwave=twave;

%zero-pad to power of 2... Speeds up fft calcs if n is large
npad=2.^ceil(log2(n));

k = 1:fix(npad/2);
k = k.*((2.*pi)/npad);
k = [0., k, -k(fix((npad-1)/2):-1:1)];

k2=k.^2;
snorm=scale./dt;
for ii=1:size(wave,1)
    F=exp(-.5*(snorm(ii)^2)*k2); %Thanks to Bing Si for finding a bug here.
    smooth=ifft(F.*fft(wave(ii,:),npad));
    twave(ii,:)=smooth(1:n);
end

if isreal(wave)
    twave=real(twave); %-------hack-----------
end

%scale smoothing (boxcar with width of .6)

%
% TODO: optimize. Because this is done many times in the monte carlo run.
%

dj0=0.6;
dj0steps=dj0/(dj*2);
% for ii=1:size(twave,1)
%     number=0;
%     for l=1:size(twave,1);
%         if ((abs(ii-l)+.5)<=dj0steps)
%             number=number+1;
%             swave(ii,:)=swave(ii,:)+twave(l,:);
%         elseif ((abs(ii-l)+.5)<=(dj0steps+1))
%             fraction=mod(dj0steps,1);
%             number=number+fraction;
%             swave(ii,:)=swave(ii,:)+twave(l,:)*fraction;
%         end
%     end
%     swave(ii,:)=swave(ii,:)/number;
% end

kernel=[mod(dj0steps,1); ones(2 * round(dj0steps)-1,1); ...
   mod(dj0steps,1)]./(2*round(dj0steps)-1+2*mod(dj0steps,1));
swave=conv2(twave,kernel,'same'); %thanks for optimization by Uwe Graichen

%swave=twave;

author: grinsted

repo: wavelet-coherence


text:
% +IOSR
%
%   Contents file for +IOSR and its subfolders.
%
%   +IOSR
%   iosr.install                         - Set search paths, and download and install dependencies
%
%   +IOSR/+ACOUSTICS
%   iosr.acoustics.irStats               - Calculate RT, DRR, Cte, and EDT for impulse response file
%   iosr.acoustics.rtEst                 - Estimate reverberation time based on room size and absorption
%
%   +IOSR/+AUDITORY
%   iosr.auditory.azimuth2itd            - Convert azimuth in degrees to ITD
%   iosr.auditory.binSearch              - Conduct a binary search
%   iosr.auditory.calcIld                - Calculate normalised interaural level difference
%   iosr.auditory.chXcorr                - Calculate cross-correlograms with a wide range of options
%   iosr.auditory.chXcorr2               - Calculate cross-correlograms with a range of options
%   chXcorr2_c.c
%   chXcorr_c.c
%   iosr.auditory.createWindow           - Create a Hann or exp. window with specified onsets/offsets
%   iosr.auditory.dupWeight              - Calculate duplex weighting coefficients for ITD and ILD
%   iosr.auditory.erbRate2hz             - Convert ERB rate to Hz
%   iosr.auditory.freqMulti              - Calculate frequency coefficient for ITD-azimuth warping
%   iosr.auditory.gammatoneFast          - Produce an array of responses from gammatone filters via FFT
%   iosr.auditory.hz2erbRate             - Convert Hz to ERB rate
%   iosr.auditory.iso226                 - ISO 226:2003 Normal equal-loudness-level contours
%   iosr.auditory.itd2azimuth            - Convert ITD to azimuth
%   iosr.auditory.lindemannInh           - Signal pre-processing for Lindemann's cross-correlation
%   iosr.auditory.loudWeight             - Calculate loudness weighting coefficients
%   iosr.auditory.makeErbCFs             - Make a series of center frequencies equally spaced in ERB-rate
%   iosr.auditory.meddisHairCell         - Calculate Ray Meddis' hair cell model for a number of channels
%   iosr.auditory.perceptualCentroid     - Perceptual spectral centroid
%   iosr.auditory.perceptualCentroid2    - Alternative perceptual spectral centroid
%   iosr.auditory.xcorrLindemann         - Cross-correlation based on Lindemann's precedence model
%   xcorrLindemann_c.c
%
%   +IOSR/+BSS
%   iosr.bss.applyIdealMasks             - Calculate and apply ideal masks via STFT
%   iosr.bss.applyMask                   - Apply a time-frequency mask to an STFT
%   iosr.bss.calcImr                     - Calculates the Ideal Mask Ratio (IMR)
%   iosr.bss.calcSnr                     - Calculate the separation SNR
%   iosr.bss.cfs2fcs                     - Calculate gammatone crossover frequencies
%   iosr.bss.example                     - Determine STFT parameters
%   iosr.bss.generateMixtures            - Generate arrays of mixtures from targets and interferers
%   iosr.bss.getFullMask                 - Convert frame rate mask to a sample-by-sample mask
%   iosr.bss.idealMasks                  - Calculate ideal time-frequency masks from STFTs
%   iosr.bss.mixture                     - Class of sound source separation mixture
%   iosr.bss.resynthesise                - Resynthesise a target from a time-frequency mask
%   iosr.bss.source                      - Class of sound source separation source
%
%   +IOSR/+DSP
%   iosr.dsp.audio                       - Abstract superclass providing audio-related properties and methods
%   iosr.dsp.autocorr                    - Perform autocorrelation via FFT
%   iosr.dsp.convFft                     - Convolve two vectors using FFT multiplication
%   iosr.dsp.istft                       - Calculate the Inverse Short-Time Fourier Transform
%   iosr.dsp.lapwin                      - Laplace window
%   iosr.dsp.localpeaks                  - Find local peaks and troughs in a vector
%   iosr.dsp.ltas                        - Calculate the long-term average spectrum of a signal
%   iosr.dsp.matchEQ                     - Match the LTAS of a signal to an arbitrary spectral magnitude
%   iosr.dsp.rcoswin                     - Raised cosine window
%   iosr.dsp.rms                         - Calculate the rms of a vector or matrix
%   iosr.dsp.sincFilter                  - Apply a near-ideal low-pass or band-pass brickwall filter
%   iosr.dsp.smoothSpectrum              - Apply 1/N-octave smoothing to a frequency spectrum
%   iosr.dsp.stft                        - Calculate the short-time Fourier transform of a signal
%   iosr.dsp.vsmooth                     - Smooth a vector using mathematical functions
%
%   +IOSR/+FIGURES
%   iosr.figures.chMap                   - Create a monochrome-compatible colour map
%   iosr.figures.cmrMap                  - Create a monochrome-compatible colour map
%   iosr.figures.multiwaveplot           - Stacked line plots from a matrix or vectors
%   iosr.figures.subfigrid               - Create axis positions for subfigures
%
%   +IOSR/+GENERAL
%   iosr.general.cell2csv                - Output a cell array to a CSV file
%   iosr.general.checkMexCompiled        - Check if mex file is compiled for system
%   iosr.general.getContents             - Get the contents of a specified directory
%   iosr.general.updateContents          - Create a Contents.m file including subdirectories
%   iosr.general.urn                     - Generate random number sequence without duplicates
%
%   +IOSR/+STATISTICS
%   iosr.statistics.boxPlot              - Draw a box plot
%   iosr.statistics.functionalBoxPlot    - Draw a functional boxplot
%   iosr.statistics.functionalPlot       - Abstract superclass for functional plots
%   iosr.statistics.functionalSpreadPlot - Draw a functional plot showing data spread
%   iosr.statistics.getRmse              - Calculate the root-mean-square error between input data
%   iosr.statistics.laprnd               - Pseudorandom numbers drawn from the Laplace distribution
%   iosr.statistics.qqPlot               - Quantile-quantile plot with patch option
%   iosr.statistics.quantile             - Quantiles of a sample via various methods
%   iosr.statistics.statsPlot            - An abstract superclass for classes that plot statistics
%   iosr.statistics.tab2box              - Prepare tabular data for boxPlot function
%   iosr.statistics.trirnd               - Pseudorandom numbers drawn from the triangular distribution
%
%   +IOSR/+SVN
%   iosr.svn.buildSvnProfile             - Read data from files tagged with SVN keywords
%   iosr.svn.headRev                     - Retrieve the head revision for specified files
%   iosr.svn.readSvnKeyword              - Read data from a file tagged with an SVN keyword
%
%   This file was generated by updateContents.m on 31 May 2017 at 17:06:32.

author: IoSR-Surrey

repo: MatlabToolbox


text:
function [x,f,funEvals] = minConF_BC(funObj,x,LB,UB,options)
% function [x,f] = minConF_BC(funObj,x,LB,UB,options)
%
% Function for using Two-Metric Projection to solve problems of the form:
%   min funObj(x)
%   s.t. LB_i <= x_i <= UB_i
%
%   @funObj(x): function to minimize (returns gradient as second argument)
%
%   options:
%       verbose: level of verbosity (0: no output, 1: final, 2: iter (default), 3:
%       debug)
%       optTol: tolerance used to check for progress (default: 1e-7)
%       maxIter: maximum number of calls to funObj (default: 250)
%       numDiff: compute derivatives numerically (0: use user-supplied
%       derivatives (default), 1: use finite differences, 2: use complex
%       differentials)
%       method: 'sd', 'lbfgs', 'newton'

nVars = length(x);

% Set Parameters
if nargin < 5
    options = [];
end
[verbose,numDiff,optTol,maxIter,suffDec,interp,method,corrections,damped] = ...
    myProcessOptions(...
    options,'verbose',3,'numDiff',0,'optTol',1e-6,'maxIter',500,'suffDec',1e-4,...
    'interp',1,'method','lbfgs','corrections',100,'damped',0);

% Output Log
if verbose >= 2
    fprintf('%10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
end

% Make objective function (if using numerical derivatives)
funEvalMultiplier = 1;
if numDiff
    if numDiff == 2
        useComplex = 1;
    else
        useComplex = 0;
    end
    funObj = @(x)autoGrad(x,useComplex,funObj);
    funEvalMultiplier = nVars+1-useComplex;
end

% Evaluate Initial Point
x = projectBounds(x,LB,UB);
if strcmp(method,'newton')
    [f,g,H] = funObj(x);
    secondOrder = 1;
else
    [f,g] = funObj(x);
    secondOrder = 0;
end
funEvals = 1;
if verbose >= 1
  fprintf('Initial function value: %15.5e\n', f);
end

% Compute Working Set
working = ones(nVars,1);
working((x < LB+optTol*2) & g >= 0) = 0;
working((x > UB-optTol*2) & g <= 0) = 0;
working = find(working);

% Check Optimality
if isempty(working)
    if verbose >= 1
        fprintf('All variables are at their bound and no further progress is possible at initial point\n');
    end
    return;
elseif norm(g(working)) <= optTol
    if verbose >=1
        fprintf('All working variables satisfy optimality condition at initial point\n');
    end
    return;
end

if verbose >= 3
    switch method
        case 'sd'
            fprintf('Steepest Descent\n');
        case 'lbfgs'
            fprintf('L-BFGS\n');
        case 'bfgs'
            fprintf('BFGS\n');
        case 'newton'
            fprintf('Newton\n');
    end
end

i = 1;
while funEvals <= maxIter

    % Compute Step Direction
    d = zeros(nVars,1);
    switch(method)
        case 'sd'
            d(working) = -g(working);
        case 'lbfgs'
            if i == 1
                d(working) = -g(working);
                old_dirs = zeros(nVars,0);
                old_stps = zeros(nVars,0);
                Hdiag = 1;
            else
                if damped
                    [old_dirs,old_stps,Hdiag] = dampedUpdate(g-g_old,x-x_old,corrections,verbose==3,old_dirs,old_stps,Hdiag);
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,x-x_old,corrections,verbose==3,old_dirs,old_stps,Hdiag);
                end
                    curvSat = sum(old_dirs(working,:).*old_stps(working,:)) > 1e-10;
               d(working) = lbfgsC(-g(working),old_dirs(working,curvSat),old_stps(working,curvSat),Hdiag);
            end
            g_old = g;
            x_old = x;
        case 'bfgs'
            if i == 1
                d(working) = -g(working);
                B = eye(nVars);
            else
                y = g-g_old;
                s = x-x_old;

                ys = y'*s;

                if i == 2
                    if ys > 1e-10
                        B = ((y'*y)/(y'*s))*eye(nVars);
                    end
                end
                if ys > 1e-10
                    B = B + (y*y')/(y'*s) - (B*s*s'*B)/(s'*B*s);
                else
                    if verbose == 2
                        fprintf('Skipping Update\n');
                    end
                end
                d(working) = -B(working,working)\g(working);
            end
            g_old = g;
            x_old = x;

        case 'newton'
            [R,posDef] = chol(H(working,working));

            if posDef == 0
                d(working) = -R\(R'\g(working));
            else
                if verbose == 3
                    fprintf('Adjusting Hessian\n');
                end
                H(working,working) = H(working,working) + eye(length(working)) * max(0,1e-12 - min(real(eig(H(working,working)))));
                d(working) = -H(working,working)\g(working);
            end
        otherwise
            fprintf('Unrecognized Method: %s\n',method);
            break;
    end

    % Check that Progress can be made along the direction
    f_old = f;
    gtd = g'*d;
    if gtd > -optTol
        if verbose >= 2
            fprintf('Directional Derivative below optTol\n');
        end
        break;
    end

    % Select Initial Guess to step length
    if i == 1 && ~secondOrder
        t = min(1,1/sum(abs(g(working))));
    else
        t = 1;
    end

    % Evaluate the Objective and Projected Gradient at the Initial Step
    x_new = projectBounds(x+t*d,LB,UB);
    if secondOrder
        [f_new,g_new,H] = funObj(x_new);
    else
        [f_new,g_new] = funObj(x_new);
    end
    funEvals = funEvals+1;

    % Backtracking Line Search
    lineSearchIters = 1;
    while f_new > f + suffDec*g'*(x_new-x) || ~isLegal(f_new)
        temp = t;
        if interp == 0 || ~isLegal(f_new) || ~isLegal(g_new)
            if verbose == 3
                fprintf('Halving Step Size\n');
            end
            t = .5*t;
        else
            if verbose == 3
                fprintf('Cubic Backtracking\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d]);
        end

        % Adjust if change is too small
        if t < temp*1e-3
            if verbose == 3
                fprintf('Interpolated value too small, Adjusting\n');
            end
            t = temp*1e-3;
        elseif t > temp*0.6
            if verbose == 3
                fprintf('Interpolated value too large, Adjusting\n');
            end
            t = temp*0.6;
        end

        % Check whether step has become too small
        if sum(abs(t*d)) < optTol
            if verbose == 3
                fprintf('Line Search failed\n');
            end
            t = 0;
            f_new = f;
            g_new = g;
            break;
        end

        % Evaluate New Point
        x_new = projectBounds(x+t*d,LB,UB);
        [f_new,g_new] = funObj(x_new);
        funEvals = funEvals+1;
        lineSearchIters = lineSearchIters+1;

    end

    % Take Step
    x = x_new;
    f = f_new;
    g = g_new;

    % Compute Working Set
    working = ones(nVars,1);
    working((x < LB+optTol*2) & g >= 0) = 0;
    working((x > UB-optTol*2) & g <= 0) = 0;
    working = find(working);

    % Output Log
    if verbose >= 2
        fprintf('%10d %10d %15.5e %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g(working))),norm(g(working)));
    end

    % Check Optimality
    if isempty(working)
        if verbose >= 1
            fprintf('All variables are at their bound and no further progress is possible\n');
        end
        break;
    elseif norm(g(working)) <= optTol
        if verbose >=1
            fprintf('All working variables satisfy optimality condition\n');
        end
        break;
    end

    % Check for lack of progress
    if sum(abs(t*d)) < optTol
        if verbose >= 1
            fprintf('Step size below optTol\n');
        end
        break;
    end

    if abs(f-f_old) < optTol
        if verbose >= 1
            fprintf('Function value changing by less than optTol\n');
        end
        break;
    end

    if funEvals*funEvalMultiplier > maxIter
        if verbose >= 1
            fprintf('Function Evaluations exceeds maxIter\n');
        end
        break;
    end

    % If necessary, compute Hessian
    if secondOrder && lineSearchIters > 1
        [f_new,g_new,H] = funObj(x);
    end

    i = i + 1;
end
end

function [x] = projectBounds(x,LB,UB)
x(x < LB) = LB(x < LB);
x(x > UB) = UB(x > UB);
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [R, G, B] = lab2rgb(L, a, b)
%LAB2RGB Convert an image from CIELAB to RGB
%
% function [R, G, B] = Lab2RGB(L, a, b)
% function [R, G, B] = Lab2RGB(I)
% function I = Lab2RGB(...)
%
% Lab2RGB takes L, a, and b double matrices, or an M x N x 3 double
% image, and returns an image in the RGB color space.  Values for L are in
% the range [0,100] while a* and b* are roughly in the range [-110,110].
% If 3 outputs are specified, the values will be returned as doubles in the
% range [0,1], otherwise the values will be uint8s in the range [0,255].
%
% This transform is based on ITU-R Recommendation BT.709 using the D65
% white point reference. The error in transforming RGB -> Lab -> RGB is
% approximately 10^-5.
%
% See also RGB2LAB.

% By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.
% Fixed a bug in conversion back to uint8 9 September 1999.
% Updated for MATLAB 7 30 March 2009.

if nargin == 1
  b = L(:,:,3);
  a = L(:,:,2);
  L = L(:,:,1);
end

% Thresholds
T1 = 0.008856;
T2 = 0.206893;

[M, N] = size(L);
s = M * N;
L = reshape(L, 1, s);
a = reshape(a, 1, s);
b = reshape(b, 1, s);

% Compute Y
fY = ((L + 16) / 116) .^ 3;
YT = uint8(fY > T1);
fY = uint8(~YT) .* (L / 903.3) + YT .* fY;
Y = fY;

% Alter fY slightly for further calculations
fY = uint8(double(YT) .* (double(fY) .^ (1/3)) + double(~YT) .* (7.787 .* double(fY) + 16/116));

% Compute X
fX = a / 500 + fY;
XT = fX > T2;
X = uint8((double(XT) .* (double(fX) .^ 3) + double(~XT) .* ((double(fX) - 16/116) / 7.787)));

% Compute Z
fZ = fY - b / 200;
ZT = fZ > T2;
Z = uint8((double(ZT) .* (double(fZ) .^ 3) + double(~ZT) .* ((double(fZ) - 16/116) / 7.787)));

% Normalize for D65 white point
X = X * 0.950456;
Z = Z * 1.088754;

% XYZ to RGB
MAT = [ 3.240479 -1.537150 -0.498535;
       -0.969256  1.875992  0.041556;
        0.055648 -0.204043  1.057311];

RGB = max(min(uint8(MAT * double([X; Y; Z])), 1), 0);

R = reshape(RGB(1,:), M, N);
G = reshape(RGB(2,:), M, N);
B = reshape(RGB(3,:), M, N);

if nargout < 2
  R = uint8(round(cat(3,R,G,B) * 255));
end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
%demo_SingleRoomSRIR - Demonstrates the usage of the SingleRoomSRIR conventions.

% #Author: Michael Mihocic
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% SOFA Toolbox - demo script
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Get an empy conventions structure
clear;
conventions='SingleRoomSRIR';
disp(['Creating SOFA file with ' conventions 'conventions...']);
Obj = SOFAgetConventions(conventions);

%% Fill random data...
Obj.Data.IR=rand(4800,1);
Obj.ListenerPosition=zeros(4800,3); Obj.ListenerPosition(:,1)=1;
Obj.SourcePosition=zeros(4800,3); Obj.SourcePosition(:,2)=1;

% Add ReceiverDescriptions as string array
str={};
for ii=1:Obj.API.R
  str{ii,1}=['String' num2str(round(rand(1,1)*10000))];
end
Obj = SOFAaddVariable(Obj,'ReceiverDescriptions','RS',str);

%% Update dimensions
Obj=SOFAupdateDimensions(Obj);

%% Fill with attributes
Obj.GLOBAL_ListenerShortName = 'dummy';
Obj.GLOBAL_History = 'created with a demo script';
Obj.GLOBAL_DatabaseName = 'none';
Obj.GLOBAL_ApplicationName = 'Demo of the SOFA Toolbox';
Obj.GLOBAL_ApplicationVersion = SOFAgetVersion('API');
Obj.GLOBAL_Organization = 'Acoustics Research Institute';
Obj.GLOBAL_AuthorContact = 'michael.mihocic@oeaw.ac.at';

%% save the SOFA file
SOFAfn=fullfile(SOFAdbPath,'sofatoolbox_test',[conventions '_' Obj.GLOBAL_SOFAConventionsVersion '.sofa']);
disp(['Saving:  ' SOFAfn]);
Obj=SOFAsave(SOFAfn, Obj);

author: sofacoustics

repo: SOFAtoolbox


text:
function r8_cube_root_test ( )

%*****************************************************************************80
%
%% R8_CUBE_ROOT_TEST tests R8_CUBE_ROOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_CUBE_ROOT_TEST\n' );
  fprintf ( 1, '  R8_CUBE_ROOT computes the cube root of an R8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       X               Y               Y^3\n' );
  fprintf ( 1, '\n' );

  a = -10.0;
  b = +10.0;
  seed = 123456789;

  for i = 1 : 10
    [ x1, seed ] = r8_uniform_ab ( a, b, seed );
    y = r8_cube_root ( x1 );
    x2 = y ^ 3;
    fprintf ( 1, '  %14.6g  %14.6g  %14.6g\n', x1, y, x2 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [gradf, err] = tapas_riddersgradient(f, x, varargin)
% Calculates the gradient of the function f at point x according to Ridders' method:
%
% Ridders, CJF. (1982). Accurate computation of F'(x) and F'(x) F''(x). Advances in Engineering
%     Software, 4(2), 75-6.
%
% INPUT:
%    f             Function handle of a real function of n real variables which are passed as
%                  *one* vector with n elements
%    x             Point at which to differentiate f
%
% OUTPUT:
%    gradf         Gradient of f at x (row vector)
%    err           Error estimates (row vector)
%
% OPTIONS:
%    Optionally, the third argument of the function can be a structure containing further
%    settings for Ridder's method.
%
%    varargin{1}.init_h      Initial finite difference (default: 1)
%    varargin{1}.div         Divisor used to reduce h on each step (default: 1.2)
%    varargin{1}.min_steps   Minimum number of steps in h (default: 3)
%    varargin{1}.max_steps   Maximum number of steps in h (default: 100)
%    varargin{1}.tf          Terminate if last step worse than preceding by a factor of tf
%                            (default: 2)
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is released under the terms of the GNU General Public Licence (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either version 3 or, at your option,
% any later version). For further details, see the file COPYING or <http://www.gnu.org/licenses/>.

    n = length(x);

    % Defaults
    options.init_h     = 1;
    options.div        = 1.2;
    options.min_steps  = 3;
    options.max_steps  = 100;
    options.tf         = 2;
    gradf              = NaN(1,n);
    err                = NaN(1,n);

    % Overrides
    if nargin > 2
        options_ovr = varargin{1};

        if isfield(options_ovr,'init_h')
            options.init_h = options_ovr.init_h;
        end

        if isfield(options_ovr,'div')
            options.div = options_ovr.div;
        end

        if isfield(options_ovr,'min_steps')
            options.min_steps = options_ovr.min_steps;
        end

        if isfield(options_ovr,'max_steps')
            options.max_steps = options_ovr.max_steps;
        end

        if isfield(options_ovr,'tf')
            options.tf = options_ovr.tf;
        end
    end

    % Check if f and x match
    try
        f(x);
    catch err
        error('tapas:hgf:ridders:CannotEvalFun', 'Function cannot be evaluated at differentiation point.');
    end

    % Loop through argument variables
    for i = 1:n

        % Construct filehandle to be passed to riddersdiff
        fxih = @(xi) fxi(f,x,i,xi);

        % Calculate derivative
        [gradf(i), err(i)] = tapas_riddersdiff(fxih,x(i),options);
    end
end

function fxi = fxi(f,x,i,xi)
    xx    = x;
    xx(i) = xi;
    fxi   = f(xx);
end

author: translationalneuromodeling

repo: tapas


text:
function points = intersectCircles(circle1, circle2)
%INTERSECTCIRCLES Intersection points of two circles.
%
%   POINTS = intersectCircles(CIRCLE1, CIRCLE2)
%   Computes the intersetion point of the two circles CIRCLE1 and CIRCLE1.
%   Both circles are given with format: [XC YC R], with (XC,YC) being the
%   coordinates of the center and R being the radius.
%   POINTS is a 2-by-2 array, containing coordinate of an intersection
%   point on each row.
%   In the case of tangent circles, the intersection is returned twice. It
%   can be simplified by using the 'unique' function.
%
%   Example
%     % intersection points of two distant circles
%     c1 = [0  0 10];
%     c2 = [10 0 10];
%     pts = intersectCircles(c1, c2)
%     pts =
%         5   -8.6603
%         5    8.6603
%
%     % intersection points of two tangent circles
%     c1 = [0  0 10];
%     c2 = [20 0 10];
%     pts = intersectCircles(c1, c2)
%     pts =
%         10    0
%         10    0
%     pts2 = unique(pts, 'rows')
%     pts2 =
%         10    0
%
%   References
%   http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
%   http://mathworld.wolfram.com/Circle-CircleIntersection.html
%
%   See also
%   circles2d, intersectLineCircle, radicalAxis
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-01-20, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% adapt sizes of inputs
n1 = size(circle1, 1);
n2 = size(circle2, 1);
if n1 ~= n2
    if n1 > 1 && n2 == 1
        circle2 = repmat(circle2, n1, 1);
    elseif n2 > 1 && n1 == 1
        circle1 = repmat(circle1, n2, 1);
    else
        error('Both input should have same number of rows');
    end
end

% extract center and radius of each circle
center1 = circle1(:, 1:2);
center2 = circle2(:, 1:2);
r1 = circle1(:,3);
r2 = circle2(:,3);

% allocate memory for result
nPoints = length(r1);
points = NaN * ones(2*nPoints, 2);

% distance between circle centers
d12 = distancePoints(center1, center2, 'diag');

% get indices of circle couples with intersections
inds = d12 >= abs(r1 - r2) & d12 <= (r1 + r2);

if sum(inds) == 0
    return;
end

% angle of line from center1 to center2
angle = angle2Points(center1(inds,:), center2(inds,:));

% position of intermediate point, located at the intersection of the
% radical axis with the line joining circle centers
d1m  = d12(inds) / 2 + (r1(inds).^2 - r2(inds).^2) ./ (2 * d12(inds));
tmp = polarPoint(center1(inds, :), d1m, angle);

% distance between intermediate point and each intersection point
h   = sqrt(r1(inds).^2 - d1m.^2);

% indices of valid intersections
inds2 = find(inds)*2;
inds1 = inds2 - 1;

% create intersection points
points(inds1, :) = polarPoint(tmp, h, angle - pi/2);
points(inds2, :) = polarPoint(tmp, h, angle + pi/2);

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
filename='Cantilever_hexahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function r = isreal(a)
%ISREAL       returns 1 if c is real
%
%   r = isreal(a)
%

% written  05/21/09     S.M. Rump
%

  r = isreal(a.t);

repo: OpenMAS


text:
function [ hScatterObj ] = ScatterSparse( mS, varargin )
% ----------------------------------------------------------------------------------------------- %
% [ hScatterObj ] = ScatterSparse( mS, varargin )
% Visualize sparsity pattern of matrix. Comparable to `spy()` with all the
% controls given by `scatter()`.
% Input:
%   - mS                -   Input Sparse Matrix.
%                           Structure: Matrix (numRows x numCols).
%                           Type: 'Single' / 'Double' (Sprase).
%                           Range: (-inf, inf).
%   - varargin          -   Scatter Plot Parameters.
%                           Set of parameters accepeted by `scatter()`.
%                           Structure: NA.
%                           Type: NA.
%                           Range: NA.
% Output:
%   - hScatterObj       -   Scatter Plot Object.
%                           Structure: Scalar.
%                           Type: Handler / Object.
%                           Range: NA.
% References:
%   1.  A
% Remarks:
%   1.  Supports any functionality by `scatter()` plot.
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     19/07/2021  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    mS (:, :) {mustBeNumeric, mustBeReal, mustBeSparse}
end
arguments (Repeating)
    varargin
end

[vI, vJ, ~] = find(mS);
hScatterObj = scatter(vJ, vI, varargin{:});

end

function [ ] = mustBeSparse( mS )

if(~issparse(mS))
    eid = 'mustBeSparse:mSMustBeSparse';
    msg = 'The input matrix must be sparse.';
    throwAsCaller(MException(eid, msg));
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function selectp(in_or_out)
    %  This .m file selects the earthquakes within a polygon
    %  and plots them. Sets "primeCatalog" equal to the catalogue produced after the
    %  general parameter selection. Operates on "storedcat", replaces "primeCatalog"
    %  with new data and makes "primeCatalog" equal to ZG.newcat
    %
    %   operates on main map window
    % plot tags:
    %  'poly_selected_events' : earthquakes in/out of polygon
    %  'mouse_points_overlay' : polygon outline

    ZG=ZmapGlobal.Data;
    echo on
    % ___________________________________________________________
    %  Please use the left mouse button or the cursor to select
    %  the polygon vertexes.
    %
    %  Use the right mouse button to select the final point.
    %_____________________________________________________________
    echo off
    report_this_filefun();
    %zoom off
    ZG.newt2 = [ ];           % reset catalogue variables
    %a=storedcat;              % uses the catalogue with the pre-selected main
    % general parameters
    ZG.newcat = copy(ZG.primeCatalog);

    delete(findobj('Tag','mouse_points_overlay'));
    delete(findobj('Tag','poly_selected_events'));

    % pick polygon points,
    ax=findobj(gcf,'Tag','mainmap_ax');
    [x,y, mouse_points_overlay] = select_polygon(ax);

    if ~exist('in_or_out','var')
        in_or_out = 'inside';
    end
    if isnumeric(ZG.primeCatalog)
        error('old catalog');
    else
        mask = polygon_filter(x,y, ZG.primeCatalog.Longitude, ZG.primeCatalog.Latitude, in_or_out);
        ZG.newt2 = ZG.primeCatalog.subset(mask);

        % Plot of new catalog
        washeld=ishold(ax); ax.NextPlot='add';
        hMap=mainmap();
        hMap.plotOtherEvents(ZG.newt2,0,...
            'Marker','.',...
            'LineStyle','none',...
            'MarkerEdgeColor','g',...
            'MarkerFaceColor','none',...
            ...'Linewidth',1.5,...
            'DisplayName','Selected Events');

        %change the polygon characteristics
        set(mouse_points_overlay,'LineStyle','--','LineWidth',2,'Color',[.5 .5 .5],'Marker','none');
        if ~washeld; ax.NextPlot='replace';end
    end
    % ask for a name for this new catalog
    %ploc = get(0,'PointerLocation');
    prompt = 'Please provide a catalog name:';
    dlgname = 'CatalogSelection';
    numlines = 1;
    defaultans = {'new_catalog'};
    sel_nm = inputdlg('Please provide a catalog name:','Name of Selected Catalog',1, defaultans);
    if ~isempty(sel_nm)
        ZG.newt2.Name = sel_nm{1};
    end
    xy = [x y];

    %save polcor.dat xy -ascii
    [file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.Directories.data, '*.txt'),'Save Polygon ? (yes/cancel)');
    if length(file1) > 1
        if length(file1)>3
            if strcmp(file1(length(file1)-3:length(file1)),'.txt')==0
                file1=[file1 '.txt']
            end
        end
        save([path1 file1],'xy', '-ascii');
    end

    %++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %
    %   The new catalog (ZG.newcat) with points only within the
    %   selected Polygon is created and resets the original
    %   "primeCatalog" .
    %
    %++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    ZG.newcat = ZG.newt2;                   % resets ZG.newcat and ZG.newt2

    ctp=CumTimePlot(ZG.newt2);
    ctp.plot();
end

author: CelsoReyes


text:
function [GLSZM] = getGLSZM(ROIOnly,levels)
% -------------------------------------------------------------------------
% [GLSZM] = getGLSZM(ROIOnly,levels)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the Gray-Level Size Zone Matrix (GLSZM) of the
% region of interest (ROI) of an input volume. The input volume is assumed
% to be isotropically resampled. The zones of different sizes are computed
% using 26-voxel connectivity.
%
% --> This function is compatible with 2D analysis (language not adapted in the text)
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Thibault, G., Fertil, B., Navarro, C., Pereira, S., Cau, P., Levy,
%     N., Mari, J.-L. (2009). Texture Indexes and Gray Level Size Zone
%     Matrix. Application to Cell Nuclei Classification. In Pattern
%     Recognition and Information Processing (PRIP) (pp. 140–145).
% -------------------------------------------------------------------------
% INPUTS:
% - ROIonly: Smallest box containing the ROI, with the imaging data ready
%            for texture analysis computations. Voxels outside the ROI are
%            set to NaNs.
% - levels: Vector containing the quantized gray-levels in the tumor region
%           (or reconstruction levels of quantization).
%
% ** 'ROIonly' and 'levels' should be outputs from 'prepareVolume.m' **
% -------------------------------------------------------------------------
% OUTPUTS:
% - GLSZM: Gray-Level Size Zone Matrix of 'ROIOnly'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2013
% - Revision: May 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% PRELIMINARY
nLevel = length(levels);
if nLevel > 100
    adjust = 10000;
else
    adjust = 1000;
end
levelTemp = max(levels) + 1;
ROIOnly(isnan(ROIOnly)) = levelTemp;
levels = [levels,levelTemp];

% QUANTIZATION EFFECTS CORRECTION
% In case (for example) we initially wanted to have 64 levels, but due to
% quantization, only 60 resulted.
uniqueVect = round(levels*adjust)/adjust;
ROIOnly = round(ROIOnly*adjust)/adjust;
NL = length(levels) - 1;

% INITIALIZATION
nInit = numel(ROIOnly);
GLSZM = zeros(NL,nInit);

% COMPUTATION OF GLSZM
temp = ROIOnly;
for i = 1:NL
    temp(ROIOnly~=uniqueVect(i)) = 0;
    temp(ROIOnly==uniqueVect(i)) = 1;
    connObjects = bwconncomp(temp,26);
    nZone = length(connObjects.PixelIdxList);
    for j = 1:nZone
        col = length(connObjects.PixelIdxList{j});
        GLSZM(i,col) = GLSZM(i,col) + 1;
    end
end

% REMOVE UNECESSARY COLUMNS
stop = find(sum(GLSZM),1,'last');
GLSZM(:,(stop+1):end) = [];

end

author: mvallieres

repo: radiomics


text:
function [AC, nmi_value, error_cnt] = CalcMetrics(label, result)

result = bestMap(label, result);
error_cnt = sum(label ~= result);
AC = length(find(label == result))/length(label);

nmi_value = nmi(label, result);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% STD Standard deviation.
%    For vectors, Y = STD(X) returns the standard deviation.  For matrices,
%    Y is a row vector containing the standard deviation of each column.  For
%    N-D arrays, STD operates along the first non-singleton dimension of X.
%
%    STD normalizes Y by N-1 if N>1, where N is the sample size.  This is
%    the sqrt of an unbiased estimator of the variance of the population
%    from which X is drawn, as long as X consists of independent,
%    identically distributed samples. For N=1, Y is normalized by N.
%
%    Y = STD(X,1) normalizes by N and produces the square root of the second
%    moment of the sample about its mean.  STD(X,0) is the same as STD(X).
%
%    Y = STD(X,W) computes the standard deviation using the weight vector W.
%    W typically contains either counts or inverse variances.  The length of
%    W must equal the length of the dimension over which STD operates, and
%    its elements must be nonnegative.  If X(I) is assumed to have standard
%    deviation proportional to 1/SQRT(W(I)), then Y * SQRT(MEAN(W)/W(I)) is
%    an estimate of the standard deviation of X(I).  In other words, Y *
%    SQRT(MEAN(W)) is an estimate of standard deviation for an observation
%    given weight 1.
%
%    Y = STD(X,W,DIM) takes the standard deviation along the dimension DIM
%    of X.  Pass in 0 for W to use the default normalization by N-1, or 1 to
%    use N.
%
%    The standard deviation is the square root of the variance (VAR).
%
%    STD(...,NANFLAG) specifies how NaN (Not-A-Number) values are treated.
%    The default is 'includenan':
%
%    'includenan' - the standard deviation of a vector containing NaN
%                   values is also NaN.
%    'omitnan'    - elements of X or W containing NaN values are ignored.
%                   If all elements are NaN, the result is NaN.
%
%    Example:
%        X = [4 -2 1; 9 5 7]
%        std(X,0,1)
%        std(X,0,2)
%
%    Class support for inputs X, W:
%       float: double, single
%
%    See also COV, MEAN, VAR, MEDIAN, CORRCOEF.
%
%    Reference page in Doc Center
%       doc std
%
%    Other functions named std
%
%       datetime/std       duration/std    tall/std          ts/std
%       distributed/std    fints/std       timeseries/std
%

author: jmaih


text:
function airy_ai_int_values_test ( )

%*****************************************************************************80
%
%% AIRY_AI_INT_VALUES_TEST demonstrates the use of AIRY_AI_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AIRY_AI_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  AIRY_AI_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the integral of the Airy Ai function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = airy_ai_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end