text:
function [Rchordal, tchordal] = chordal_initialization(measurements)
% function [Rchordal, tchordal] = chordal_initialization(measurements)
%
% This helper function accepts a MATLAB struct containing the raw
% measurements specifying a special Euclidean synchronization problem, and
% returns the corresponding chordal initialization for this problem.
%
% INPUT: A MATLAB struct 'measurements' containing the following fields
% (see eq. (11) in the long-form version of the paper for details):
% edges:  An (mx2)-dimension encoding the edges in the measurement network;
%     edges(k, :) = [i,j] means that the kth measurement is of the
%     relative transform from pose i to pose j.  NB:  This indexing scheme
%     requires that the states x_i are numbered sequentially as
%     x_1, ... x_n.
% R:  An m-dimensional cell array whose kth element is the rotational part
%     of the kth measurement
% t:  An m-dimensional cell array whose kth element is the translational
%     part of the kth measurement
% kappa:  An m-dimensional cell array whose kth element gives the precision
%     of the rotational part of the kth measurement.
% tau:  An m-dimensional cell array whose kth element gives the precision
%     of the translational part of the kth measurement.

% OUTPUT:
%  -Rchordal:  A d x dn block matrix Rchordal containing the estimating
%     orientations.
%  -tchordal [optional]:  A d x n block matrix containing the estimated
%     positions.

% Copyright (C) 2016 by David M. Rosen

d = length(measurements.t{1});
n = max(max(measurements.edges));
m = size(measurements.edges, 1);

if nargout > 1
    [B3, B2, B1] = construct_B_matrices(measurements);
else
    B3 = construct_B_matrices(measurements);
end

% First, estimate the rotations using *only* the rotational observations
Id = eye(d);
Id_vec = Id(:);

% Compute the constant vector cR induced by fixing the first orientation
% estimate to be the identity Id.
cR = B3(:, 1:d^2) * Id_vec;

% Compute an estimate of the remaining rotations by solving the resulting
% least squares problem *without* enforcing the constraint the the
% estimates lie in SO(d)
r2vec = - B3(:, d^2 + 1 : end) \ cR;
rvec = vertcat(Id_vec, r2vec);
R_LS = reshape(rvec, d, d*n);

% Now reproject these estimates onto SO(d)
Rchordal = zeros(d, d*n);
for i = 1:n
    Rchordal(:, d*(i-1) + 1 : d*i) = project_to_SOd(R_LS(:, d*(i-1) + 1 : d*i));
end

if nargout > 1

    % Solve for the translations in terms of the rotations

    % Constant vector induced by fixing R
    cT = B2 * Rchordal(:);

    % Solve for t_2, ... t_n assuming that t_1 = 0.

    t2 = - B1(:, d + 1 : end) \ cT;
    tvec = vertcat(zeros(d, 1), t2);
    tchordal = reshape(tvec, d, n);
end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [newROIcoords,tmpNewROIcoords] = xformROIcoords(ROIcoords,Xform,inputVoxSize,outputVoxSize,sampRate)
%
% newROIcoords = xformROIcoords(ROIcoords,Xform,inputVoxSize,outputVoxSize,[sampRate])
%
% Transforms ROI coords using Xform, supersampling in each
% dimension to accumulate partial volumes, then keeping only
% those voxels with partial volumes above thresh.
%
% ROIcoords: 3xN matrix of coordinates (y,x,z).
% Xform: 4x4 homogeneous transform
% inputVoxSize: 3-vector, size of voxels (mm) in ROIcoords
% outputVoxSize: 3-vector, size of voxels (mm) in newROIcoords
% sampRate: 3-vector, supersampling rate for each dimension
%           default is odd number >= 4x ratio of inputVoxSize/outputVoxSize
%
% newROIcoords: 3xN matrix of (y,x,z)
%
% djh, 8/98.  Modified from ROIcoords/transformROI.m in mrLoadRet-1
% 7/19/02 djh, modified to maintain equal volumes
% 11/02/05 ARW now returns the 'raw' transformed coordinates - useful for
% various applications including interfacing with SPM
if ~isa(ROIcoords, 'double'), ROIcoords = double(ROIcoords); end

if ~exist('sampRate','var')
	sampRate = ceil(inputVoxSize ./ outputVoxSize) .* [4,4,4];
	sampRate = 2*floor(sampRate/2) + 1;
end

% Convert ROI coords to homogenous coordinates, by adding a fourth
% row of 1's, and transform.
%
tmpNewROIcoords = ones(4,size(ROIcoords,2));
tmpNewROIcoords(1:3,:) = ROIcoords;
tmpNewROIcoords = Xform * tmpNewROIcoords;

% Find bounding (min and max) volume.
%
minPos = [min(tmpNewROIcoords(1,:));min(tmpNewROIcoords(2,:));min(tmpNewROIcoords(3,:));1];
maxPos = [max(tmpNewROIcoords(1,:));max(tmpNewROIcoords(2,:));max(tmpNewROIcoords(3,:));1];
minPos = floor(minPos)-[2,2,2,0]';
maxPos = ceil(maxPos)+[2,2,2,0]';
dims = maxPos(1:3)-minPos(1:3)+ones(3,1);

% Initialize accumulator for partial volume calculation, a vector
% of length appropriate to index the bounding volume.
%
accum = zeros(1,prod(dims));

% Calculate offsets that will be added within the loop to do the
% partial voluming.
%
xoffsets=[-.5+1/(2*sampRate(1)):1/sampRate(1):.5-1/(2*sampRate(1))];
yoffsets=[-.5+1/(2*sampRate(2)):1/sampRate(2):.5-1/(2*sampRate(2))];
zoffsets=[-.5+1/(2*sampRate(3)):1/sampRate(3):.5-1/(2*sampRate(3))];
% xoffsets=[0:1/sampRate(1):1-1/sampRate(1)];
% yoffsets=[0:1/sampRate(2):1-1/sampRate(2)];
% zoffsets=[0:1/sampRate(3):1-1/sampRate(3)];

% Divide alpha by prod(sampRate) to get partial volume for the
% supersampled voxels.
%
alpha = repmat(1/prod(sampRate),[1 size(ROIcoords,2)]);

% Set up a mrvWaitbar if needed.
verbose = prefsVerboseCheck;
if verbose > 1
	waitHandle = mrvWaitbar(0,'Transforming ROI coordinates.  Please wait...');
end

% Loop through supersamples, transform them, and accumulate
% partial volume.
%
ROIcoords = double(ROIcoords);
for ioff=1:length(xoffsets)
	xoff=xoffsets(ioff);

	if verbose>1, mrvWaitbar(ioff/length(xoffsets)); end

	for yoff=yoffsets
		for zoff=zoffsets
			% Add offset
			tmpNewROIcoords(1:3,:) = ROIcoords + ...
				repmat([xoff;yoff;zoff],[1,size(ROIcoords,2)]);

			% Transform
			tmpNewROIcoords = Xform * tmpNewROIcoords;

			% Round and subtract minPos
			coords = round(tmpNewROIcoords(1:3,:)) - ...
				repmat(minPos(1:3),[1,size(tmpNewROIcoords,2)]);

			% Convert to indices and remove duplicates
			indices = coords2Indices(coords,dims);

			% Accumulate partial volume.  Need to do it in a loop
			% instead of:
			%    accum(indices) = accum(indices) + alpha;
			% because an index can appear twice in indices and we want
			% to accumulate them both.
			for jj=1:length(indices)
				accum(indices(jj)) = accum(indices(jj)) + alpha(jj);
			end
		end
	end
end

if verbose>1, close(waitHandle); end

% Build newROIcoords
%
[sortedAccum,indices] = sort(accum);
nonZeroSize = length(find(accum > 0));
newROIsize = round(prod(inputVoxSize)*size(ROIcoords,2) / prod(outputVoxSize));
newROIsize = max(1,newROIsize); % we don't want an ROI size of 0. This will cause an error.
newROIsize = min(nonZeroSize,newROIsize);
indices = indices(length(indices)-newROIsize+1:length(indices));
if ~isempty(indices)
  newROIcoords = indices2Coords(indices,dims);
  newROIcoords = newROIcoords + repmat(minPos(1:3),[1,length(indices)]);
else
  newROIcoords = [];
end

return;

%%%%%%%%%%%%%%
% Debug/test %
%%%%%%%%%%%%%%

ROIcoords = [0; 0; 0];
ROIcoords = [1 2 3 4;
	         1 1 1 1;
	         1 1 1 1];
Xform = [1 0 0 0;
	     0 1 0 0;
	     0 0 1 0.5;
	     0 0 0 1];
inputVoxSize = [1,1,1];
outputVoxSize = [1,1,1/2];
xformROIcoords(ROIcoords,Xform,inputVoxSize,outputVoxSize)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function g = estimate_time_constant(y, p, sn, lags, fudge_factor)
%% Estimate noise standard deviation and AR coefficients if they are not present

%% inputs:
%   y: N X T matrix, fluorescence trace
%   p: positive integer, order of AR system
%   sn: scalar, noise standard deviation, estimated if not provided
%   lags: positive integer, number of additional lags where he autocovariance is computed
%    fudge_factor: float (0< fudge_factor <= 1) shrinkage factor to reduce bias

%% outputs
%   g: 1 x p vector, AR coefficient

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% adapted from the MATLAB implemention by Eftychios Pnevmatikakis and the
% Python implementation from Johannes Friedrich

%% References
% Pnevmatikakis E. et.al., Neuron 2016, Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data

%% input arguments
if ~exist('p', 'var') || isempty(p)
    p = 2;
end
if ~exist('sn', 'var') || isempty(sn)
    sn = GetSn(y);
end
if ~exist('lags', 'var') || isempty(lags)
    lags = 5;
end
if ~exist('fudge_factor', 'var') || isempty(fudge_factor)
    fudge_factor = 1;
end

%% estimate time constants
lags = lags + p;
if ~isempty(which('xcov')) %signal processing toolbox
    xc = xcov(y,lags,'biased');
else
    ynormed = (y - mean(y));
    xc = nan(lags + 1, 1);
    for k = 0:lags
        xc(k + 1) = ynormed(1 + k:end)' * ynormed(1:end - k);
    end
    xc = [flipud(xc(2:end)); xc] / numel(y);
end
xc = xc(:);
A = toeplitz(xc(lags+(1:lags)),xc(lags+(1:p))) - sn^2*eye(lags,p);
g = pinv(A)*xc(lags+2:end);

% while max(abs(roots([1,-g(:)']))>1) && p < 3
%     warning('No stable AR(%i) model found. Checking for AR(%i) model \n',p,p+1);
%     p = p + 1;
%     g = estimate_time_constants(y,p,sn,lags);
% end
% if p == 5
%     g = 0;
% end

% re-adjust time constant values
rg = roots([1;-g(:)]);
if ~isreal(rg); rg = real(rg) + .001*randn(size(rg)); end
rg(rg>1) = 0.95 + 0.001*randn(size(rg(rg>1)));
rg(rg<0) = 0.15 + 0.001*randn(size(rg(rg<0)));
pg = poly(fudge_factor*rg);
g = -pg(2:end);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
classdef PTKLungROI < PTKPlugin
    % PTKLungROI. Plugin for finding the lung region of interest.
    %
    %     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.
    %
    %     PTKLungROI runs the library function PTKGetLungROI to find the region of
    %     interest for the lung image.
    %
    %
    %     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 = 'Lungs & Airways<BR>ROI'
        ToolTip = 'Change the context to display the lungs and airways'
        Category = 'Context'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceImage'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
        Version = 3

        MemoryCachePolicy = 'Temporary'
        DiskCachePolicy = 'Permanent'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            if dataset.IsGasMRI
                results = PTKGetLungROIForGasMRI(dataset.GetResult('PTKOriginalImage', PTKContext.OriginalImage), reporting);
            elseif strcmp(dataset.GetImageInfo.Modality, 'MR')
                lung_threshold = dataset.GetResult('PTKMRILungThreshold', PTKContext.OriginalImage);
                results = dataset.GetResult('PTKOriginalImage', PTKContext.OriginalImage);
                results.ResizeToMatch(lung_threshold.LungMask);
            else
                results = PTKGetLungROIForCT(dataset.GetResult('PTKOriginalImage', PTKContext.OriginalImage), reporting);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function computeAllFinalRF_Aerts_Old(pathRF,nBoot,seed)

startpath = pwd;

cd(pathRF), load('training')
nameOutcomes = {'DeathSign'}; nOutcomes = numel(nameOutcomes);
fSetNames = {'CT'}; nFset = numel(fSetNames);

for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = training.outcomes.Death;
    else
        outcome = training.outcomes.(nameOutcomes{o});
    end
    for f = 1:nFset
        text = training.textures.(nameOutcomes{o}).(fSetNames{f}); nText = size(text,2);
        indClinic = training.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f});
        cost = training.cost.(nameOutcomes{o}).(fSetNames{f});
        tableTrain = [text,training.clinical.table(:,indClinic)];
        cat = logical([zeros(1,nText),training.clinical.categories(indClinic)]);
        rng(seed), [RF] = trainRF_table(tableTrain,outcome,cat,nBoot,cost);
        RF = compact(RF); % Compact version
        save(['RF_',fSetNames{f},'_',nameOutcomes{o}],'RF')
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function options = optiset(varargin)
%OPTISET  Create or alter the options for Optimization with OPTI
%
% options = optiset('param1',value1,'param2',value2,...) creates an OPTI
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the OPTI
% default.
%
% options = optiset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = optiset() creates an options structure with all fields set to
% OPTI defaults.
%
% optiset() prints a list of all possible fields and their function.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

% Print out possible values of properties.
if((nargin == 0) && (nargout == 0))
    printfields();
    return
end

%Names and Defaults
Names = {'solver';'maxiter';'maxfeval';'maxnodes';'maxtime';'tolrfun';'tolafun';'tolint';'solverOpts';'dynamicOpts';'iterfun';'warnings';'display';'derivCheck'};
Defaults = {'auto';1500;1e4;1e4;1000;1e-7;1e-7;1e-5;[];[];[];'critical';'off';'off'};

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

%Ensure backwards compatible
if(strcmpi(options.warnings,'on'))
    options.warnings = 'all';
elseif(strcmpi(options.warnings,'off'))
    options.warnings = 'none';
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar non negative double
    case {'tolafun','tolrfun','tolint'}
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar non zero double
    case 'maxtime'
        err = opticheckval.checkScalarGrtZ(value,field);
    %Scalar non zero integer
    case {'maxiter','maxfeval','maxnodes'}
        err = opticheckval.checkScalarIntGrtZ(value,field);
    %Struct
    case {'solveropts','dynamicopts'}
        err = opticheckval.checkStruct(value,field);
    %Function Handle
    case 'iterfun'
        err = opticheckval.checkFunHandle(value,field);
    %Other misc
    case 'display'
        err = opticheckval.checkValidString(value,field,{'off','iter','final'});
    case 'warnings'
        err = opticheckval.checkValidString(value,field,{'on','critical','off','all','none'});
    case 'derivcheck'
        err = opticheckval.checkValidString(value,field,{'on','off'});
    case 'solver'
        err = opticheckval.checkOptiSolver(value,field);

    otherwise
        err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults

solvers = [{'AUTO'} checkSolver('all')];
len = length(solvers);
str = '';
for i = 1:len
    if(i < len)
        str = [str solvers{i} ', ']; %#ok<AGROW>
    else
        str = [str solvers{i}]; %#ok<AGROW>
    end
end
str = regexprep(str,'AUTO','{AUTO}');

fprintf('\n OPTISET Fields:\n');
fprintf(['           solver: [ ' str ' ] \n']);
fprintf('           maxiter: [ Maximum Solver Iterations {1.5e3} ] \n');
fprintf('          maxfeval: [ Maximum Function Evaluations {1e4} ] \n');
fprintf('          maxnodes: [ Maximum Integer Solver Nodes {1e4} ] \n');
fprintf('           maxtime: [ Maximum Solver Evaluation Time {1e3s} ] \n');
fprintf('           tolrfun: [ Relative Function Tolerance {1e-7} ] \n');
fprintf('           tolafun: [ Absolute Function Tolerance {1e-7} ] \n');
fprintf('            tolint: [ Absolute Integer Tolerance {1e-5} ] \n');
fprintf('        solverOpts: [ Solver Specific Options Structure ] \n');
fprintf('       dynamicOpts: [ Dynamic Optimization Options Structure ] \n');
fprintf('           iterfun: [ Iteration Callback Function, stop = iterfun(iter,fval,x) {} ] \n');
fprintf('          warnings: [ ''all'' or {''critical''} or ''none'' ] \n');
fprintf('           display: [ {''off''}, ''iter'', ''final'' ] \n');
fprintf('        derivCheck: [ Derivative Checker {''off''}, ''on'' ] \n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [D, DGC] = plotConcVSdGft0GroupContUncertainty(modelT)
% Compares the difference between minimum & maximum concentration, on a
% logarithmic scale, and the group contribution uncertainty for each
% metabolite.
%
% USAGE:
%
%    [D, DGC] = plotConcVSdGft0GroupContUncertainty(modelT)
%
% INPUT:
%    modelT:    structure with fields:
%
%                 * modelT.concMax
%                 * modelT.concMin
%                 * modelT.dfGt0GroupContUncertainty
%
% OUTPUTS:
%    D:
%    DGC:

[nMet,nRxn]=size(modelT.mets);

RT=modelT.T*modelT.gasConstant;

D=zeros(nMet,1);
DGC=zeros(nMet,1);
q=0;
q2=0;
r=0;
for m=1:nMet
    d=(RT*log(modelT.concMax(m)/modelT.concMin(m)))/2;
    dgc=modelT.dfGt0GroupContUncertainty(m);
    if isempty(d) || isnan(dgc)
        D(m)=NaN;
        DGC(m)=NaN;
    else
        D(m)=d;
        DGC(m)=dgc;
        if dgc>(d/2)
            q=q+1;
        end
        if dgc>(d)
            q2=q2+1;
        end
        r=r+1;
    end
end

fraction=q/r;
fraction2q=q2/r

q2=0;
r2=0;
for m=1:nMet
    d=(RT*log(modelT.concMax(m)/modelT.concMin(m)))/2;
    dgc=modelT.dfGt0GroupContUncertainty(m);
    if strcmp(modelT.dfGt0Source(m),'Keq')
        dgc=NaN;
    end
    if ~(isempty(d) || isnan(dgc))
        if dgc>(d/2)
            q2=q2+1;
        end
        r2=r2+1;
    end
end

fraction2r=q2/r2;

fprintf('%s\n',['Fraction of metabolites where GC uncertainty is more significant: ' num2str(fraction)]);

d=(RT*log(20/0.2))/2;

% figure
% plot(D,DGC,'.')
% plot(DGC,'.')

figure1=figure;
axes1 = axes('Parent',figure1,'FontSize',14,'CLim',[1 2],'Layer','top');
hist(DGC,100)
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
line([d d],[0 300],'Color','b','LineWidth',4,'LineStyle','--')

text('Position',[d,175],'String','$$\longleftarrow\;5.9=\frac{1}{2}RT\ln\left(\frac{20}{0.2}\right)$$','Interpreter','latex','FontSize',18,'Color','b');

xlim([0, 100])
set(gca,'XTick',0:10:100,'TickDir','out')
ylabel('# Reactants','FontSize',14)
xlabel('Standard Error in {\Delta_{f}G_{est}^{0}} (kJ/mol)','FontSize',14)

author: opencobra

repo: cobratoolbox


text:
function prob = brunbirthday ( nw )

% For this script the default configuration should be 'ithaca_2009b

% nw is the number of workers use 4 as the default

  if nargin < 1 || nw > 4 || nw < 1; nw = 4; end

% Since the batch command is inside a function, we wait for the job to
% complete before retrieving the results and setting the return value
% of prob.
% If you want to continue working in Matlab while the job runs, issue the
% my_job = batch(..); from the Command Window. Check for 'finished'
% by using state = get(my_job,'state')

% Create a batch job for the birthday simulation
%
  my_job = batch ( 'birthdayscript', ...
    'FileDependencies', 'birthday.m', ...
    'matlabpool', nw );
%
% birthdayscript.m uses a parfor loop (matlabpool)
% The job requires one worker to run the script and two
% additional workers for a matlabool
%
% Wait for the job to complete
%
  wait ( my_job );
%
%  Retrieve output from the job
%
  results = getAllOutputArguments(my_job);
  prob = results{1}.prob;
%
% If the cell array results is empty (say 1 x 0), the script submitted
% by 'batch' has failed.

% Could alternatively use the load command
% prob = load(my_job, 'prob')

%
% Clean up after the job
%
  destroy(my_job);

author: johannesgerer

repo: jburkardt-m


text:
classdef TestDescriptorMatcher
    %TestDescriptorMatcher

    properties (Constant)
        fields = {'queryIdx', 'trainIdx', 'imgIdx', 'distance'};
    end

    methods (Static)
        function test_1
            X = randn(20,3);
            Y = randn(5,3);
            matchers = { ...
                cv.DescriptorMatcher('BruteForce'), ...
                cv.DescriptorMatcher('BruteForce-L1'), ...
                cv.DescriptorMatcher('BruteForce-SL2'), ...
                cv.DescriptorMatcher('FlannBased'), ...
                cv.DescriptorMatcher('BFMatcher', ...
                    'NormType','L2', 'CrossCheck',false), ...
                cv.DescriptorMatcher('FlannBasedMatcher', ...
                    'Index',  {'KDTree', 'Trees',4}, ...
                    'Search', {'Checks',32, 'EPS',0, 'Sorted',true})
            };
            for i = 1:numel(matchers)
                typename = matchers{i}.typeid();
                matchers{i}.add(X);
                matchers{i}.train();

                m = matchers{i}.match(Y);
                validateattributes(m, {'struct'}, {'vector'});
                assert(all(ismember(TestDescriptorMatcher.fields, ...
                    fieldnames(m))));

                k = 3;
                m = matchers{i}.knnMatch(Y, k);
                validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});
                cellfun(@(mm) validateattributes(mm, {'struct'}, ...
                    {'vector', 'numel',k}), m);
                cellfun(@(mm) assert(all(ismember(...
                    TestDescriptorMatcher.fields, fieldnames(mm)))), m);

                maxDist = 2.0;
                m = matchers{i}.radiusMatch(Y, maxDist);
                validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});
                cellfun(@(mm) validateattributes(mm, ...
                    {'struct'}, {'vector'}), m);
                cellfun(@(mm) assert(all(ismember(...
                    TestDescriptorMatcher.fields, fieldnames(mm)))), m);
                assert(all(cellfun(@(mm) all([mm.distance] <= maxDist), m)));
            end
        end

        function test_2
            X = randi([0,255], [20,3], 'uint8');
            Y = randi([0,255], [5,3], 'uint8');
            matchers = { ...
                cv.DescriptorMatcher('BruteForce-Hamming'), ...
                cv.DescriptorMatcher('BruteForce-HammingLUT'), ...
                ... cv.DescriptorMatcher('BruteForce-Hamming(2)')  %TODO
            };
            for i = 1:numel(matchers)
                typename = matchers{i}.typeid();
                matchers{i}.add(X);
                matchers{i}.train();

                m = matchers{i}.match(Y);
                validateattributes(m, {'struct'}, {'vector'});
                assert(all(ismember(TestDescriptorMatcher.fields, ...
                    fieldnames(m))));

                k = 3;
                m = matchers{i}.knnMatch(Y, k);
                validateattributes(m, {'cell'}, ...
                    {'vector', 'numel',size(Y,1)});
                cellfun(@(mm) validateattributes(mm, {'struct'}, ...
                    {'vector', 'numel',k}), m);
                cellfun(@(mm) assert(all(ismember(...
                    TestDescriptorMatcher.fields, fieldnames(mm)))), m);

                maxDist = 10.0;
                m = matchers{i}.radiusMatch(Y, maxDist);
                validateattributes(m, {'cell'}, ...
                    {'vector', 'numel',size(Y,1)});
                cellfun(@(mm) validateattributes(mm, ...
                    {'struct'}, {'vector'}), m);
                cellfun(@(mm) assert(all(ismember(...
                    TestDescriptorMatcher.fields, fieldnames(mm)))), m);
                assert(all(cellfun(@(mm) all([mm.distance] <= maxDist), m)));
            end
        end

        function test_3
            X = randn(20,3);
            Y = randn(5,3);
            matcher = cv.DescriptorMatcher('BruteForce');
            assert(matcher.isMaskSupported());
            mask = rand(size(Y,1), size(X,1)) > 0.5;
            m = matcher.match(Y, X, 'Mask',mask);
            validateattributes(m, {'struct'}, {'vector'});
            m = matcher.knnMatch(Y, X, 3, 'Mask',mask);
            validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});
            m = matcher.radiusMatch(Y, X, 2.0, 'Mask',mask);
            validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});
        end

        function test_4
            files = dir(fullfile(mexopencv.root(),'test','shape0*.png'));
            if isempty(files)
                error('mexopencv:testskip', 'missing data');
            end
            query = [];

            % compute descriptors for all training images
            obj = cv.ORB();
            matcher = cv.DescriptorMatcher('BruteForce-Hamming');
            idxQuery = 1;
            idxTrain = setdiff(1:numel(files), idxQuery);
            for i=1:numel(files)
                img = cv.imread(fullfile(mexopencv.root(),'test',files(i).name), 'Flags',0);
                [~,desc] = obj.detectAndCompute(img);
                assert(~isempty(desc));
                if i == idxQuery
                    query = desc;
                else
                    matcher.add(desc);
                end
            end

            D = matcher.getTrainDescriptors();
            validateattributes(D, {'cell'}, ...
                {'vector', 'numel',numel(idxTrain)});
            cellfun(@(desc) validateattributes(desc, ...
                {obj.descriptorType()}, ...
                {'size',[NaN obj.descriptorSize()]}), D);

            % find matches between query and train set
            matcher.train();
            m = matcher.match(query);

            % naive approach to find nearest image
            [~,ord] = sort([m.distance]);         % retain best 100 (min dist)
            m = m(ord(1:min(100,end)));
            ii = idxTrain(mode([m.imgIdx]) + 1);  % most frequest image index
            query = files(idxQuery).name;
            closest = files(ii).name;
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function fighandle=dockfig(varargin)
% fighandle=dockfig(figs)      DOCK one or more FIGures.
%   "fighandle" is an array of handle numbers for the figures.
%   "figs" is an array of figure handles.
%       If "figs" is empty then all figures are docked.
%   Examples:
%       figh=dockfig([4,3,5]), dockfig, dockfig 4 3 5
% Note that the command line format is supported.

%   Copyright 2006 Mirtech, Inc.
%   created 09/20/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.       email: mhrovat@email.com
figs=[];
if nargin>0,
    for n=1:length(varargin),
        arg=varargin{n};
        if ischar(arg),
            tmp=str2num(arg);   %#ok
            if isempty(tmp),    % since it is empty then it may be a variable
                tmp = evalin('base',arg);
                if isempty(tmp),
                    error ('  Syntax error in argument!'),
                end
            end  % isempty(tmp)
        else
            tmp=arg;
        end  % ischar(arg)
        figs=[figs,tmp];
    end  % for
end  % if nargin
if isempty(figs),       figs=0;     end

% set this to "on" so that all figures with "HandleVisibility" set to off can be seen.
set(0,'ShowHiddenHandles','on')
% Below are two alternative methods to find all of the figures.
% Second method should be more foolproof.
% allhands=get(0,'Children');
allhands=findobj('Type','figure');
set(0,'ShowHiddenHandles','off')
if figs==0,
    figs=sort(allhands)';
else
    [tmp,indx]=intersect(figs,allhands);    %#ok tmp not used, view figures that exist
    figs=figs(sort(indx));                  % order the figures as specified
end
numfigs=numel(figs);

if isempty(figs),
    warning('MATLAB:EmptyValue','There are no figures to dock.')
else
    for n=numfigs:-1:1,
        set(figs(n),'WindowStyle','docked');
        figure(figs(n))             % bring figure to front.
    end
end

if nargout==1,
    fighandle=figs;
end

author: Sable

repo: mcbench-benchmarks


text:
% rejstatepoch() - reject bad eeg trials based a statistical measure. Can
%                 be applied either to the raw eeg data or the ICA
%                 component activity. This is an interactive function.
%
% Usage:
%   >> [ Irej, Irejdetails, n, threshold, thresholdg] = ...
%                         rejstatepoch( signal, rej, 'key1', value1...);
%
% Inputs:
%   signal     - 3 dimensional data array channel x points x trials
%                (instead of channels, one might also use independent
%                components).
%   rej        - rejection array, one value per trial and per channel or
%                component (size channel x trials).
%                By default, values are normalized by this
%                function and trehshold is expressed in term of standard
%                deviation of the mean.
%
% Optional inputs:
%   'plot'       - ['on'|'off'] interactive mode or just rejection.
%                  In the interactive mode, it plots the normalized
%                  entropy of original signal (blue) and the limits
%                  (red) (default:'on')
%   'threshold'  - percentage error threshold (default 1-0.25/nb_trials)
%                  for individual trials of individual channel/component.
%                  This treshold is expressed in term of standart
%                  deviation from the mean (default 5).
%   'global'     - ['on'|'off'], also perform threshold on the global
%                  measure (by default, the mean over all channel or
%                  electrodes).
%   'rejglob'    - rejection array, one value per trials. Use this
%                  argument when the global measure for all channels or
%                  does not correspond to their mean.
%   'thresholdg' - global threshold for the reunion of all channels
%                  or components. By default, it is equal to 'threshold'.
%   'normalize'  - ['on'|'off'], normalize values before applying the
%                  treshold. Default is 'on'.
%   'plotcom'    - sting command to plot single trials. Default
%                  is none.
%   'title'      - title of the graph. Default is none.
%   'labels'     - labels for electrodes (array not cell).
%
% Outputs:
%   Irej        - indexes of trials to be rejected
%   Irejdetails - array for rejected components or channel (nb_rejected x
%                 nb_channel or nb_rejected x nb_components)
%   n           - number of trials rejected
%   threshold   - percentage error threshold
%   thresholdg  - percentage error threshold for global rejection
%
% See also: eeglab()

% Algorithm:
%   normalise the measure given as input and reject trials based on
%   an uniform distribution of the data

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% userdata
% gcf : plotsig - signal to plot in the pop_out window
%       pnts - number of points per epoch
%       Irej - rejection of trials
%       Irejelec - rejection of trials and electrodes
%       labels - labels for curves
% plotwin : rej1 - rejection array (with all electrodes x trials)
%           rej2 - global rejection array (size trials)
%           thr1 - threshold (electrodes)
%           thr2 - threshold global

function [ Irej, Irejdetails, n, threshold, thresholdg] = rejstatepoch( signal, ...
		rej, varargin); % pnts, th_E, th_rejg, command, commandplot, typetitle, E, rejg);

if nargin < 1
	help rejstatepoch;
	return;
end;

if ~isstr( signal )

	if nargin < 2
		help rejstatepoch;
		return;
	end;

	if ~isempty( varargin ), g=struct(varargin{:});
	else g= []; end;
	try, g.plot; 			catch, g.plot='on'; end;
	try, g.threshold; 		catch, g.threshold=5; end;
	try, g.thresholdg; 		catch, g.thresholdg=5; end;
	try, g.global; 			catch, g.global='on'; end;
	try, g.rejglob; 		catch, g.rejglob=[]; end;
	try, g.normalize; 		catch, g.normalize='on'; end;
	try, g.plotcom; 		catch, g.plotcom=''; end;
	try, g.title;	 		catch, g.title=''; end;
	try, g.labels;	 		catch, g.labels=''; end;

	g.rej = rej;
	clear rej
	switch lower(g.plot)
		case {'on', 'off'} ;
		otherwise disp('Error: Plot must be either ''on'' or ''off'''); return;
	end;
	switch lower(g.global)
		case {'on', 'off'} ;
		otherwise disp('Error: Global must be either ''on'' or ''off'''); return;
	end;
	switch lower(g.normalize)
		case {'on', 'off'} ;
		otherwise disp('Error: Normalize must be either ''on'' or ''off'''); return;
	end;
	if ~isstr(g.plotcom)
		disp('Error: Plotcom must be a string to evaluate'); return;
	end;
	if ~isstr(g.title)
		disp('Error: Title must be a string'); return;
	end;
	try, g.threshold*2;
		catch, disp('Error: Threhsold must be a number'); return;
	end;
	try, g.thresholdg*2;
		catch, disp('Error: Threhsoldg must be a number'); return;
	end;
	if length(g.threshold(:)) > 1
		disp('Error: Threhsold must be a single number'); return;
	end;
	if length(g.thresholdg(:)) > 1
		disp('Error: Threhsoldg must be a single number'); return;
	end;
	if ~isempty(g.rejglob)
		if length(g.rejglob) ~= size(g.rej,2)
			disp('Error: Rejglob must be have the same length as rej columns'); return;
		end;
	else
		switch lower(g.global), case 'on', g.rejglob = sum(g.rej,1); end;
	end;
	if size(signal,3) ~= size(g.rej,2)
		disp('Error: Signal must be have the same number of element in 3rd dimension as rej have columns'); return;
	end;
	if isempty(g.labels)
		for index = 1:size(g.rej,1)
			g.labels(index,:) = sprintf('%3d', index);
		end;
		if ~isempty(g.rejglob)
			g.labels(index+2,:) = 'g. ';
		end;
	end;

	switch lower(g.normalize),
		case 'on',
			g.rej = (g.rej-mean(g.rej,2)*ones(1, size(g.rej,2)))./ (std(g.rej, 0, 2)*ones(1, size(g.rej,2)));
			switch lower(g.global), case 'on', g.rejglob = (g.rejglob(:)-mean(g.rejglob(:)))./ std(g.rejglob(:)); end;
	end;
	switch lower(g.global), case 'off',g.rejglob = []; end;

	% plot the buttons
	% ----------------
	try, icadefs; catch, GUIBUTTONCOLOR = [0.8 0.8 0.8]; BACKCOLOR = [0.8 0.8 0.8]; end;
	figure('color', BACKCOLOR);
	set(gcf, 'name', 'Rejectrials');
	pos = get(gca,'position'); % plot relative to current axes
	set( gca, 'tag', 'mainaxis');
	q = [pos(1) pos(2) 0 0];
	s = [pos(3) pos(4) pos(3) pos(4)]./100; % allow to use normalized position [0 100] for x and y
	axis('off');

	plotsig = sum(abs(signal(:,:)),1);
	set(gcf, 'userdata', { plotsig, size(signal, 2), [], [] }); % the two last arguments are the rejection

	% Create axis
	% -----------
	h6 = axes('Units','Normalized', 'tag',  'Plotwin', 'Position',[-10 12 120 84].*s+q);
	title(g.title);
	set( h6, 'userdata', { g.rej g.rejglob g.threshold g.thresholdg g.labels }); % g.rej was put twice because it is used to compute the global entropy

	% CANCEL button
	% -------------
	h  = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'Cancel', 'Units','Normalized','Position',...
			[-15 -6 15 6].*s+q, 'callback', 'close(gcf);', 'backgroundcolor', GUIBUTTONCOLOR);

	% Entropy component text and threshold
	% ------------------------------------
    makebutton( 'Single-channel', 'Plotwin'  ,  [3 -6 27 6].*s+q, [30 -6 15 6].*s+q, 3, g.threshold, GUIBUTTONCOLOR);
    makebutton( 'All-channel', 'Plotwin'     ,  [50 -6  27 6].*s+q, [77 -6  15 6].*s+q, 4, g.thresholdg, GUIBUTTONCOLOR);

	% EEGPLOT button
	% --------------
	if ~isempty(g.plotcom)
			h  = uicontrol(gcf, 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'pushbutton', 'string', 'EEGPLOT', 'Units','Normalized','Position',[95 -2 15 6].*s+q, ...
					'callback',['TMPEEG = get(gcbf, ''userdata'');' ...
								'TMPREJ = TMPEEG{3};' ...
								'TMPREJELEC = TMPEEG{4};' ...
								g.plotcom ] );
			posaccept = [95 -10 15 6];
	else	posaccept = [95 -6 15 6];
	end;

	% ACCEPT button
	% -------------
	command = 'fprintf(''Rejection indices has been put in the matrix TMPREJ\n'')';
	haccept  = uicontrol(gcf, 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'pushbutton', 'string', 'UPDATE', 'Units','Normalized','Position', posaccept.*s+q, ...
					'callback', [	'set(gcbo, ''userdata'', 1);' ... %signal to signify termination
									'TMPEEG = get(gcbf, ''userdata'');' ...
									'TMPREJ = TMPEEG{3};' ...
									'TMPREJELEC = TMPEEG{4};' ...
									command ] );

	command = [ 'entwin = get(gcbf, ''currentobject'');' ...
				'tmptype = get(entwin,''type'');' ...
				'if tmptype(1) == ''l'' entwin = get(entwin, ''parent''); end;' ...
				'tagwin = get(entwin,''tag'');' ... % either entropy or kurtosis
				'switch tagwin,' ...
				' case ''Plotwin'',' ... % check if the user clicked on the right window
				'   alldata = get(gcf, ''userdata'');' ...
				'   plotsig = alldata{1};' ...
				'   pnts = alldata{2};' ...
				'   fig = figure(''position'', [100 300 600 400],''color'', [1 1 1]);' ...
				'   I = alldata{3};' ...
				'   sweeps = size(plotsig,2) /pnts;' ...
				'   h1 = axes(''parent'', fig, ''Units'',''Normalized'', ''Position'',[0.6 0.11 0.27 0.815]);' ...
				'   pos = get(entwin, ''currentpoint'');' ...
				'   component = round(pos(1) / sweeps + 0.5);' ... % determine the component number
				'   alldata = get(entwin, ''userdata'');' ...
		    	'   rej = alldata{1};' ...
				'   if component <= size(rej,1)' ... % component
				'        rej_threshold = alldata{3};' ...
				'        component = max(1,min(component, size(rej,1))); ' ...
		    	'   	 rej = rej(component, :);' ...
				'        titlegraph = sprintf(''' g.title ' #%d'', component);' ...
				'        colorgraph = ''b'';' ...
				'   else' ...                        % global
		    	'   	 rej = alldata{2};' ...
				'        rej_threshold = alldata{4};' ...
				'        titlegraph = sprintf(''' g.title ' global'');' ...
				'        colorgraph = ''g'';' ...
				'   end;' ...
				'   plot([1:length(rej)], rej, ''color'', colorgraph);' ...
				'   ss = get(h1, ''xlim'');' ...
				'   set(h1, ''view'', [90 90]);' ...
				'   set(h1, ''xdir'', ''reverse'');' ...
				'   set(h1, ''XLim'', ss);' ...
				'   hold on;' ...  % plot component
				'   yl = get(h1, ''ylim'');' ...
				'   set(h1, ''xtickmode'', ''manual'', ''xtick'', sweeps/2, ''xticklabel'', component, ''xlim'', [ 0 sweeps ]);' ...
				'   title( titlegraph );' ...
				'   plot( get(h1, ''xlim''), [rej_threshold rej_threshold], ''r'');' ... % plot limit
				'   plot( get(h1, ''xlim''), [-rej_threshold -rej_threshold], ''r'');' ... % plot limit
				'   set(h1, ''xticklabel'', []);' ...
				'   hold on;' ...
				'   h2 = axes(''parent'', fig,''Units'',''Normalized'', ''Position'',[0.13 0.11 0.27 0.815]);' ...
				'   erpimage( plotsig, ones(1,size(plotsig,2)/pnts), [0:1/(pnts-1):1], '''', 3, 1, ''nosort'', ''noplot'');' ...
				'   title(''Currentset all chans''); xlabel(''''); ylabel(''''); ' ...
				'   set(gca, ''xticklabel'', []);' ...
				'   hold on;' ...
				'   h3 = axes(''parent'', fig,''Units'',''Normalized'', ''Position'',[0.45 0.11 0.1 0.815]);' ...
				'   if any(I ~= 0)' ...
				'      rejImage = (I'' * ones(1, 10))'';' ...
				'      imagesc( rejImage'' );' ...
				'      set(gca, ''ydir'', ''normal'');' ...
				'   end;' ...
				'   title(''Rejected (all)''); xlabel(''''); ylabel('''');' ...
				'   set(gca, ''xticklabel'', [], ''yticklabel'', []);' ...
				'end;' ...
				'clear fig tmptype tagwin alldata rej rejImage plotsig sweeps pnts rej_threshold ss q s h1 h2 h3 pos component yl;' ];

%				'      erpimage( rejImage(:)'', ones(1,size(I,2)), [0:1/(10-1):1], '''', 1, 0, ''nosort'', ''noplot'');' ...
	set(gcf, 'WindowButtonDownFcn', command);

	rejstatepoch('draw');
	switch g.plot,
		case 'on', waitfor( haccept, 'userdata'); drawnow;
	end;

	threshold  = g.threshold;
	thresholdg = g.thresholdg;
	Irej = [];
	Irejdetails = [];
	n = 0;
	try
		TMPEEG = get(gcf, 'userdata');
		Irej = TMPEEG{3};
		Irejdetails = TMPEEG{4};
		n = length(find(Irej == 1));

		plothandler = findobj( 'parent', gcf, 'tag', 'Plotwin');
		TMPEEG = get(plothandler, 'userdata');
		threshold = TMPEEG{3};
		thresholdg = TMPEEG{4};
		close(gcf);
	catch, end;
else %if signal is a string draw everything

	% retreive data
	% -------------
	gcfdata = get(gcf, 'userdata');
	plotsig = gcfdata {1};
	pnts    = gcfdata {2};
    sweeps  = size(plotsig,2)/pnts;

	h6 = findobj('parent', gcf, 'tag', 'Plotwin');
	alldata = get(h6, 'userdata');
	g.rej       = alldata {1};
	g.rejg      = alldata {2};
	g.threshold   = alldata {3};
	g.thresholdg  = alldata {4};
	set(h6, 'userdata', alldata);

	nbchans = size(g.rej,1);

	% reject trials
	% -------------
	rejelec = abs(g.rej) > g.threshold;
	rej  = max(rejelec,[],1);
	n1 = sum(rej(:));
	if ~isempty( g.rejg )
		rej2 = abs(g.rejg) > g.thresholdg;
		n2 = sum(rej2(:));
		rej = rej | rej2(:)';
	end;
	fprintf('%d trials rejected (single:%d, all:%d)\n', sum(rej(:)), n1, n2);
	gcfdata {3} = rej;
	gcfdata {4} = rejelec;
	set(gcf, 'userdata', gcfdata);

	% plot the sorted entropy curve
	% -----------------------------
	plotstat( 'Plotwin');

end;
return;

function plotstat( id );

	% plot the sorted entropy curve
	% -----------------------------
	h6 = findobj('parent', gcf, 'tag', id);
	axes(h6); cla;
	ttmp = get(gca, 'title');
	oldtitle = get(ttmp, 'string');

	% get datas
    % ---------
   	alldata = get(gca, 'userdata');
	g.rej       = alldata {1};
	g.rejg      = alldata {2};
	g.threshold   = alldata {3};
	g.thresholdg  = alldata {4};
	g.labels      = alldata {5};
	nbchans = size(g.rej,1);
	sweeps  = size(g.rej,2);

	% plot datas
    % ----------
	g.rej = g.rej'; plot(g.rej(:)); g.rej = g.rej';
	hold on;
	yl = get(gca, 'ylim');

	% plot vertival bars to separate components and the trehsold
	% ----------------------------------------------------------
	set( gca, 'tag',  id, 'ylimmode', 'manual');
	set(gca, 'xtickmode', 'manual', 'xtick', [0:sweeps:(size(g.rej(:),1)-1+2*sweeps)] + sweeps/2, ...
			 'xticklabel', g.labels, 'xlim', [ 0 (size(g.rej(:),1)-1+2*sweeps)]);
	plot( [1 size(g.rej(:),1)], [-g.threshold -g.threshold], 'r');	% plot threshold
	plot( [1 size(g.rej(:),1)], [g.threshold g.threshold], 'r');	% plot threshold

	if ~isempty(g.rejg) % plot global ?
		plot([size(g.rej(:),1)+sweeps:size(g.rej(:),1)+2*sweeps-1],  g.rejg(:), 'g');
		pp = patch([size(g.rej(:),1) size(g.rej(:),1) size(g.rej(:),1)+sweeps size(g.rej(:),1)+sweeps], [yl(1)-1 yl(2)+1 yl(2)+1 yl(1)-1], get(gcf, 'color'), 'clipping', 'off');
		set(pp, 'EdgeColor',  get(gcf, 'color'));
		plot( [size(g.rej(:),1)+sweeps length(g.rejg)+size(g.rej(:),1)+sweeps], [-g.thresholdg -g.thresholdg], 'r');	% plot threshold
		plot( [size(g.rej(:),1)+sweeps length(g.rejg)+size(g.rej(:),1)+sweeps], [g.thresholdg g.thresholdg], 'r');	% plot threshold
		plot([size(g.rej(:),1)+sweeps size(g.rej(:),1)+sweeps], yl, 'k');
	else
		pp = patch([size(g.rej(:),1) size(g.rej(:),1) size(g.rej(:),1)+2*sweeps size(g.rej(:),1)+2*sweeps], [yl(1)-1 yl(2)+1 yl(2)+1 yl(1)-1], get(gcf, 'color'), 'clipping', 'off');
		set(pp, 'EdgeColor',  get(gcf, 'color'));
	end;
	for index = 0:sweeps:size(g.rej(:),1);
		plot([index index], yl, 'k');
	end;

	% restore properties
	title(oldtitle);
	set(h6, 'userdata', alldata);

return;

function makebutton( string, ax, pos1, pos2, userindex, init, GUIBUTTONCOLOR );
	h  = uicontrol(gcf , 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'radiobutton', 'string', string, 'value', fastif(init == 0, 0, 1), 'Units','Normalized', 'Position', pos1, ...
					'callback', [ 'textresh = findobj(''parent'', gcbf, ''tag'', ''' string ''');' ...
								  'checkstatus = get(gcbo, ''value'');' ...
								  'ax = findobj(''parent'', gcbf, ''tag'', ''' ax ''');' ...
								  'tmpdat = get(ax, ''userdata'');' ...
								  'if checkstatus' ... % change status of the textbox
								  '   set(textresh, ''enable'', ''on'');' ...
								  '   tmpdat{' int2str(userindex) '} = str2num(get(textresh, ''string''));' ...
								  'else' ...
								  '   set(textresh, ''enable'', ''off'');' ...
								  '   tmpdat{' int2str(userindex) '} = 0;' ...
								  'end;' ...
								  'set(ax, ''userdata'' , tmpdat);' ...
								  'rejstatepoch(''draw'');' ...
								  'clear tmpdat ax checkstatus textresh;'  ] );
	h  = uicontrol(gcf, 'Style', 'edit', 'backgroundcolor', [1 1 1], 'tag', string, 'string', num2str(init), 'enable', fastif(init == 0, 'off', 'on'), 'Units','Normalized', 'Position', pos2, ...
					'callback', [ 	'ax = findobj(''parent'', gcbf, ''tag'', ''' ax ''');' ...
								    'tmpdat = get(ax, ''userdata'');' ...
									'tmpdat{' int2str(userindex) '} = str2num(get(gcbo, ''string''));' ...
									'set(ax, ''userdata'' , tmpdat);' ...
								    'rejstatepoch(''draw'');' ...
									'clear tmpdat ax;' ] );
return;

author: PatternRecognition

repo: OpenBMI


text:
function[overlaps] = computeBlobOverlapSum(propBlobsA, propBlobsB, imageSize)
% [overlaps] = computeBlobOverlapSum(propBlobsA, propBlobsB, imageSize)
%
% Take Sel. Search regions and reconstruct which regions overlap by how many pixels.
% (result is a matrix)
%
% Copyright by Holger Caesar, 2015

if numel(imageSize) == 3,
    imageSize = imageSize(1:2);
end;
propCountA = numel(propBlobsA);
propCountB = numel(propBlobsB);
overlaps = zeros(propCountA, propCountB); %Dense is fine

% Precompute blob inds
blobsIndsA = cell(propCountA, 1);
for i = 1 : propCountA,
    blob = propBlobsA(i);
    blobsIndsA{i} = blobToImageInds(blob, imageSize);
end;

blobsIndsB = cell(propCountB, 1);
for i = 1 : propCountB,
    blob = propBlobsB(i);
    blobsIndsB{i} = blobToImageInds(blob, imageSize);
end;

for i = 1 : propCountA,
    for j = 1 : propCountB,
        % Compute the intersection / union
        % (Faster than ismember())
        overlaps(i, j) = sum(builtin('_ismemberhelper', blobsIndsA{i}, blobsIndsB{j}));
%         assert(overlap == sum(ismember(blobsIndsA{i}, blobsIndsB{j})));
%         assert(sum(ismember(blobsIndsA{i}, blobsIndsB{j})) == sum(ismember(blobsIndsA{j}, blobsIndsB{i})));
    end;
end;

author: nightrome

repo: matconvnet-calvin


text:
function [stats, ntmaps, hh, hhfill, table_group, multcomp_group] = hansen_neurotransmitter_maps(fmri_data_obj, varargin)
% Apply a set of PET neurotransmitter-related maps from Hansen 2022 Nature Neuroscience
% - Uses cosine similarity metric on gray matter masked images
%
% :Usage:
% ::
%
%     [stats, ntmaps, hh, hhfill, table_group, multcomp_group] = hansen_neurotransmitter_maps(fmri_data_obj, [optional inputs])
%
% :Inputs:
%
%   **fmri_data_obj:**
%        An fmri_data object with one or more images; can be a signature
%
% :Optional Inputs:
%   **colors:**
%        Followed by plot colors in cell array
%
%   **nofigure:**
%        Do not create a new figure; use existing
%
%   **noplot:**
%        Do not plot polar plot, return tables only
%
%   **correlation:**
%        Use correlation similarity metric instead of cosine similarity
%
%   **dofixrange:**
%        Followed by range of values [min max] for plot (see
%        image_similarity_plot)
%
% :Outputs:
%
%   **stats:**
%        Stats structure with associations
%
%   **ntmaps:**
%        fmri_data object with neurotransmitter maps and meta-data, in
%        sorted order
%
%   **'hh'**
%        line handles
%
%   **'hhfill'**
%        fill handles
%
%    **table_group**
%              multiple one-way ANOVA tables (one for each
%              spatial basis) with group as column factor (requires
%              'average' to be specified)
%
%    **multcomp_group**
%              mutiple comparisons of means across groups, one output
%              cell for each spatial basis, critical value determined
%              by Tukey-Kramer method (see multcompare)
%
% :Examples:
% ::
%
% % Load and plot profile for a negative affect neuromarker from Chang et al. 2015, Plos Biology:
% % ------------------------------------------------------------------
% pines = load_image_set('pines');
% [stats, ntmaps] = hansen_neurotransmitter_maps(pines);
%
% figure;
% [stats, ntmaps] = hansen_neurotransmitter_maps(pines, 'nofigure', 'colors', [0 0 1]);
% [stats, ntmaps] = hansen_neurotransmitter_maps(pines, 'noplot');
%
% % Load a drug craving neuromarker from Koban et al. 2022, Nature
% % Neuroscience, and plot the profiles for drug-only and drug+food signatures on top of
% % one another:
% % ------------------------------------------------------------------
% ncs_all = load_image_set('ncs'); % craving
% ncsdrug = get_wh_image(ncs_all, 2);
% ncs = get_wh_image(ncs_all, 1);
% [stats, ntmaps, hh, fillh, table_group] = hansen_neurotransmitter_maps(ncsdrug);
% [stats, ntmaps, hh, fillh, table_group] = hansen_neurotransmitter_maps(ncs, 'colors', {[0 0 1]}, 'nofigure');

% :References:
%   Hansen et al. 2022, Nature Neuroscience
%
% :See also:
%   - image_similarity_plot
%

%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2023 Mijin Kwon and Tor Wager
%
%     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/>.
% ..

% ..
%    Programmers' notes:
%    Tor: Created from Mijin's script -- plus resort rows of image names and
%    fullpath fields
% ..

% -------------------------------------------------------------------------
% DEFAULT ARGUMENT VALUES
% -------------------------------------------------------------------------

colors = [1 0 0];
dofigure = true;
doplot = true;
similarity_metric = 'cosine_similarity';
dofixrange = [];

% -------------------------------------------------------------------------
% OPTIONAL INPUTS
% -------------------------------------------------------------------------

% This is a compact way to assign multiple variables. The input argument
% names and variable names must match, however:

allowable_inputs = {'colors' 'doplot' 'similarity_metric' 'dofixrange'};

keyword_inputs = {'noplot' 'nofigure' 'correlation'};

% optional inputs with default values - each keyword entered will create a variable of the same name

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case allowable_inputs

                eval([varargin{i} ' = varargin{i+1}; varargin{i+1} = [];']);

            case keyword_inputs
                % Skip deal with this below

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% 2nd pass: Keyword inputs. These supersede earlier inputs
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'noplot'
                doplot = false;

            case 'nofigure'
                dofigure = false;

            case 'correlation'
                similarity_metric = 'correlation';

        end
    end
end

% -------------------------------------------------------------------------
% MAIN FUNCTION
% -------------------------------------------------------------------------

ntmaps = load_image_set('hansen22');

ntmaps = reorder_and_add_metadata(ntmaps);

% These are already gray-matter masked in repo, but make sure:
ntmaps = apply_mask(ntmaps, which('gray_matter_mask.nii'));

% This may not be masked...so mask with gray matter:
fmri_data_obj = apply_mask(fmri_data_obj, which('gray_matter_mask.nii'));

if dofigure
    create_figure('Neurotransmitter polar plot')
end

if doplot
    if ~iscell(colors), colors = {colors}; end

    if isempty(dofixrange)

        [stats, hh, hhfill, table_group, multcomp_group] = image_similarity_plot(fmri_data_obj, 'mapset', ntmaps, similarity_metric, 'plotstyle', 'polar', 'networknames', ntmaps.metadata_table.target, 'colors', colors, 'nofigure');

    else % we have fixed range

        [stats, hh, hhfill, table_group, multcomp_group] = image_similarity_plot(fmri_data_obj, 'mapset', ntmaps, similarity_metric, 'plotstyle', 'polar', 'networknames', ntmaps.metadata_table.target, 'colors', colors, 'nofigure', 'dofixrange', dofixrange);

    end

else
    [stats, hh, hhfill, table_group, multcomp_group] = image_similarity_plot(fmri_data_obj, 'mapset', ntmaps, similarity_metric, 'noplot');

end

end % main function

% ---------------------------------------------------------------
% Subfunctions
% ---------------------------------------------------------------

function ntmaps = reorder_and_add_metadata(ntmaps)

nrows = height(ntmaps.metadata_table);

transmitter = cell(nrows, 1);
ntmaps.metadata_table = addvars(ntmaps.metadata_table, transmitter,'After','target');
ntmaps.metadata_table.Properties.VariableNames{2} = 'transmitter';

for i = 1:nrows

    if ismember(ntmaps.metadata_table.target{i}, {'5HT1a', '5HT1b', '5HT2a', '5HT4', '5HT6', '5HTT'})
        ntmaps.metadata_table.transmitter{i} = 'serotonin';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'CB1'})
        ntmaps.metadata_table.transmitter{i} = 'cannabinoid';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'D1', 'D2', 'DAT'})
        ntmaps.metadata_table.transmitter{i} = 'dopamine';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'GABAa', 'GABAabz'})
        ntmaps.metadata_table.transmitter{i} = 'GABA';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'H3'})
        ntmaps.metadata_table.transmitter{i} = 'histamine';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'M1', 'a4b2', 'VAChT'})
        ntmaps.metadata_table.transmitter{i} = 'acetylcholine';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'MOR'})
        ntmaps.metadata_table.transmitter{i} = 'opioid';
    end
    if ismember(ntmaps.metadata_table.target{i}, {'NET'})
        ntmaps.metadata_table.transmitter{i} = 'norepinephrine';
    end

    if ismember(ntmaps.metadata_table.target{i}, {'mGluR5'})
        ntmaps.metadata_table.transmitter{i} = 'glutamate';
    end
end

% display(ntmaps.metadata_table)

[ntmaps.metadata_table, wh_indx_all] = sortrows(ntmaps.metadata_table,'transmitter','ascend');

ntmaps.dat = ntmaps.dat(:, wh_indx_all);

% Save original image names in sorted order
ntmaps.metadata_table = addvars(ntmaps.metadata_table, ntmaps.image_names(wh_indx_all, :), 'NewVariableName', {'image_names'}, 'Before', 'modeling_notes');

% re-sort
ntmaps.fullpath = ntmaps.fullpath(wh_indx_all, :);

end % end subfunction

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Calibrating kl, ku, mu and nu for a sabr model
% the prices are calculated from known sabr parameters

function y = of_sabr(a, b, r, n, f, k, t, mu, nu, kl, ku, call, put)
% objective function for SABR calibration
    [kl, cl, bl, al, ku, cu, bu, au] = ...
        psabr_param_2(a, b, r, n, f, k, t,mu,nu, kl, ku);

    callkl = interp1(k,call,kl);        % get value at kl
    putku = interp1(k,put,ku);          % get value at ku

    % using only call or put prices is enough
    % can use put left of atm and calls right!!!
    % that is put for [kl,f) and calls (f,ku]
    y = (sprice_5(a, b, r, n, f, kl, t, kl, ku, ...
        mu, cl, bl,al, nu, cu,bu,au, 1) ...
        - callkl)^2;           % call prics
    y = y + (sprice_5(a, b, r, n, f, 0, t, kl, ku, ...
        mu, cl, bl,al, nu, cu,bu,au, 0) ...
        - f)^2;                % atm price
    y = y + (sprice_5(a, b, r, n, f, ku, t, kl, ku, ...
        mu, cl, bl,al, nu, cu,bu,au, 0) ...
        - putku)^2;            % put prices


text:
function err = getHcurlerror3NE(node,elem,curlE,Eh,markedElem)
%% GETHCURLERROR3NE Hcurl norm of approximation error for the lowest order Nedelect element in 3-D.
%
% err = getHcurlerror3NE(node,elem,curlE,Eh,markedElem);
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     maxIt = 4;
%     HcurlErr = zeros(maxIt,1);
%     N = zeros(maxIt,1);
%     for k = 1:maxIt
%         [node,elem] = uniformbisect3(node,elem);
%         [elem2dof,edge] = dof3edge(elem);
%         pde = Maxwelldata2;
%         uI = edgeinterpolate(pde.exactu,node,edge);
%         HcurlErr(k) = getHcurlerror3NE(node,elem,pde.curlu,uI);
%         N(k) = length(uI);
%     end
%     r = showrate(N,HcurlErr,1,'b-+');
%     legend('||u-u_I||_{curl}',['N^{' num2str(r) '}'],'LOCATION','Best');
%
% See also getHcurlerror3NE1, getHcurlerror3NE2, getL2error3NE
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Sort elem to ascend ordering
elem = sort(elem,2);

%% Construct Data Structure
elem2dof = dof3edge(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
% compute curl phi
[Dlambda,volume] = gradbasis3(node,elem);
% curl phi = 2*Dlambda_i cross Dlambda_j;
curlPhi(:,:,6) = 2*cross(Dlambda(:,:,3),Dlambda(:,:,4),2);
curlPhi(:,:,1) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,2),2);
curlPhi(:,:,2) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,3),2);
curlPhi(:,:,3) = 2*cross(Dlambda(:,:,1),Dlambda(:,:,4),2);
curlPhi(:,:,4) = 2*cross(Dlambda(:,:,2),Dlambda(:,:,3),2);
curlPhi(:,:,5) = 2*cross(Dlambda(:,:,2),Dlambda(:,:,4),2);
curlEhp = zeros(NT,3);
for k = 1:6
    curlEhp = curlEhp + ...
    repmat(Eh(elem2dof(:,k)),1,3).*curlPhi(:,:,k);
end

%% compute Hcurl error element-wise
[lambda,w] = quadpts3(2);
nQuad = size(lambda,1);
err = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    if isnumeric(curlE) % a constant vector
        curlEp = repmat(curlE,NT,1);
    else % function handel
        curlEp = curlE(pxy);
    end
%     curlEp = curlE(pxy);
    % compute Ehp at quadrature points
    err = err + w(p)*sum((curlEp - curlEhp).^2,2);
end
err = err.*volume;
% modify the error
err(isnan(err)) = 0; % remove the singular part
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % error on some marked region
end
err = sqrt(sum(err));
%% TODO write more M-lint

repo: Aero-matlab


text:
%
%   Copyright (C) 2016  Starsky Wong <sununs11@gmail.com>
%
%   Note: The SIFT algorithm is patented in the United States and cannot be
%   used in commercial products without a license from the University of
%   British Columbia.  For more information, refer to the file LICENSE
%   that accompanied this distribution.

function [ descrs, locs ] = getFeatures( input_img )
% Function: Get sift features and descriptors
global gauss_pyr;
global dog_pyr;
global init_sigma;
global octvs;
global intvls;
global ddata_array;
global features;
if(size(input_img,3)==3)
    input_img = rgb2gray(input_img);
end
input_img = im2double(input_img);

%% Build DoG Pyramid
% initial sigma
init_sigma = 1.6;
% number of intervals per octave
intvls = 3;
s = intvls;
k = 2^(1/s);
sigma = ones(1,s+3);
sigma(1) = init_sigma;
sigma(2) = init_sigma*sqrt(k*k-1);
for i = 3:s+3
    sigma(i) = sigma(i-1)*k;
end
% default cubic method
input_img = imresize(input_img,2);
% assume the original image has a blur of sigma = 0.5
input_img = gaussian(input_img,sqrt(init_sigma^2-0.5^2*4));
% smallest dimension of top level is about 8 pixels
octvs = floor(log( min(size(input_img)) )/log(2) - 2);

% gaussian pyramid
[img_height,img_width] =  size(input_img);
gauss_pyr = cell(octvs,1);
% set image size
gimg_size = zeros(octvs,2);
gimg_size(1,:) = [img_height,img_width];
for i = 1:octvs
    if (i~=1)
        gimg_size(i,:) = [round(size(gauss_pyr{i-1},1)/2),round(size(gauss_pyr{i-1},2)/2)];
    end
    gauss_pyr{i} = zeros( gimg_size(i,1),gimg_size(i,2),s+3 );
end
for i = 1:octvs
    for j = 1:s+3
        if (i==1 && j==1)
            gauss_pyr{i}(:,:,j) = input_img;
        % downsample for the first image in an octave, from the s+1 image
        % in previous octave.
        elseif (j==1)
            gauss_pyr{i}(:,:,j) = imresize(gauss_pyr{i-1}(:,:,s+1),0.5);
        else
            gauss_pyr{i}(:,:,j) = gaussian(gauss_pyr{i}(:,:,j-1),sigma(j));
        end
    end
end
% dog pyramid
dog_pyr = cell(octvs,1);
for i = 1:octvs
    dog_pyr{i} = zeros(gimg_size(i,1),gimg_size(i,2),s+2);
    for j = 1:s+2
    dog_pyr{i}(:,:,j) = gauss_pyr{i}(:,:,j+1) - gauss_pyr{i}(:,:,j);
    end
end
% for i = 1:size(dog_pyr,1)
%     for j = 1:size(dog_pyr{i},3)
%         imwrite(im2bw(im2uint8(dog_pyr{i}(:,:,j)),0),['dog_pyr\dog_pyr_',num2str(i),num2str(j),'.png']);
%     end
% end

%% Accurate Keypoint Localization
% width of border in which to ignore keypoints
img_border = 5;
% maximum steps of keypoint interpolation
max_interp_steps = 5;
% low threshold on feature contrast
contr_thr = 0.04;
% high threshold on feature ratio of principal curvatures
curv_thr = 10;
prelim_contr_thr = 0.5*contr_thr/intvls;
ddata_array = struct('x',0,'y',0,'octv',0,'intvl',0,'x_hat',[0,0,0],'scl_octv',0);
ddata_index = 1;
for i = 1:octvs
    [height, width] = size(dog_pyr{i}(:,:,1));
    % find extrema in middle intvls
    for j = 2:s+1
        dog_imgs = dog_pyr{i};
        dog_img = dog_imgs(:,:,j);
        for x = img_border+1:height-img_border
            for y = img_border+1:width-img_border
                % preliminary check on contrast
                if(abs(dog_img(x,y)) > prelim_contr_thr)
                    % check 26 neighboring pixels
                    if(isExtremum(j,x,y))
                        ddata = interpLocation(dog_imgs,height,width,i,j,x,y,img_border,contr_thr,max_interp_steps);
                        if(~isempty(ddata))
                            if(~isEdgeLike(dog_img,ddata.x,ddata.y,curv_thr))
                                 ddata_array(ddata_index) = ddata;
                                 ddata_index = ddata_index + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

function [ flag ] = isExtremum( intvl, x, y)
% Function: Find Extrema in 26 neighboring pixels
    value = dog_imgs(x,y,intvl);
    block = dog_imgs(x-1:x+1,y-1:y+1,intvl-1:intvl+1);
    if ( value > 0 && value == max(block(:)) )
        flag = 1;
    elseif ( value == min(block(:)) )
        flag = 1;
    else
        flag = 0;
    end
end

%% Orientation Assignment
% number of detected points
n = size(ddata_array,2);
% determines gaussian sigma for orientation assignment
ori_sig_factr = 1.5;
% number of bins in histogram
ori_hist_bins = 36;
% orientation magnitude relative to max that results in new feature
ori_peak_ratio = 0.8;
% array of feature
features = struct('ddata_index',0,'x',0,'y',0,'scl',0,'ori',0,'descr',[]);
feat_index = 1;
for i = 1:n
    ddata = ddata_array(i);
    ori_sigma = ori_sig_factr * ddata.scl_octv;
    % generate a histogram for the gradient distribution around a keypoint
    hist = oriHist(gauss_pyr{ddata.octv}(:,:,ddata.intvl),ddata.x,ddata.y,ori_hist_bins,round(3*ori_sigma),ori_sigma);
    for j = 1:2
        smoothOriHist(hist,ori_hist_bins);
    end
    % generate feature from ddata and orientation hist peak
    % add orientations greater than or equal to 80% of the largest orientation magnitude
    feat_index = addOriFeatures(i,feat_index,ddata,hist,ori_hist_bins,ori_peak_ratio);
end

%% Descriptor Generation
% number of features
n = size(features,2);
% width of 2d array of orientation histograms
descr_hist_d = 4;
% bins per orientation histogram
descr_hist_obins = 8;
% threshold on magnitude of elements of descriptor vector
descr_mag_thr = 0.2;
descr_length = descr_hist_d*descr_hist_d*descr_hist_obins;
local_features = features;
local_ddata_array = ddata_array;
local_gauss_pyr = gauss_pyr;
clear features;
clear ddata_array;
clear gauss_pyr;
clear dog_pyr;
parfor feat_index = 1:n
    feat = local_features(feat_index);
    ddata = local_ddata_array(feat.ddata_index);
    gauss_img = local_gauss_pyr{ddata.octv}(:,:,ddata.intvl);
% computes the 2D array of orientation histograms that form the feature descriptor
    hist_width = 3*ddata.scl_octv;
    radius = round( hist_width * (descr_hist_d + 1) * sqrt(2) / 2 );
    feat_ori = feat.ori;
    ddata_x = ddata.x;
    ddata_y = ddata.y;
    hist = zeros(1,descr_length);
    for i = -radius:radius
        for j = -radius:radius
            j_rot = j*cos(feat_ori) - i*sin(feat_ori);
            i_rot = j*sin(feat_ori) + i*cos(feat_ori);
            r_bin = i_rot/hist_width + descr_hist_d/2 - 0.5;
            c_bin = j_rot/hist_width + descr_hist_d/2 - 0.5;
            if (r_bin > -1 && r_bin < descr_hist_d && c_bin > -1 && c_bin < descr_hist_d)
                mag_ori = calcGrad(gauss_img,ddata_x+i,ddata_y+j);
                if (mag_ori(1) ~= -1)
                    ori = mag_ori(2);
                    ori = ori - feat_ori;
                    while (ori < 0)
                        ori = ori + 2*pi;
                    end
                    % i think it's theoretically impossible
                    while (ori >= 2*pi)
                        ori = ori - 2*pi;
                        disp('###################what the fuck?###################');
                    end
                    o_bin = ori * descr_hist_obins / (2*pi);
                    w = exp( -(j_rot*j_rot+i_rot*i_rot) / (2*(0.5*descr_hist_d*hist_width)^2) );
                    hist = interpHistEntry(hist,r_bin,c_bin,o_bin,mag_ori(1)*w,descr_hist_d,descr_hist_obins);
                end
            end
        end
    end
    local_features(feat_index) = hist2Descr(feat,hist,descr_mag_thr);
end
% sort the descriptors by descending scale order
features_scl = [local_features.scl];
[~,features_order] = sort(features_scl,'descend');
% return descriptors and locations
descrs = zeros(n,descr_length);
locs = zeros(n,2);
for i = 1:n
    descrs(i,:) = local_features(features_order(i)).descr;
    locs(i,1) = local_features(features_order(i)).x;
    locs(i,2) = local_features(features_order(i)).y;
end

end

repo: sw-sift


text:
function obj = compile(obj, export_path, varargin)
    % export the symbolic expressions of the system dynamics matrices and
    % vectors and compile as MEX files.
    %
    % Parameters:
    %  export_path: the path to export the file @type char
    %  varargin: variable input parameters @type varargin
    %   StackVariable: whether to stack variables into one @type logical
    %   File: the (full) file name of exported file @type char
    %   ForceExport: force the export @type logical
    %   BuildMex: flag whether to MEX the exported file @type logical
    %   Namespace: the namespace of the function @type char

    input_categories = {'Control','ConstraintWrench','External'};
    for k=1:3
        category = input_categories{k};
        % export the input map
        gmap_funcs = fieldnames(obj.Gmap.(category));
        if ~isempty(gmap_funcs)
            for i=1:length(gmap_funcs)
                fun = gmap_funcs{i};
                if ~isempty(obj.Gmap.(category).(fun))
                    export(obj.Gmap.(category).(fun),export_path,varargin{:});
                end
            end
        end

        % export the input vector fields
        gvec_funcs = fieldnames(obj.Gvec.(category));
        if ~isempty(gvec_funcs)
            for i=1:length(gvec_funcs)
                fun = gvec_funcs{i};
                if ~isempty(obj.Gvec.(category).(fun))
                    export(obj.Gvec.(category).(fun),export_path,varargin{:});
                end
            end
        end
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function results = report_results(base_dir,out_full_name)

files = dir(fullfile(base_dir,'*.mat'));

PSNR = zeros(length(files),1);
MAE = zeros(size(PSNR));
names = {files(:).name};
for i = 1 : length(files)
    load(fullfile(base_dir,files(i).name));
    PSNR(i) = result.PSNR;
end

mean_PSNR = mean(PSNR);
q1_PSNR = q1(PSNR);
q2_PSNR = q2(PSNR);
q3_PSNR = q3(PSNR);

results.PSNR = PSNR; results.names = names;

results.mean_PSNR = mean_PSNR;

results.q1_PSNR = q1_PSNR;

results.q2_PSNR = q2_PSNR;

results.q3_PSNR = q3_PSNR;

save(out_full_name,'results','-v7.3');

metrics = {'PSNR'};
for e = 1 : length(metrics)
    fprintf('%s:\n',metrics{e});
    mean_ = eval(sprintf('mean(results.%s);',metrics{e}));
    q1_ = eval(sprintf('q1(results.%s);',metrics{e}));
    q2_ = eval(sprintf('q2(results.%s);',metrics{e}));
    q3_ = eval(sprintf('q3(results.%s);',metrics{e}));
    fprintf('mean = %.2f, Q1 = %.2f, Q2 = %.2f, Q3 = %.2f\n',...
        mean_,q1_,q2_,q3_);
end

end

function a = q1 (error)
a = quantile(error,0.25);
end

function a = q2 (error)
a = quantile(error,0.5);
end

function a = q3 (error)
a = quantile(error,0.75);
end

author: mahmoudnafifi


text:
function [ degrees, minutes, seconds] = latlonvec(decimal_degrees)
%LATLONVEC Convert decimal degrees into degrees/minutes/seconds
%
% Description:
% Return degrees, minutes, seconds given decimal degrees. Sign of
% coordinate will be returned in degrees portion of coordinate.
%
% Sample fuction call:
% [ degrees, minutes, seconds] = latlonvec(decimal_degrees)
%
% Author: Scott Lee
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

degrees = fix(abs(decimal_degrees));	% integer degrees
frac = abs(decimal_degrees) - degrees;	% fractional degrees, used to compute minutes
minutes = fix(frac*60);	    % integer minutes
frac = frac - minutes/60;	% fractional minutes, used to compute seconds
seconds = frac*3600;		% decimal seconds

% Handle sign.  Degrees portion will contain the sign of the coordinate.
% Minutes and seconds will always be positive.
% sign function returns -1, 0, +1 for x < 0, x == 0, x > 0, respectively
degrees = sign(decimal_degrees).*degrees;

if nargout<2
  degrees = [degrees minutes seconds];
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%              ==== VERY IMPORTANT ====
% This code needs the support of Tensor Toolbox developed by Tamara Kolda
%  which is available at:
%         http://www.sandia.gov/~tgkolda/TensorToolbox/index-2.5.html
%
clear;
clc;
I=[50,50,50];
R=[5,6,7];
N=numel(I);

% Generate data;
A=cell(N,1);
for n=1:N
    A{n}=rand(I(n),R(n));
end
Y=ttensor(tensor(rand(R)),A);
Y=tensor(Y);

opts=struct('NumOfComp',R,'nlssolver','hals','maxiter',100,'maxiniter',20,'tdalgFile','call_tucker_als_opts.mat');
tic;
[Ydec]=lraNTD_ANLS(Y,opts);
toc;
fprintf('Complete. Fit=%f\n',fitness(Y,Ydec));

author: andrewssobral

repo: lrslibrary


text:
function dat  =  testing(algo,dat)

%  if isempty(algo.child.dat) % maybe used a custom ker to train and then changed!
%    kerMaTemp=get_kernel(algo.child,dat,algo.xSV);
%  else
%    kerMaTemp=test(algo.child,dat);
%    kerMaTemp=kerMaTemp.X;
%  end

  %% <<---To avoid large kernel matrices, we test in batches--->
  sz=get_dim(algo.Xsv);   %% <---500x500 point are the maximum for one batch
  if sz==0 sz=1; end;
  batch_size=round((500^2)/sz);

  yEst=[];
  for i=[1:batch_size:get_dim(dat)]
    take= [i:min(i+batch_size-1,get_dim(dat))];
    if ~isempty(algo.alpha)
     kerMaTemp=get_kernel(algo.child,get(dat,take),algo.Xsv);
     yEst=[yEst; ((algo.alpha'* kerMaTemp)+algo.b0)'];
    else
        yEst=[yEst ; algo.b0*ones(length(take),1)];
    end
  end

  if algo.algorithm.use_signed_output==1
    yEst=sign(yEst);
  end

  dat=set_x(dat,yEst);
  dat=set_name(dat,[get_name(dat) ' -> ' get_name(algo)]);


text:
function value = r8_mop ( i )

%*****************************************************************************80
%
%% R8_MOP returns the I-th power of -1 as an R8 value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the power of -1.
%
%    Output, real R8_MOP, the I-th power of -1.
%
  if ( mod ( i, 2 ) == 0 )
    value = + 1.0;
  else
    value = - 1.0;
  end

  return
end


text:
function NodeValue = hlp_computeGraphMeasure(varargin)
% NodeValue = hlp_computeGraphMeasure(causality,ch1,selectedvars,graphMetric)
%
% Compute a univariate graph-theoretic measure for a given node of a graph
%
% Inputs:
%
%   causality:      [num_vars_to x num_vars_from x <num_times> x <num_freqs>] causal matrix
%                   obtained from est_mvarConnectivity().
%   srcNodes:       indices of source nodes
%   selectedvars:   indices of target nodes
%   graphMetric:   which graph measure to compute (see below)
%
% Outputs:
%
%   NodeValue:      [srcNodes x <num_times> x <num_freqs>] graph measures for the source nodes
%
%
% See Also: vis_causalBrainMovie3D(), est_mvarConnectivity(),
%           hlp_collapseFrequencies()
%
%
% References:
%
%   Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 6.
%   Available at: http://www.sccn.ucsd.edu/wiki/SIFT
%
% Author: Tim Mullen, 2010-2013, SCCN/INC/UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

arg_define([0 Inf],varargin, ...
    arg_norep({'cmatrix','CausalMatrix'}, mandatory,[],'Causality array. Shape is [num_vars_to x num_vars_from x <num_times> x <num_freqs>]. Last two dims are optional','shape','matrix'), ...
    arg({'srcNodes','SourceNodes'},  [],[],'Indices of source nodes. Default is all nodes'), ...
    arg({'targNodes','TargetNodes'}, [],[],'Indices of target nodes. Default is all nodes'), ...
    arg({'graphMetric','GraphMetric'},'none',{'none','outflow','mag_outflow','inflow', ...
                                              'causalflow','outdegree','indegree',     ...
                                              'causaldegree','asymmetryratio'}, ...
                                              'Graph measure to compute'), ...
    arg({'ignoreSelfConn','IgnoreSelfConn'},true,[],'Ignore self connectivity in graph measure'));

% handle defaults
if isempty(srcNodes)
    srcNodes = 1:size(cmatrix,2);
end
if isempty(targNodes)
    targNodes = 1:size(cmatrix,1);
end

% remove source node from target nodes
if ignoreSelfConn
    cmatrix=hlp_setdiags(cmatrix,0);
%     targNodes = setdiff_bc(targNodes,srcNodes);
end

% compute graph metrics
switch lower(graphMetric)
    case 'none'
        NodeValue = zeros(1,size(cmatrix,3));
    case 'outflow'
        % Compute outflow from srcNodes in each freq
        NodeValue = squeeze(sum(cmatrix(targNodes,srcNodes,:,:),1));
    case 'mag_outflow'
        % Compute outflow from srcNodes in each freq, ignoring sign
        NodeValue = squeeze(sum(abs(cmatrix(targNodes,srcNodes,:,:)),1));
    case 'inflow'
        % Compute inflow to srcNodes in each freq
        NodeValue = squeeze(sum(cmatrix(srcNodes,targNodes,:,:),2));
    case 'causalflow'
        outflow =   squeeze(sum(cmatrix(targNodes,srcNodes,:,:),1));
        inflow  =   squeeze(sum(cmatrix(srcNodes,targNodes,:,:),2));
        NodeValue = outflow - inflow;
    case 'outdegree'
        % Compute number of outgoing edges from srcNodes in each freq
        NodeValue = squeeze(sum(logical(cmatrix(targNodes,srcNodes,:,:)),1));
    case 'indegree'
        % number of incoming edges to srcNodes in each freq
        NodeValue = squeeze(sum(logical(cmatrix(srcNodes,targNodes,:,:)),2));
    case 'causaldegree'
        % outdegree - indegree
        outflow =   squeeze(sum(logical(cmatrix(targNodes,srcNodes,:,:)),1));
        inflow  =   squeeze(sum(logical(cmatrix(srcNodes,targNodes,:,:)),2));
        NodeValue = outflow - inflow;
    case 'asymmetryratio'
        % 1 if all edges are outgoing, -1 if all edges are incoming.
        % 0 if balanced
        outflow =   squeeze(sum(cmatrix(targNodes,srcNodes,:,:),1));
        inflow  =   squeeze(sum(cmatrix(srcNodes,targNodes,:,:),2));
        NodeValue = (outflow - inflow)./(outflow+inflow);
    otherwise
        % user wants to map a different Conn measure to this
        % (e.g., ERSP)
end

% enforce row vector output
if iscolumn(NodeValue)
    NodeValue = NodeValue';
end

repo: REST


text:
function encmodel = standard_encoding(X,Y,T,options,binsize)
% Compute maps representing the "encoding" model for each
% time point (or window, see below) in the trial (i.e. a GML).
% The reported statistic is the explained variance from
% regressing Y on the data at this sensor/voxel.
%
% INPUT
% X: Brain data, (time by regions) or (time by trials by regions)
% Y: Stimulus, (time by q); q is no. of stimulus features OR
%              (no.trials by q), meaning that each trial has a single
%              stimulus value
% T: Length of series
% options: structure with the preprocessing options - see documentation in
%                       https://github.com/OHBA-analysis/HMM-MAR/wiki
% binsize: how many consecutive time points will be used for each estiamtion.
%           By default, this is 1, such that one encoding model
%               is estimated using 1 time point
%
% OUTPUT
% encmodel:  (time points by regions) maps of activation (explained variances)
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

options.Nfeatures = 0;
options.pca = 0;
options.embeddedlags = 0;
options.K = 1;
[X,Y,T] = preproc4hmm(X,Y,T,options);

if (binsize/2)==0
    warning(['binsize must be an odd number, setting to ' num2str(binsize+1)])
    binsize = binsize + 1;
end

N = length(T); ttrial = T(1);
p = size(X,2); q = size(Y,2);
X = reshape(X,[ttrial N p]);
Y = reshape(Y,[ttrial N q]);

if any(T~=ttrial), error('All trials must have the same length'); end
if nargin<5, binsize = 1; end

encmodel = NaN(ttrial,p);

for t = halfbin+1 : ttrial-halfbin
    r = t-halfbin:t+halfbin;
    for j = 1:p
        Xj = X(t,:,j);
        Xj = zscore(Xj(:));
        Yj = reshape(Y(r,:,:),binsize*N,q);
        beta = (Yj' * Yj) \ (Yj' * Xj);
        res = (Xj - Yj * beta).^2;
        res0 = Xj.^2;
        encmodel(t,j) = 1 - sum(res(:))/sum(res0(:));
    end
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function p = bk_kernel_q(k, l, K, L)
% This program generates probability for each experts
% The defaulty value in BK is uniform, thus, p is straightforward.
% In reality, this can be much more complex
%
% function p = bk_kernel_q(k, l, K, L)
%
% p: vector of weights for individual expert.
%
% k, l: specific parameters to tune the probability distribution.
% K, L: maximum values of k, l, respectively.
%
% Example: p = bk_kernel_q(k, l, K, L);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = 1/(K*L+1);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
% Test file for @chebfun/unwrap.m.

function pass = test_unwrap(pref)

if (nargin < 1)
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% Test empty case.
pass(1) = isempty(unwrap(chebfun()));

% Check that chebfuns with only one fun aren't unwrapped.
f = chebfun(@(x) exp(x), [-1 1], pref);
pass(2) = isequal(f, unwrap(f));

pref.techPrefs.extrapolate = 1;
f = chebfun(@mysawtooth, [0 2*pi 4*pi 6*pi 8*pi], pref);
uf = unwrap(f);
pass(3) = norm(feval(uf, x) - (x - pi), inf) < 10*vscale(uf)*eps;

% Check tol input.
g = 10*f;
ug = unwrap(g, 10*pi);
pass(4) = norm(feval(ug, x) - 10*(x - pi), inf) < 10*vscale(ug)*eps;

% Check error conditions.
try
    f = chebfun(@(x) [sin(x) cos(x)], [-1 0 1]);
    g = unwrap(f);
    pass(5) = false;
catch ME
    pass(5) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:unwrap:array');
end

end

function y = mysawtooth(x)
    y = mod(x, 2*pi) - pi;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nodes, edges] = createTestGraph02(varargin)
%CREATETESTGRAPH02  Test graph for geodesic functions
%
%   [nodes edges] = createTestGraph02;
%
%   Example
%     [nodes edges] = createTestGraph02;
%     figure;
%     axis([0 100 10 90]);
%     axis equal;
%     drawGraph(nodes, edges);
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

nodes = [ ...
   10 40; ...
   10 60; ...
   20 50; ...
   40 50; ...
   50 20; ...
   50 40; ...
   50 60; ...
   50 80; ...
   60 50; ...
   80 50; ...
   90 40; ...
   90 60];

edges = [...
    1  3; ...
    2  3; ...
    3  4; ...
    4  6; ...
    4  7; ...
    5  6; ...
    6  9; ...
    7  8; ...
    7  9; ...
    9 10; ...
   10 11; ...
   10 12];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function Test_fixedrankembeddedfactory()
% Given partial observation of a low rank matrix, attempts to complete it.
%
% function low_rank_matrix_completion()
%
% This example demonstrates how to use the geometry factory for the
% embedded submanifold of fixed-rank matrices, fixedrankembeddedfactory.
% This geometry is described in the paper
% "Low-rank matrix completion by Riemannian optimization"
% Bart Vandereycken - SIAM Journal on Optimization, 2013.
%
% This can be a starting point for many optimization problems of the form
%
% minimize f(X) such that rank(X) = k, size(X) = [m, n].
%
% Input:  None. This example file generate random data.
%
% Output: None.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 15, 2014
% Contributors:
%
% Change log:
%

    % Choose the size of the problem:
    % We will complete a matrix of size mxn of rank k.
    m = 2000;
    n = 5000;
    k = 10;
    % Generate a random mxn matrix A of rank k
    L = randn(m, k);
    R = randn(n, k);
    A = L*R';
    % Generate a random mask for observed entries: P(i, j) = 1 if the entry
    % (i, j) of A is observed, and 0 otherwise.
    fraction = 4 * k*(m+n-k)/(m*n);
    P = sparse(rand(m, n) <= fraction);
    % Hence, we know the nonzero entries in PA:
    PA = P.*A;

    % Pick the manifold of matrices of size mxn of fixed rank k.
    problem.M = fixedrankembeddedfactory(m, n, k);

    % Define the problem cost function. The input X is a structure with
    % fields U, S, V representing a rank k matrix as U*S*V'.
    % f(X) = 1/2 * || P.*(X-A) ||^2
    problem.cost = @cost;
    function f = cost(X)
        % Note that it is very much inefficient to explicitly construct the
        % matrix X in this way. Seen as we only need to know the entries
        % of Xmat corresponding to the mask P, it would be far more
        % efficient to compute those only.
        Xmat = X.U*X.S*X.V';
        f = .5*norm( P.*Xmat - PA , 'fro')^2;
    end

    % Define the Euclidean gradient of the cost function, that is, the
    % gradient of f(X) seen as a standard function of X.
    % nabla f(X) = P.*(X-A)
    problem.egrad = @egrad;
    function G = egrad(X)
        % Same comment here about Xmat.
        Xmat = X.U*X.S*X.V';
        G = P.*Xmat - PA;
    end

    % Define the Euclidean Hessian of the cost at X, along H, where H is
    % represented as a tangent vector.
    % nabla^2 f(X)[Xdot] = P.*Xdot
    % (This is the directional derivative of nabla f(X) at X along Xdot.)
    problem.ehess = @euclidean_hessian;
    function ehess = euclidean_hessian(X, H)
        % The function tangent2ambient transforms H (a tangent vector) into
        % its equivalent ambient vector representation. The output is a
        % structure with fields U, S, V such that U*S*V' is an mxn matrix
        % corresponding to the tangent vector H. Note that there are no
        % additional guarantees about U, S and V. In particular, U and V
        % are not orthonormal.
        ambient_H = problem.M.tangent2ambient(X, H);
        Xdot = ambient_H.U*ambient_H.S*ambient_H.V';
        % Same comment here about explicitly constructing the ambient
        % vector as an mxn matrix Xdot: we only need its entries
        % corresponding to the mask P, and this could be computed
        % efficiently.
        ehess = P.*Xdot;
    end

    % Check consistency of the gradient and the Hessian. Useful if you
    % adapt this example for a now cost function and you would like to make
    % sure there is no mistake.
    % warning('off', 'manopt:fixedrankembeddedfactory:exp');
    % checkgradient(problem); pause;
    % checkhessian(problem); pause;

    % Compute an initial guess. Points on the manifold are represented as
    % structures with three fields: U, S and V. U and V need to be
    % orthonormal, S needs to be diagonal.
    [U, S, V] = svds(PA, k);
    X0.U = U;
    X0.S = S;
    X0.V = V;

    % Minimize the cost function using Riemannian trust-regions, starting
    % from the initial guess X0.
    options.stopfun = @stopifclosedfigure;
    X = trustregions(problem, X0, options);

    % The reconstructed matrix is X, represented as a structure with fields
    % U, S and V.
    Xmat = X.U*X.S*X.V';
    fprintf('||X-A||_F = %g\n', norm(Xmat - A, 'fro'));

    % If the problem has a small enough dimension, we may (for analysis
    % purposes) compute the spectrum of the Hessian at a point X. This may
    % help in studying the conditioning of a problem.
    if problem.M.dim() < 100
        fprintf('Computing the spectrum of the Hessian...');
        s = hessianspectrum(problem, X);
        hist(s);
    end

end

author: NicolasBoumal

repo: manopt


text:
function nWritten = export_MR_IOD(planC, filenameRoot, filenumber)
%"export_MR_IOD"
%   Writes to disk all scans contained in planC by constructing the modules
%   specified in part A.3.3 of PS 3.3 of the 2006 DICOM standard.  The
%   planC MUST have first been run through the function:
%
%   generate_DICOM_UID_Relationships.m.
%
%   filenameRoot is the the base path + filename to have file numbers
%   appended to, ie, "C:\tmp\DCM".  filenumber is the starting filenumber
%   to append, with results such as "C:\tmp\DCM00000", "C:\tmp\DCM00001"
%   etc.
%
%DK 07/25/08
%
%Usage:
%   nWritten = export_MR_IOD(planC, filenameRoot, filenumber);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS  = planC{end};

%Init number of files written to zero.
nWritten = 0;

for scanNum = 1:length(planC{indexS.scan})

    %Get the scan we are currently working with
    scanS   = planC{indexS.scan}(scanNum);

    %Build patient module.
    patientobj      = export_module('patient', 'scan', scanS);

    %Build general study module.
    studyobj        = export_module('general_study', scanS);

    %Build general series module.
    seriesobj       = export_module('general_series', scanS);

    %Build a frame of reference module.
    frameobj        = export_module('frame_of_reference', scanS);

    %Build a general equipment module.
    equipobj        = export_module('general_equipment', 'scan', scanS);

    %For slice-specific modules iterate over scaninfo.
    for i=1:length(planC{indexS.scan}.scanInfo)

        %Create a dcmobj to hold a single slice.
        dcmobj = org.dcm4che2.data.BasicDicomObject;

        %Get info for the slice we are handling.
        scanInfoS = scanS.scanInfo(i);

        %Build an image module from this particular slice (scanInfoS)
        imgobj          = export_module('general_image', 'scan', scanInfoS);

        %Build an image plane module from this particular slice (scanInfoS)
        imgplaneobj     = export_module('image_plane', 'scan', scanInfoS, scanS);

        %Build an image pixel module from this particular slice (scanInfoS)
        imgpixelobj     = export_module('image_pixel', 'scan', scanInfoS, scanS);

        %Build an MR image module from this particular slice (scanInfoS)
        mrimageobj      = export_module('MR_image', scanInfoS, scanS);

        SOPobj = export_module('SOP_common', 'scanInfo', scanInfoS);

        %Combine all modules into a single dcmobj.
        patientobj.copyTo(dcmobj);
        studyobj.copyTo(dcmobj);
        seriesobj.copyTo(dcmobj);
        frameobj.copyTo(dcmobj);
        equipobj.copyTo(dcmobj);
        imgobj.copyTo(dcmobj);
        imgplaneobj.copyTo(dcmobj);
        imgpixelobj.copyTo(dcmobj);
        mrimageobj.copyTo(dcmobj);
        SOPobj.copyTo(dcmobj);

        clear imgobj imgplaneobj imgpixelobj mrimageobj SOPobj

        fileNum  = num2str(filenumber + nWritten);
        filename = [filenameRoot repmat('0', [1 5-length(fileNum)]) fileNum];

        writefile_mldcm(dcmobj, filename);

        nWritten = nWritten + 1;

        clear dcmobj;
    end

    clear patientobj studyobj seriesobj frameobj equipobj

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function h = disp_cubes(cubes, img, mode)

% cubes = drop_thin_cuboids(cubes, vp);

if nargin<3
    mode = 1;
end

if ~isempty(img)
    figure; imshow(img); hold on;
end

h = [];
for i = 1:length(cubes)
    h = [h disp_cube_2(cubes(i))];
    if mode==3
        pause; delete(h); h = [];
    elseif mode==2
        pause(0.01); delete(h); h = [];
    end
end

%%
function h = disp_cube_2(cube)

% color = rand(1,3);
color = [1 0.5 0];
% color = [1 1 0];
linewidth = 5;

h = zeros(1,12);
h(1) = plotline(cube.junc3(1).pt, cube.junc3(2).pt, color, linewidth);
h(2) = plotline(cube.junc3(1).pt, cube.junc3(3).pt, color, linewidth);
h(3) = plotline(cube.junc3(1).pt, cube.junc3(5).pt, color, linewidth);
h(4) = plotline(cube.junc3(2).pt, cube.junc3(4).pt, color, linewidth);
h(5) = plotline(cube.junc3(2).pt, cube.junc3(6).pt, color, linewidth);
h(6) = plotline(cube.junc3(3).pt, cube.junc3(4).pt, color, linewidth);
h(7) = plotline(cube.junc3(3).pt, cube.junc3(7).pt, color, linewidth);
h(8) = plotline(cube.junc3(4).pt, cube.junc3(8).pt, color, linewidth);
h(9) = plotline(cube.junc3(5).pt, cube.junc3(6).pt, color, linewidth);
h(10) = plotline(cube.junc3(5).pt, cube.junc3(7).pt, color, linewidth);
h(11) = plotline(cube.junc3(6).pt, cube.junc3(8).pt, color, linewidth);
h(12) = plotline(cube.junc3(7).pt, cube.junc3(8).pt, color, linewidth);

% for j = 1:8
%     h(12+j) = text(cube.junc3(j).pt(1), cube.junc3(j).pt(2), num2str(j), 'Color', [1 0 0]);
% end

author: zouchuhang

repo: LayoutNet


text:
function resu = ismemberclique(v,cliques)
% b = ismemberclique(v,cliques)
%

finiclique = 0 ; resu=0 ;
cl=1;  ncl=length(cliques) ;

while (~finiclique) & (cl<=ncl);
    if ismember(v,cliques{cl})
        resu=1;
        finiclique=1 ;
    end
    cl=cl+1;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt