text:
function  [results] = VerifyError(expr, results)
wasError = 0;
try
    eval(expr);
catch ex
    wasError=1;
end
if wasError == 0
    disp ([expr ' <-- expected error, but was no error'])
    count = 0;
    if(isfield(results, 'failedTests'))
        count = size(results.failedTests, 1);
    end
    results.failedTests{count+1,1} = expr;
    results.allTestsPassed = 0;
end
return

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out=zeros(in,varargin)

if nargin==1
 out=mp(zeros(double(in)));
else
 out=mp(zeros([double(in),double([varargin{:}])]));
end

author: opencobra

repo: cobratoolbox


text:
 function y = mtimes(ob, x)
%function y = mtimes(ob, x)
% A * x
% A' * y
% y' * A
% x' * A'
% B * A		(calls private/fatrix2_mtimes2)
% etc

if isnumeric(ob)
	if isscalar(ob) % scalar * object
		y = fatrix2_scalar_times(ob, x);
	else % row_vector(s) * object
		y = (x' * ob')'; % trick: x <-> ob
	end
return
end

if ~isnumeric(x) % object * object
	if islogical(x)
		warn 'fatrix2 * logical is illogical; reconsider!'
		x = single(x);
%	elseif iscell(x) % untested feature for mai le
		% do nothing here; continue to "ordinary" multiplication below
	else
		y = fatrix2_mtimes2(ob, x);
	return
	end
end

y = fatrix2_mtimes_vector(ob, x); % "ordinary" fatrix * vector multiplication

end % mtimes()

% fatrix2_mtimes_vector()
% ordinary matrix-vector multiplication, but generalized to support arrays
% full multiplication
% y = scale * diag(odiag) * ob * diag(idiag) * x;
function y = fatrix2_mtimes_vector(ob, x)

fun_forw = @(ob, x) fatrix2_do_forw(ob, x);

[nd np] = size(ob);
dimx = size(x);
if dimx(1) == np % [np (L)] vector mode
	diml = [dimx(2:end) 1];
	LL = prod(diml);
	if LL == 1 % usual single column vector case
		x = iembed(ob, x); % [idim]
		y = fun_forw(ob, x); % [odim]
		y = y(:); % [*odim 1]
		y = oselect(ob, y); % [nd 1]
	elseif ob.does_many % multiple column vector without loop
		x = reshape(x, [np LL]); % [np *L]
		x = iembed(ob, x); % [idim *L]
		y = fun_forw(ob, x); % [odim *L]
		y = reshape(y, [prod(ob.odim) LL]); % [*odim *L]
		y = oselect(ob, y); % [nd *L]
		y = reshape(y, [nd diml]); % [nd (L)]
	else % multiple column vector case with loop
		x = reshape(x, [np LL]); % [np *L]
		tmp = iembed(ob, x(:,1)); % [idim]
		tmp = fun_forw(ob, tmp); % [odim]
		y = zeros([prod(ob.odim) LL], class(tmp)); % [*odim *L]
		y(:,1) = tmp(:);
		for ll=2:LL
			tmp = iembed(ob, x(:,ll)); % [idim]
			tmp = fun_forw(ob, tmp); % [odim]
			y(:,ll) = tmp(:);
		end
		y = oselect(ob, y); % [nd *L]
		y = reshape(y, [nd diml]); % [nd (L)]
	end

elseif first_dim_match(dimx, ob.idim)

	diml = [dimx((numel(ob.idim)+1):end) 1];
	LL = prod(diml);
	if LL == 1 % a single array
		y = fun_forw(ob, x); % [odim]
	elseif ob.does_many
		x = reshape(x, [ob.idim LL]); % [idim *L]
		y = fun_forw(ob, x); % [odim *L]
		y = reshape(y, [ob.odim diml]); % [odim (L)]
	else % loop
		x = reshape(x, [prod(ob.idim) LL]); % [*idim *L]
		tmp = reshape(x(:,1), [ob.idim 1]); % [idim]
		tmp = fun_forw(ob, tmp); % [odim]
		y = zeros([prod(ob.odim) LL], class(tmp)); % [*odim *L]
		y(:,1) = tmp(:);
		for ll=2:LL
			tmp = reshape(x(:,ll), [ob.idim 1]); % [idim]
			tmp = fun_forw(ob, tmp); % [odim]
			y(:,ll) = tmp(:);
		end
		y = reshape(y, [ob.odim diml]); % [odim (L)]
	end

else
	pr size(x)
	pr ob.idim
	fail 'input dimension size mismatch'
end

end % fatrix2_mtimes_vector()

% first_dim_match()
% see if first dimensions of dimx match idim properly
% basically, dimx = [idim (L)]
function yn = first_dim_match(dimx, idim)
yn = numel(dimx) >= numel(idim) && isequal(dimx(1:numel(idim)), idim);
if ~yn && idim(end) == 1 % handle case [... 1]
	yn = first_dim_match(dimx, idim(1:end-1));
end
end % first_dim_match()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%rate_a=ach_simo_nocsi(nn,P,error,rx,K) computes the achievability bound on
%the maximal channel coding rate over a single-input multiple-output (SIMO)
%quasi-static Rician fading channel with no channel state information at
%transmitter or receiver.
%
%The function takes the following inputs:
%
%nn: blocklength, scalar/vector;
%P: input power (in linear scale), scalar;
%error: target block error probability, scalar, should be greater than 10^(-6);
%rx: number of receive antennas, scalar;
%K: rician K-factor, scalar; the default value is 0.
function rate_a=ach_simo_nocsi(nn,P,error,rx,K)

if (nargin < 5) || isempty(K)
	K = 0;
end

loop=min(1000/error, 10^7);

rate_a=[];
hh = sqrt(0.5/(K+1))*(randn(rx,loop)+1i*randn(rx,loop)) + sqrt(K/(K+1)); %samples of Rician fading channels
w1 = sqrt(0.5) * (randn(rx, loop) + 1i*randn(rx,loop)); %AWGN noise

for n=nn
    data_a=zeros(1,loop);

    %By spherical symmetry, it suffices to consider x0=[(nP)^{-1/2},0,...,0]
    y1 = hh.*sqrt(n*P) + w1; %The received vector during the first channel use
    for ii=1:loop
        y1_sample = y1(:,ii);
        W = sqrt(0.5) * (randn(rx, n-1) + 1i*randn(rx,n-1));

        %projection of x0 onto Y
        data_a(ii) = abs(y1_sample'* inv(y1_sample*y1_sample' + W*W')*y1_sample);
        %
    end

    data_sort = sort(data_a);

    %maximize over tau...
    tau_array = [floor(error*loop)-1:-1:1]/loop ;

    gamma_tau = 1 - data_sort(1:floor(error*loop)-1);

    log_Fn = log(betainc(gamma_tau, n-rx,rx));
    index_set= find (log_Fn < -10^100);%detect if the output of betainc is zero

    if length(index_set)>0
       log_Fn(index_set) = (n-rx) * log(gamma_tau(index_set))  - sum(log(n-rx+1:1:n-1)) + log(gamma(rx)) ;
    end

    %
    log_M = max(log(tau_array) - log_Fn );

    rate_a =[rate_a, log_M/n/log(2)];
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function imResult = blendMode_SoftLight(A, B, offsetW, offsetH)
%% Soft Light blending mode: combines Multiply and Screen blend modes.
%   This is a softer version of Hard Light.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_SoftLight));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

ind = B < 0.5;
indCompl = abs(ind - 1);

tmp1 = 2 * A .* B + A.^2 .* (1 - 2 * B);
tmp2 = 2 * A .* (1 - B) + sqrt(A) .* (2 * B - 1);
C = ind .* tmp1 + indCompl .* tmp2;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end

author: Sable

repo: mcbench-benchmarks


text:
%io_readpta.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% [rf,info]=io_readpta(filename)
%
% DESCRIPTION:
% Read a Siemens .pta file into matlab.  The resulting RF matrix will have
% 2 columns specifying magnitude and phase.
%
% INPUTS:
% filename   = filename of the .pta file to read in.
%
% OUTPUTS:
% rf        = Input rf pulse waveform saved as a matlab array with 2
%               columns (magnitude and phase).
% info      = Not used.

function [rf,info]=io_readpta(filename)

%try to incorporate the header information into a structure called 'info'
fid=fopen(filename);
line=fgets(fid);
colon_index=findstr(line,':');
if isempty(colon_index)
    line=fgets(fid);
    colon_index=findstr(line,':');
end

while ~isempty(colon_index)
    fieldname=line(1:colon_index-1);
    info.(genvarname(fieldname))=line(colon_index+1:end);
    line=fgets(fid);
    colon_index=findstr(line,':');
end

% Now begin to read the data.  PTA files have a semicolon marking each
% line of Data.  Search for the semicolon on each line and read only the
%data that preceeds it.
RF=zeros(2,0);
linenum=1;
semicol_index=findstr(line,';');
if isempty(semicol_index)
    line=fgets(fid);
    semicol_index=findstr(line,';');
end

% If the line is empty skip it
while ~isempty(semicol_index)
    dataline=line(1:semicol_index-2);
    [A,count, errmsg, nextindex] = sscanf(dataline, '%f', inf);
    % If read failed, output the error
    if ~isempty(errmsg)
       fclose(fid)
       error('READPTA failed with read error: %s', errmsg);
    end
    % Store the read values into rf array
    RF(:,linenum) = A;
    linenum = linenum + 1;
    line=fgets(fid);
    semicol_index=findstr(line,';');
end
fclose(fid);

RF=RF';
rf(:,1)=RF(:,2)*180/pi;
rf(:,2)=RF(:,1);
rf(:,3)=ones(length(RF(:,1)),1);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  June 2013
% ------------------------------------------------------------------------

function [combined_stats,conv_hull_stats,params]= combine_masks(stats1, stats2, measure)
if nargin<3
    measure = 'jaccard_class';
end

n_res1 = length(stats1.mean_n_masks);
n_res2 = length(stats2.mean_n_masks);

%% Combinations
% Single results
for kk=1:n_res1
    combined_stats{kk,n_res2+1} = extract_one_from_all(stats1,kk); %#ok<AGROW>
end
for ii=1:n_res2
    combined_stats{n_res1+1,ii} = extract_one_from_all(stats2,ii); %#ok<AGROW>
end
% Combined results
for kk=1:n_res1
    for ii=1:n_res2
        combined_stats{kk,ii} = join_masks(combined_stats{kk,n_res2+1}, combined_stats{n_res1+1,ii}); %#ok<AGROW>
    end
end
% No results
combined_stats{n_res1+1,n_res2+1}.mean_n_masks = 0;
combined_stats{n_res1+1,n_res2+1}.(measure) = 0;

%% Get convex hull
nn = 1;
all_n_masks = zeros((n_res1+1)*(n_res2+1),1);
xx = zeros((n_res1+1)*(n_res2+1),1);
yy = zeros((n_res1+1)*(n_res2+1),1);
for kk=1:n_res1+1
    for ii=1:n_res2+1
        all_n_masks(nn) = combined_stats{kk,ii}.mean_n_masks;
        xx(nn) = kk;
        yy(nn) = ii;
        nn = nn+1;
    end
end
all_n_masks(end) = [];
xx(end) = [];
yy(end) = [];

nn = 1;
all_J = zeros((n_res1+1)*(n_res2+1),1);
for kk=1:n_res1+1
    for ii=1:n_res2+1
        all_J(nn)   = combined_stats{kk,ii}.(measure);
        nn = nn+1;
    end
end
all_J(end) = [];

if length(all_n_masks)==1
    sel_stats = 1;
else
    % Pareto front
    sel_stats = find(paretoGroup([all_n_masks, -all_J]));
    [~,tmp_ids] = sort(all_n_masks(sel_stats));
    sel_stats = sel_stats(tmp_ids);

    % Convex hull to simplify
    sel_stats_hull = convhull(all_n_masks(sel_stats), all_J(sel_stats),'simplify',true);
    sel_stats = sel_stats(sel_stats_hull);
end

for ii=1:length(sel_stats)
    c_hull_stats(ii) = combined_stats{xx(sel_stats(ii)), yy(sel_stats(ii))}; %#ok<AGROW>
    params(ii).params1 = xx(sel_stats(ii)); %#ok<AGROW>
    params(ii).params2 = yy(sel_stats(ii)); %#ok<AGROW>
end

%% Put the convex hull in the same struct
% Sorted in order of n_masks
[~, sel_stats_ids] = sort([c_hull_stats.mean_n_masks]);
params = params(sel_stats_ids);
conv_hull_stats = c_hull_stats(sel_stats_ids(1));
conv_hull_stats.gt_set = stats1.gt_set;
for ii=2:length(sel_stats_ids)
    conv_hull_stats.all_n_masks  = [conv_hull_stats.all_n_masks c_hull_stats(sel_stats_ids(ii)).all_n_masks];
    conv_hull_stats.max_J        = [conv_hull_stats.max_J c_hull_stats(sel_stats_ids(ii)).max_J];
    conv_hull_stats.max_fp       = [conv_hull_stats.max_fp c_hull_stats(sel_stats_ids(ii)).max_fp];
    conv_hull_stats.max_fn       = [conv_hull_stats.max_fn c_hull_stats(sel_stats_ids(ii)).max_fn];
    conv_hull_stats.max_inters   = [conv_hull_stats.max_inters c_hull_stats(sel_stats_ids(ii)).max_inters];
    conv_hull_stats.jaccard_object = [conv_hull_stats.jaccard_object c_hull_stats(sel_stats_ids(ii)).jaccard_object];
    conv_hull_stats.mean_n_masks = [conv_hull_stats.mean_n_masks c_hull_stats(sel_stats_ids(ii)).mean_n_masks];
    conv_hull_stats.jaccard_class  = [conv_hull_stats.jaccard_class c_hull_stats(sel_stats_ids(ii)).jaccard_class];

    for kk=1:length(c_hull_stats(sel_stats_ids(ii)).per_class_results)
        conv_hull_stats.per_class_results{kk}.max_J         = [conv_hull_stats.per_class_results{kk}.max_J c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.max_J];
        conv_hull_stats.per_class_results{kk}.max_fp        = [conv_hull_stats.per_class_results{kk}.max_fp c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.max_fp];
        conv_hull_stats.per_class_results{kk}.max_fn        = [conv_hull_stats.per_class_results{kk}.max_fn c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.max_fn];
        conv_hull_stats.per_class_results{kk}.max_inters    = [conv_hull_stats.per_class_results{kk}.max_inters c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.max_inters];
        conv_hull_stats.per_class_results{kk}.meanmax       = [conv_hull_stats.per_class_results{kk}.meanmax c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.meanmax];

        conv_hull_stats.per_class_results{kk}.global_fp     = [conv_hull_stats.per_class_results{kk}.global_fp c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.global_fp];
        conv_hull_stats.per_class_results{kk}.global_fn     = [conv_hull_stats.per_class_results{kk}.global_fn c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.global_fn];
        conv_hull_stats.per_class_results{kk}.global_inters = [conv_hull_stats.per_class_results{kk}.global_inters c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.global_inters];

        conv_hull_stats.per_class_results{kk}.global_J      = [conv_hull_stats.per_class_results{kk}.global_J c_hull_stats(sel_stats_ids(ii)).per_class_results{kk}.global_J];
    end
end

%% Some sanity checks
% Compute global mean
tmp = [];
for ii=1:length(conv_hull_stats.per_class_results)
    tmp = [tmp; conv_hull_stats.per_class_results{ii}.global_J]; %#ok<AGROW>
end
assert(isequal(conv_hull_stats.jaccard_class,mean(tmp)));
% Per-class means
for ii=1:length(conv_hull_stats.per_class_results)
    tmp =   conv_hull_stats.per_class_results{ii}.global_inters ./...
            (conv_hull_stats.per_class_results{ii}.global_inters+...
             conv_hull_stats.per_class_results{ii}.global_fp    +...
             conv_hull_stats.per_class_results{ii}.global_fn);
    tmp(isnan(tmp)) = 0;
    assert(isequal(tmp, conv_hull_stats.per_class_results{ii}.global_J))
end

%% Debug
% subplot(1,2,1)
% plot(all_n_masks,all_J,'b+')
% hold on
% plot(all_n_masks(sel_stats),all_J(sel_stats),'k')

end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function p = createsdpcut(p,x)
if any(p.K.s)
    % Cuts are added from violated relaxations
    % Might improve bound propagation.
    top = startofSDPCone(p.K);
    newcuts = 1;
    newF = [];
    newSOCP = [];
    for i = 1:length(p.K.s)
        n = p.K.s(i);
        X = p.F_struc(top:top+n^2-1,:)*[1;x];
        X = full(reshape(X,n,n));
        [d,v] = eig(X);
        if any(diag(v)<0)
             for m = 1:length(v)
                if v(m,m)<0
                    for j = 1:length(x)+1
                        newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(top:top+n^2-1,j),n,n)*d(:,m);
                    end
                    newF(newcuts,1)=newF(newcuts,1)+1e-6;
                    newcuts = newcuts + 1;
                    if size(p.lpcuts,1)>0
                        dist = p.lpcuts*newF(newcuts-1,:)'/(newF(newcuts-1,:)*newF(newcuts-1,:)');
                        if any(abs(dist-1)<1e-3)
                            newF = newF(1:end-1,:);
                            newcuts = newcuts - 1;
                        end
                    end
                end
            end
        end
        top = top+n^2;
    end

    if ~isempty(newF)
        m = size(newF,2);
        p.lpcuts = [newF;p.lpcuts];
        p.cutState = [ones(size(newF,1),1);p.cutState];
    end
    if ~isempty(newSOCP)
        p.socpcuts.F_struc = [p.socpcuts.F_struc;newSOCP];
        p.socpcuts.K.q = [p.socpcuts.K.q 3];
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function log_warning(str,varargin)
% Sends its inputs to the logger if it exists; otherwise does
% nothing.  The logging level is set to Warning.
% Inputs:
%   str: A string (possibly containing format specifiers).
%   varargin: An optional parameter list of values for the format
%     specifiers in the string.
% See documentation for fprintf for a description of the
% parameters 'str' and 'varargin'.

assert(isa(str,'char'))

outstr = ['Warning: ',str];

logprint(Logger.Warning,outstr,varargin{:});

author: nesl


text:
classdef Stitcher < handle
    %STITCHER  High level image stitcher
    %
    % It's possible to use this class without being aware of the entire
    % stitching pipeline. However, to be able to achieve higher stitching
    % stability and quality of the final images at least being familiar with
    % the theory is recommended.
    %
    % ![image](https://docs.opencv.org/3.3.1/StitchingPipeline.jpg)
    %
    % This figure illustrates the stitching module pipeline implemented in the
    % cv.Stitcher class. Using that class it's possible to configure/remove
    % some steps, i.e. adjust the stitching pipeline according to the
    % particular needs.
    %
    % The implemented stitching pipeline is very similar to the one proposed
    % in [BL07].
    %
    % # Camera models
    % There are currently 2 camera models implemented in stitching pipeline.
    %
    % - *Homography model* expecting perspective transformations between
    %   images implemented in `BestOf2NearestMatcher`,
    %   `HomographyBasedEstimator`, `BundleAdjusterReproj`, and
    %   `BundleAdjusterRay`.
    % - *Affine model* expecting affine transformation with 6 DOF or 4 DOF
    %   implemented in `AffineBestOf2NearestMatcher`, `AffineBasedEstimator`,
    %   `BundleAdjusterAffine`, `BundleAdjusterAffinePartial`, and
    %   `AffineWarper`.
    %
    % Homography model is useful for creating photo panoramas captured by
    % camera, while affine-based model can be used to stitch scans and object
    % captured by specialized devices. Use cv.Stitcher.Stitcher to get
    % preconfigured pipeline for one of those models.
    %
    % Note: Certain detailed settings of cv.Stitcher might not make sense.
    % Especially you should not mix classes implementing affine model and
    % classes implementing Homography model, as they work with different
    % transformations.
    %
    % ## Example
    % * A basic example on image stitching can be found in the
    %   `stitching_demo.m` sample
    % * A detailed example on image stitching can be found in the
    %   `stitching_detailed_demo.m` sample.
    %
    % ## References
    % [BL07]:
    % > Matthew Brown and David G Lowe.
    % > "Automatic panoramic image stitching using invariant features".
    % > International journal of computer vision, 74(1):59-73, 2007.
    %
    % See also: cv.Stitcher.Stitcher
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Registration resolution. default 0.6
        RegistrationResol
        % Seam estimation resolution. default 1
        SeamEstimationResol
        % Compositing resolution. default 'Orig'
        CompositingResol
        % Panorama confidence threshold. default 1
        PanoConfidenceThresh
        % Whether it should try to make panorama more horizontal (or vertical).
        % default true
        WaveCorrection
        % Correction kind. default 'Horiz'
        WaveCorrectKind
    end

    methods
        function this = Stitcher(varargin)
            %STITCHER  Creates a Stitcher configured in one of the stitching modes
            %
            %     obj = cv.Stitcher()
            %     obj = cv.Stitcher('OptionName',optionValue, ...)
            %
            % ## Options
            % * __Mode__ Scenario for stitcher operation. This is usually
            %   determined by source of images to stitch and their
            %   transformation. Default parameters will be chosen for
            %   operation in given scenario. Default 'Panorama'. One of:
            %   * __Panorama__ Mode for creating photo panoramas. Expects
            %     images under perspective transformation and projects
            %     resulting pano to sphere. See also `BestOf2NearestMatcher`,
            %     `SphericalWarper`.
            %   * __Scans__ Mode for composing scans. Expects images under
            %     affine transformation does not compensate exposure by
            %     default. See also `AffineBestOf2NearestMatcher`,
            %     `AffineWarper`
            % * __TryUseGPU__ Flag indicating whether GPU should be used
            %   whenever it's possible. default false
            %
            % See also: cv.Stitcher.stitch
            %
            this.id = Stitcher_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.Stitcher
            %
            if isempty(this.id), return; end
            Stitcher_(this.id, 'delete');
        end
    end

    methods
        function pano = stitch(this, images, varargin)
            %STITCH  Tries to stitch the given images
            %
            %     pano = obj.stitch(images)
            %     pano = obj.stitch(images, rois)
            %     [pano, status] = obj.stitch(...)
            %
            % ## Input
            % * __images__ Input cell-array of images.
            % * __rois__ Optional region of interest rectangles, a cell-array
            %   of cell-arrays of 4-elements vectors `{{[x,y,w,h], ...}, ...}`
            %   or `{[x,y,w,h; ...], ...}`.
            %
            % ## Output
            % * __pano__ Final pano.
            % * __status__ optional output status code. If not requested, the
            %   function throws an error if the operation fails. A string one
            %   of:
            %   * __OK__
            %   * **ERR_NEED_MORE_IMGS**
            %   * **ERR_HOMOGRAPHY_EST_FAIL**
            %   * **ERR_CAMERA_PARAMS_ADJUST_FAIL**
            %
            % The function throws an error if the stitch function returns a
            % non-OK status code.
            %
            % See also: cv.Stitcher.estimateTransform,
            %  cv.Stitcher.composePanorama
            %
            pano = Stitcher_(this.id, 'stitch', images, varargin{:});
        end

        function estimateTransform(this, images, varargin)
            %ESTIMATETRANSFORM  Estimate transformation
            %
            %     obj.estimateTransform(images)
            %     obj.estimateTransform(images, rois)
            %     status = obj.estimateTransform(...)
            %
            % ## Input
            % * __images__ Input cell-array of images.
            % * __rois__ Optional region of interest rectangles, a cell-array
            %   of cell-arrays of 4-elements vectors `{{[x,y,w,h], ...}, ...}`
            %   or `{[x,y,w,h; ...], ...}`.
            %
            % ## Output
            % * __status__ optional output status code. If not requested, the
            %   function throws an error if the operation fails. A string one
            %   of:
            %   * __OK__
            %   * **ERR_NEED_MORE_IMGS**
            %   * **ERR_HOMOGRAPHY_EST_FAIL**
            %   * **ERR_CAMERA_PARAMS_ADJUST_FAIL**
            %
            % This function tries to match the given images and to estimate
            % rotations of each camera.
            %
            % Note: Use the function only if you're aware of the stitching
            % pipeline, otherwise use cv.Stitcher.stitch.
            %
            % See also: cv.Stitcher.composePanorama
            %
            Stitcher_(this.id, 'estimateTransform', images, varargin{:});
        end

        function pano = composePanorama(this, varargin)
            %COMPOSEPANORAMA  Compose panorama
            %
            %     pano = obj.composePanorama()
            %     pano = obj.composePanorama(images)
            %     [pano, status] = obj.composePanorama(...)
            %
            % ## Input
            % * __images__ Input cell-array of images.
            %
            % ## Output
            % * __pano__ Final pano.
            % * __status__ optional output status code. If not requested, the
            %   function throws an error if the operation fails. A string one
            %   of:
            %   * __OK__
            %   * **ERR_NEED_MORE_IMGS**
            %   * **ERR_HOMOGRAPHY_EST_FAIL**
            %   * **ERR_CAMERA_PARAMS_ADJUST_FAIL**
            %
            % This function tries to compose the given images (or images
            % stored internally from the other function calls) into the final
            % `pano` under the assumption that the image transformations were
            % estimated before.
            %
            % Note: Use the function only if you're aware of the stitching
            % pipeline, otherwise use cv.Stitcher.stitch.
            %
            % See also: cv.Stitcher.estimateTransform
            %
            pano = Stitcher_(this.id, 'composePanorama', varargin{:});
        end

        function indices = component(this)
            %COMPONENT  Image indices
            %
            %     indices = obj.component()
            %
            % ## Output
            % * __indices__ Vector of integer indices (0-based).
            %
            % See also: cv.Stitcher.cameras, cv.Stitcher.workScale
            %
            indices = Stitcher_(this.id, 'component');
        end

        function params = cameras(this)
            %CAMERAS  Estimates camera parameters
            %
            %     params = obj.cameras()
            %
            % ## Output
            % * __params__ Describes camera parameters, a struct-array with
            %   the following fields:
            %   * __aspect__ Aspect ratio.
            %   * __focal__ Focal length.
            %   * __ppx__ Principal point X.
            %   * __ppy__ Principal point Y.
            %   * __R__ 3x3 camera rotation matrix.
            %   * __t__ 3x1 camera translation vector.
            %   * __K__ 3x3 camera intrinsic parameters.
            %
            % Note: Translation is assumed to be zero during the whole
            % stitching pipeline.
            %
            % See also: cv.Stitcher.component, cv.Stitcher.workScale
            %
            params = Stitcher_(this.id, 'cameras');
        end

        function wscale = workScale(this)
            %WORKSCALE  Work scale
            %
            %     wscale = obj.workScale()
            %
            % ## Output
            % * __wscale__ scalar double value.
            %
            % See also: cv.Stitcher.cameras, cv.Stitcher.component
            %
            wscale = Stitcher_(this.id, 'workScale');
        end
    end

    methods
        function mask = getMatchingMask(this)
            %GETMATCHINGMASK
            %
            %     mask = obj.getMatchingMask()
            %
            % ## Output
            % * __mask__
            %
            % See also: cv.Stitcher.setMatchingMask
            %
            mask = Stitcher_(this.id, 'getMatchingMask');
        end

        function setMatchingMask(this, mask)
            %SETMATCHINGMASK
            %
            %     obj.setMatchingMask(mask)
            %
            % ## Input
            % * __mask__
            %
            % See also: cv.Stitcher.getMatchingMask
            %
            Stitcher_(this.id, 'setMatchingMask', mask);
        end

        function value = getFeaturesFinder(this)
            %GETFEATURESFINDER  Get the features finder
            %
            %     value = obj.getFeaturesFinder()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setFeaturesFinder
            %
            value = Stitcher_(this.id, 'getFeaturesFinder');
        end

        function setFeaturesFinder(this, finderType, varargin)
            %SETFEATURESFINDER  Set the features finder
            %
            %     obj.setFeaturesFinder(finderType)
            %     obj.setFeaturesFinder(finderType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __finderType__ Feature finder type. One of:
            %   * __OrbFeaturesFinder__ ORB features finder. See cv.ORB
            %   * __AKAZEFeaturesFinder__ AKAZE features finder. See cv.AKAZE
            %   * __SurfFeaturesFinder__ SURF features finder. See cv.SURF
            %     (requires `xfeatures2d` module)
            %   * __SurfFeaturesFinderGpu__ (requires CUDA and `xfeatures2d`
            %     module)
            %
            % ## Options
            % The following are options for the various finders:
            %
            % ### `OrbFeaturesFinder`
            % * __GridSize__ default [3,1]
            % * __NFeatures__ default 1500
            % * __ScaleFactor__ default 1.3
            % * __NLevels__ default 5
            %
            % ### `AKAZEFeaturesFinder`
            % * __DescriptorType__ default 'MLDB'
            % * __DescriptorSize__ default 0
            % * __DescriptorChannels__ default 3
            % * __Threshold__ default 0.001
            % * __NOctaves__ default 4
            % * __NOctaveLayers__ default 4
            % * __Diffusivity__ default `PM_G2`
            %
            % ### `SurfFeaturesFinder`
            % * __HessThresh__ default 300.0
            % * __NumOctaves__ default 3
            % * __NumLayers__ default 4
            % * __NumOctaveDescr__ default 3
            % * __NumLayersDesc__ default 4
            %
            % The class uses `OrbFeaturesFinder` by default or
            % `SurfFeaturesFinder` if `xfeatures2d` module is available.
            %
            % See also: cv.Stitcher.getFeaturesFinder, cv.FeaturesMatcher
            %
            Stitcher_(this.id, 'setFeaturesFinder', finderType, varargin{:});
        end

        function value = getFeaturesMatcher(this)
            %GETFEATURESMATCHER  Get the features matcher
            %
            %     value = obj.getFeaturesMatcher()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setFeaturesMatcher
            %
            value = Stitcher_(this.id, 'getFeaturesMatcher');
        end

        function setFeaturesMatcher(this, matcherType, varargin)
            %SETFEATURESMATCHER  Set the features matcher
            %
            %     obj.setFeaturesMatcher(matcherType)
            %     obj.setFeaturesMatcher(matcherType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __matcherType__ Feature matcher type. One of:
            %   * __BestOf2NearestMatcher__ A "best of 2 nearest" matcher.
            %     Features matcher which finds two best matches for each
            %     feature and leaves the best one only if the ratio between
            %     descriptor distances is greater than the threshold
            %     `MatchConf`.
            %   * __BestOf2NearestRangeMatcher__
            %   * __AffineBestOf2NearestMatcher__ A "best of 2 nearest"
            %     matcher that expects affine trasformation between images.
            %     Features matcher similar to `BestOf2NearestMatcher` which
            %     finds two best matches for each feature and leaves the best
            %     one only if the ratio between descriptor distances is
            %     greater than the threshold `MatchConf`. Unlike
            %     `BestOf2NearestMatcher` this matcher uses affine
            %     transformation (affine trasformation estimate will be placed
            %     in `matches_info`).
            %
            % ## Options
            % The following are options accepted by all matchers:
            %
            % * __TryUseGPU__ Should try to use GPU or not. default false
            % * __MatchConf__ Match distances ration threshold. default 0.3
            % * __NumMatchesThresh1__ Minimum number of matches required for
            %   the 2D projective transform estimation used in the inliers
            %   classification step. default 6
            % * __NumMatchesThresh2__ Minimum number of matches required for
            %   the 2D projective transform re-estimation on inliers.
            %   default 6
            %
            % The following are options for the various algorithms:
            %
            % ### `BestOf2NearestRangeMatcher`
            % * __RangeWidth__ default 5
            %
            % ### `AffineBestOf2NearestMatcher`
            % * __FullAffine__ whether to use full affine transformation with
            %   6 degress of freedom or reduced transformation with 4 degrees
            %   of freedom using only rotation, translation and uniform
            %   scaling. default false
            %
            % The class uses `BestOf2NearestMatcher` by default.
            %
            % See also: cv.Stitcher.getFeaturesMatcher, cv.FeaturesMatcher
            %
            Stitcher_(this.id, 'setFeaturesMatcher', matcherType, varargin{:});
        end

        %{
        function value = getEstimator(this)
            %GETESTIMATOR  Get the estimator
            %
            %     value = obj.getEstimator()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setEstimator
            %
            value = Stitcher_(this.id, 'getEstimator');
        end

        function setEstimator(this, estimatorType, varargin)
            %SETESTIMATOR  Set the estimator
            %
            %     obj.setEstimator(estimatorType)
            %     obj.setEstimator(estimatorType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __estimatorType__ Estimator type. One of:
            %   * __HomographyBasedEstimator__ Homography based rotation
            %     estimator.
            %   * __AffineBasedEstimator__ Affine transformation based
            %     estimator. This estimator uses pairwise transformations
            %     estimated by matcher to estimate final transformation for
            %     each camera.
            %
            % The following are options for the various algorithms:
            %
            % ### `HomographyBasedEstimator`
            % * __IsFocalsEstimated__ default false
            %
            % See also: cv.Stitcher.getEstimator. cv.Estimator
            %
            Stitcher_(this.id, 'setEstimator', estimatorType, varargin{:});
        end
        %}

        function value = getBundleAdjuster(this)
            %GETBUNDLEADJUSTER  Get the bundle adjuster
            %
            %     value = obj.getBundleAdjuster()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setBundleAdjuster
            %
            value = Stitcher_(this.id, 'getBundleAdjuster');
        end

        function setBundleAdjuster(this, adjusterType, varargin)
            %SETBUNDLEADJUSTER  Set the bundle adjuster
            %
            %     obj.setBundleAdjuster(adjusterType)
            %     obj.setBundleAdjuster(adjusterType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __adjusterType__ camera parameters refinement method. One of:
            %   * __NoBundleAdjuster__ Stub bundle adjuster that does nothing.
            %   * __BundleAdjusterRay__ Implementation of the camera
            %     parameters refinement algorithm which minimizes sum of the
            %     distances between the rays passing through the camera center
            %     and a feature. It can estimate focal length. It ignores the
            %     refinement mask for now.
            %   * __BundleAdjusterReproj__ Implementation of the camera
            %     parameters refinement algorithm which minimizes sum of the
            %     reprojection error squares. It can estimate focal length,
            %     aspect ratio, principal point. You can affect only on them
            %     via the refinement mask.
            %   * __BundleAdjusterAffine__ Bundle adjuster that expects affine
            %     transformation represented in homogeneous coordinates in R
            %     for each camera param. Implements camera parameters
            %     refinement algorithm which minimizes sum of the reprojection
            %     error squares. It estimates all transformation parameters.
            %     Refinement mask is ignored. See also
            %     cv.AffineBasedEstimator, `AffineBestOf2NearestMatcher`.
            %   * __BundleAdjusterAffinePartial__ Bundle adjuster that expects
            %     affine transformation with 4 DOF represented in homogeneous
            %     coordinates in R for each camera param. Implements camera
            %     parameters refinement algorithm which minimizes sum of the
            %     reprojection error squares. It estimates all transformation
            %     parameters. Refinement mask is ignored.
            %
            % ## Options
            % The following are options accepted by all adjusters:
            %
            % * __ConfThresh__ default 1
            % * __RefinementMask__ default `eye(3)`
            % * __TermCriteria__ default
            %   `struct('type','Count+EPS', 'maxCount',1000, 'epsilon',eps)`
            %
            % The class uses `BundleAdjusterRay` by default.
            %
            % See also: cv.Stitcher.getBundleAdjuster, cv.BundleAdjuster
            %
            Stitcher_(this.id, 'setBundleAdjuster', adjusterType, varargin{:});
        end

        function value = getWarper(this)
            %GETWARPER  Get the warper
            %
            %     value = obj.getWarper()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setWarper
            %
            value = Stitcher_(this.id, 'getWarper');
        end

        function setWarper(this, warperType, varargin)
            %SETWARPER  Set the image warper
            %
            %     obj.setWarper(warperType)
            %     obj.setWarper(warperType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __warperType__ image warper factory class type, used to create
            %   the rotation-based warper. One of:
            %   * __PlaneWarper__ Plane warper factory class. Warper that maps
            %     an image onto the `z = 1` plane.
            %   * __PlaneWarperGpu__ (requires CUDA)
            %   * __AffineWarper__ Affine warper factory class. Affine warper
            %     that uses rotations and translations. Uses affine
            %     transformation in homogeneous coordinates to represent both
            %     rotation and translation in camera rotation matrix.
            %   * __CylindricalWarper__ Cylindrical warper factory class.
            %     Warper that maps an image onto the `x*x + z*z = 1` cylinder.
            %   * __CylindricalWarperGpu__ (requires CUDA)
            %   * __SphericalWarper__ Warper that maps an image onto the unit
            %     sphere located at the origin. Projects image onto unit
            %     sphere with origin at [0,0,0] and radius `scale`, measured
            %     in pixels. A 360 panorama would therefore have a resulting
            %     width of `2*scale*pi` pixels. Poles are located at [0,-1,0]
            %     and [0,1,0] points.
            %   * __SphericalWarperGpu__ (requires CUDA)
            %   * __FisheyeWarper__
            %   * __StereographicWarper__
            %   * __CompressedRectilinearWarper__
            %   * __CompressedRectilinearPortraitWarper__
            %   * __PaniniWarper__
            %   * __PaniniPortraitWarper__
            %   * __MercatorWarper__
            %   * __TransverseMercatorWarper__
            %
            % ## Options
            % The following are options for the various warpers:
            %
            % ### `CompressedRectilinearWarper`, `CompressedRectilinearPortraitWarper`, `PaniniWarper`, `PaniniPortraitWarper`
            % * __A__ default 1
            % * __B__ default 1
            %
            % The class uses `SphericalWarper` by default.
            %
            % See also: cv.Stitcher.getWarper, cv.RotationWarper
            %
            Stitcher_(this.id, 'setWarper', warperType, varargin{:});
        end

        function value = getExposureCompensator(this)
            %GETEXPOSURECOMPENSATOR  Get the exposire compensator
            %
            %     value = obj.getExposureCompensator()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setExposureCompensator
            %
            value = Stitcher_(this.id, 'getExposureCompensator');
        end

        function setExposureCompensator(this, compensatorType, varargin)
            %SETEXPOSURECOMPENSATOR  Set the exposure compensator
            %
            %     obj.setExposureCompensator(compensatorType)
            %     obj.setExposureCompensator(compensatorType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __compensatorType__ exposure compensator type. One of:
            %   * __NoExposureCompensator__ Stub exposure compensator which
            %     does nothing.
            %   * __GainCompensator__ Exposure compensator which tries to
            %     remove exposure related artifacts by adjusting image
            %     intensities, see [BL07] and [WJ10] for details.
            %   * __BlocksGainCompensator__ Exposure compensator which tries
            %     to remove exposure related artifacts by adjusting image
            %     block intensities, see [UES01] for details.
            %
            % ## Options
            % The following are options for the various compensators:
            %
            % ### `BlocksGainCompensator`
            % * __Width__ Block width. default 32
            % * __Heigth__ Block height. default 32
            %
            % The class uses `BlocksGainCompensator` by default.
            %
            % ## References
            % [BL07]:
            % > Matthew Brown and David G Lowe.
            % > "Automatic panoramic image stitching using invariant features".
            % > International journal of computer vision, 74(1):59-73, 2007.
            %
            % [WJ10]:
            % > Wei Xu and Jane Mulligan. "Performance evaluation of color
            % > correction approaches for automatic multi-view image and video
            % > stitching". In Computer Vision and Pattern Recognition (CVPR),
            % > 2010 IEEE Conference on, pages 263-270. IEEE, 2010.
            %
            % [UES01]:
            % > Matthew Uyttendaele, Ashley Eden, and R Skeliski.
            % > "Eliminating ghosting and exposure artifacts in image mosaics".
            % > In Computer Vision and Pattern Recognition, 2001. CVPR 2001.
            % > Proceedings of the 2001 IEEE Computer Society Conference on,
            % > volume 2, pages II-509. IEEE, 2001.
            %
            % See also: cv.Stitcher.getExposureCompensator,
            %  cv.ExposureCompensator
            %
            Stitcher_(this.id, 'setExposureCompensator', compensatorType, varargin{:});
        end

        function value = getSeamFinder(this)
            %GETSEAMFINDER  Get the seam finder
            %
            %     value = obj.getSeamFinder()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setSeamFinder
            %
            value = Stitcher_(this.id, 'getSeamFinder');
        end

        function setSeamFinder(this, seamType, varargin)
            %SETSEAMFINDER  Set the seam finder
            %
            %     obj.setSeamFinder(seamType)
            %     obj.setSeamFinder(seamType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __seamType__ seam estimator type. One of:
            %   * __NoSeamFinder__ Stub seam estimator which does nothing.
            %   * __VoronoiSeamFinder__ Voronoi diagram-based pairwise seam
            %     estimator.
            %   * __DpSeamFinder__
            %   * __GraphCutSeamFinder__ Minimum graph cut-based seam
            %     estimator. See details in [V03].
            %   * __GraphCutSeamFinderGpu__ (requires CUDA)
            %
            % ## Options
            % The following are options for the various seam finders:
            %
            % ### `DpSeamFinder`
            % * __CostFunction__ default 'Color'. One of:
            %   * __Color__
            %   * __ColorGrad__
            %
            % ### `GraphCutSeamFinder`
            % * __CostType__ default 'ColorGrad'. One of:
            %   * __Color__
            %   * __ColorGrad__
            % * __TerminalCost__ default 10000.0
            % * __BadRegionPenaly__ default 1000.0
            %
            % The class uses `GraphCutSeamFinder` by default.
            %
            % ## References
            % [V03]:
            % > Vivek Kwatra, Arno Schodl, Irfan Essa, Greg Turk, and Aaron
            % > Bobick. "Graphcut textures: image and video synthesis using
            % > graph cuts". In ACM Transactions on Graphics (ToG), volume 22,
            % > pages 277-286. ACM, 2003.
            %
            % See also: cv.Stitcher.getSeamFinder, cv.SeamFinder
            %
            Stitcher_(this.id, 'setSeamFinder', seamType, varargin{:});
        end

        function value = getBlender(this)
            %GETBLENDER  Get the blender
            %
            %     value = obj.getBlender()
            %
            % ## Output
            % * __value__ output scalar struct.
            %
            % See also: cv.Stitcher.setBlender
            %
            value = Stitcher_(this.id, 'getBlender');
        end

        function setBlender(this, blenderType, varargin)
            %SETBLENDER  Set the blender
            %
            %     obj.setBlender(blenderType)
            %     obj.setBlender(blenderType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __blenderType__ image blender type. One of:
            %   * __NoBlender__ Simple blender which puts one image over
            %     another.
            %   * __FeatherBlender__ Simple blender which mixes images at its
            %     borders.
            %   * __MultiBandBlender__ Blender which uses multi-band blending
            %     algorithm (see [BA83]).
            %
            % ## Options
            % The following are options for the various blenders:
            %
            % ### `FeatherBlender`
            % * __Sharpness__ default 0.02
            %
            % ### `MultiBandBlender`
            % * __TryGPU__ default false
            % * __NumBands__ default 5
            % * __WeightType__ One of:
            %   * __single__ (default)
            %   * __int16__
            %
            % The class uses `MultiBandBlender` by default.
            %
            % ## References
            % [BA83]:
            % > Peter J Burt and Edward H Adelson.
            % > "A multiresolution spline with application to image mosaics".
            % > ACM Transactions on Graphics (TOG), 2(4):217-236, 1983.
            %
            % See also: cv.Stitcher.getBlender, cv.Blender
            %
            Stitcher_(this.id, 'setBlender', blenderType, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.RegistrationResol(this)
            value = Stitcher_(this.id, 'get', 'RegistrationResol');
        end
        function set.RegistrationResol(this, value)
            Stitcher_(this.id, 'set', 'RegistrationResol', value);
        end

        function value = get.SeamEstimationResol(this)
            value = Stitcher_(this.id, 'get', 'SeamEstimationResol');
        end
        function set.SeamEstimationResol(this, value)
            Stitcher_(this.id, 'set', 'SeamEstimationResol', value);
        end

        function value = get.CompositingResol(this)
            value = Stitcher_(this.id, 'get', 'CompositingResol');
        end
        function set.CompositingResol(this, value)
            Stitcher_(this.id, 'set', 'CompositingResol', value);
        end

        function value = get.PanoConfidenceThresh(this)
            value = Stitcher_(this.id, 'get', 'PanoConfidenceThresh');
        end
        function set.PanoConfidenceThresh(this, value)
            Stitcher_(this.id, 'set', 'PanoConfidenceThresh', value);
        end

        function value = get.WaveCorrection(this)
            value = Stitcher_(this.id, 'get', 'WaveCorrection');
        end
        function set.WaveCorrection(this, value)
            Stitcher_(this.id, 'set', 'WaveCorrection', value);
        end

        function value = get.WaveCorrectKind(this)
            value = Stitcher_(this.id, 'get', 'WaveCorrectKind');
        end
        function set.WaveCorrectKind(this, value)
            Stitcher_(this.id, 'set', 'WaveCorrectKind', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function log_b = binomial_coeff(n_max)
    % Compute log binomial coefficients
    log_b = {};
    for n = 0:n_max
        Ks = 0:n;
        log_b{n+1} = gammaln(n+1) - gammaln(Ks+1) - gammaln(n-Ks+1);
    end
end


text:
classdef KTUEAMapFile

% UEAMAPFILE -- Defines a class that contains information on a UEA
% electrode array using CMP files provided by Blackrock Microsystems.
%
% To load a CMP file type 'MapName = UEAMapFile' where MapName is the name
% of the variable that will be created and will contain the map class.
%
% Example: myArrayMap = UEAMapFile;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% METHODS:
%
% There are several available methods for mapfiles. Here's a list and a
% description of what they perform.
%
% Electrode2Channel(Electrode):
%   Will return the Channel number corresponding to a passed Electrode.
%
% Channel2Electrode(Channel):
%   Will return the Electrode number corresponding to a passed Channel.
%
% GetChannelBankID(Channel):
%   Returns the BankID for a passed Channel.
%
% GetChannelPin(Channel):
%   Returns the Pin for a passed Channel.
%
% GetChannelLabel(Channel):
%   Returns the Label for a passed Channel.
%
% GenerateChannelSubplot(Channel):
%   Returns a subplot handle for the passed Channel to be used in plotting
%   UEA-like maps.
%
% GenerateChannelSubplotNames(Channel):
%   Plots the names of channels and electrodes on the subplot corresponding
%   to the passed Channel.
%
% GetChannelColumnRow(Channel):
%   Returns the Column and Row positions of the passed number that is used
%   to plot UEA-like maps.
%
% PlotCMP:
%   Will plot a map of the CMP Map with all the electrode and channel
%   names.
%
%   Kian Torab
%   Blackrock Microsystems
%   ktorab@blackrockmicro.com
%
%   Version 1.6.1.0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
% 1.6.1.0:
%   - Minor bug fix.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
    properties (SetAccess = private)
        ElecNum
        ChanNum
        Column
        Row
        Bank
        Pin
        Label
        pathHandle
        fileHandle
    end
    methods (Hidden)
        function obj = KTUEAMapFile(fileName)
            if ~exist('fileName', 'var')
                if exist('getFile.m', 'file') == 2
                    [obj.fileHandle obj.pathHandle] = getFile('*.cmp', 'Open a CMP Mapfile...');
                else
                    [obj.fileHandle obj.pathHandle] = uigetfile('*.cmp', 'Open a CMP Mapfile...');
                end
                if ~obj.fileHandle
                    disp('No file was selected.');
                    return;
                end
                mapfileDataCell = importdata([obj.pathHandle obj.fileHandle], ' ', 200);
            else
                mapfileDataCell = importdata(fileName, ' ', 200);
            end
            mapfileDataCellParsed = regexp(mapfileDataCell, '\t', 'split');
            if size(mapfileDataCellParsed, 2) == 1
                mapfileDataCellParsed = regexp(mapfileDataCell, '\s+', 'split');
            end
            lastJunkIDX = 1;
            while 1
                if size(mapfileDataCellParsed{lastJunkIDX}, 2) == 6 || size(mapfileDataCellParsed{lastJunkIDX}, 2) == 5
                    if ~isempty(str2num(mapfileDataCellParsed{lastJunkIDX}{1})) && ...
                    ~isempty(str2num(mapfileDataCellParsed{lastJunkIDX}{2})) && ...
                    ~isempty(str2num(mapfileDataCellParsed{lastJunkIDX}{4}))
                    lastJunkIDX = lastJunkIDX - 1;
                        break;
                    end
                end
                lastJunkIDX = lastJunkIDX + 1;
            end
            mapfileDataCellParsed(1:lastJunkIDX) = [];
            for i = 1:size(mapfileDataCellParsed, 1)
                if ~all(isspace(mapfileDataCellParsed{i}))
                    obj.Column(i) = str2num(mapfileDataCellParsed{i,:}{1});
                    obj.Row(i)    = str2num(mapfileDataCellParsed{i,:}{2});
                    obj.Bank(i)   = mapfileDataCellParsed{i,:}{3}-'@';
                    obj.Pin(i)    = str2num(mapfileDataCellParsed{i,:}{4});
                    ElectNum      = str2num(mapfileDataCellParsed{i,:}{5}(5:end));
                    if isempty(ElectNum)
                        obj.ElecNum(i) = str2num(mapfileDataCellParsed{i,:}{5}(2:end-4));
                    else
                        obj.ElecNum(i) = ElectNum;
                    end
                    obj.ChanNum(i)  = (obj.Bank(i) - 1) * 32 + obj.Pin(i);
                    obj.Label{i}    = mapfileDataCellParsed{i,:}{5};
                end
            end
        end
    end
    methods
        function validFlag = isValid(obj)
            if ~obj.getFilename
                validFlag = 0;
            else
                validFlag = 1;
            end
        end
        function FilePath = getPathName(obj)
            FilePath = obj.pathHandle;
        end
        function FileName = getFilename(obj)
            FileName = obj.fileHandle;
        end
        function Channel = Electrode2Channel(obj, Electrode)
            if ~exist('Electrode', 'var'); disp('Electrode is a required input. Refer to help for more information.'); return; end;
            if isempty(obj.ChanNum(obj.ElecNum == Electrode(1)))
                disp('Electrode number is invalid.');
                Channel = [];
                return;
            end
            for ElecIDX = 1:length(Electrode)
                Channel(ElecIDX) = obj.ChanNum(obj.ElecNum == Electrode(ElecIDX));
            end
        end
        function Electrode = Channel2Electrode(obj, Channel)
            if isempty(obj.ElecNum(obj.ChanNum == Channel(1)))
                disp('Channel number is invalid.');
                Electrode = [];
                return;
            end
            if ~exist('Channel', 'var'); disp('Channel is a required input. Refer to help for more information.'); return; end;
            for ChanIDX = 1:length(Channel)
                Electrode(ChanIDX) = obj.ElecNum(obj.ChanNum == Channel(ChanIDX));
            end
        end
        function BankID = getChannelBankID(obj, Channel)
            BankID = obj.Bank(obj.ChanNum == Channel);
        end
        function Label = getChannelLabel(obj, Channel)
            Label = char(obj.Label{obj.ChanNum == Channel});
        end
        function Pin = getChannelPin(obj, Channel)
            Pin = obj.Pin(obj.ChanNum == Channel);
        end
        function spHandle = GenerateChannelSubplot(obj, Channel)
            if ~exist('Channel', 'var'); disp('Channel is a required input. Refer to help for more information.'); return; end;
            [x, y] = getChannelColumnRow(obj, Channel);
            spHandle=subplot('position', [x/10+0.004, y/10+0.004, 1/10-0.004, 1/10-0.004]);
            axis off;
            box on;
            set(spHandle,'XTickLabel', []);
            set(spHandle,'YTickLabel', []);
            %%
        end
        function GenerateChannelSubplotNames(obj, Channel, fColor)
            if ~exist('Channel', 'var'); disp('Channel is a required input. Refer to help for more information.'); return; end;
            if ~exist('fColor', 'var'); fColor  = [1,1,1]; end
            Electrode = Channel2Electrode(obj, Channel);
            text(0.13, 0.23, ['elec ', num2str(Electrode)], 'Color', fColor, 'FontSize', 8);
            text(0.13, 0.10, ['chan ', num2str(Channel)], 'Color', fColor, 'FontSize', 8);
        end
        function [x, y] = getChannelColumnRow(obj, Channel)
            if ~exist('Channel', 'var'); disp('Channel is a required input. Refer to help for more information.'); return; end;
            x = obj.Column(obj.ChanNum == Channel);
            y = obj.Row(obj.ChanNum == Channel);
        end
        function PlotCMP(obj)
            figure;
            for Channel = 1:96
                GenerateChannelSubplot(obj, Channel);
                GenerateChannelSubplotNames(obj, Channel, [0,0,0]);
            end
        end
        function setAxisBackgroundColor(~, curColor)
            if isa(curColor, 'double') && (length(curColor) == 3)
                if all(curColor <= 1) && all(curColor >= 0)
                    set(gca, 'Color', curColor);
                else
                    disp('The color values have to be between 0 and 1.');
                end
            elseif isa(curColor, 'char')
                try
                	set(gca, 'Color', curColor);
                catch
                    disp('Color is not valid. Type ''doc color'' for more information.');
                end
            else
                disp('The input needs to be a valid color input.');
                disp('A valid color input is a vector of length 3 of values between 0 and 1.');
            end
        end
        function setXAxisColor(~, curColor)
            if isa(curColor, 'double') && (length(curColor) == 3)
                if all(curColor <= 1) && all(curColor >= 0)
                    set(gca, 'XColor', curColor);
                else
                    disp('The color values have to be between 0 and 1.');
                end
            elseif isa(curColor, 'char')
                try
                    set(gca, 'XColor', curColor);
                catch
                    disp('Color is not valid. Type ''doc color'' for more information.');
                end
            else
                disp('The input needs to be a valid color input.');
                disp('A valid color input is a vector of length 3 of values between 0 and 1.');
            end
        end
        function setYAxisColor(~, curColor)
            if isa(curColor, 'double') && (length(curColor) == 3)
                if all(curColor <= 1) && all(curColor >= 0)
                    set(gca, 'YColor', curColor);
                else
                    disp('The color values have to be between 0 and 1.');
                end
            elseif isa(curColor, 'char')
                try
                    set(gca, 'YColor', curColor);
                catch
                    disp('Color is not valid. Type ''doc color'' for more information.');
                end
            else
                disp('The input needs to be a valid color input.');
                disp('A valid color input is a vector of length 3 of values between 0 and 1.');
            end
        end
        function setAxesColor(~, curColor)
            if isa(curColor, 'double') && (length(curColor) == 3)
                if all(curColor <= 1) && all(curColor >= 0)
                    set(gca, 'YColor', curColor);
                    set(gca, 'XColor', curColor);
                else
                    disp('The color values have to be between 0 and 1.');
                end
            elseif isa(curColor, 'char')
                try
                    set(gca, 'YColor', curColor);
                    set(gca, 'XColor', curColor);
                catch
                    disp('Color is not valid. Type ''doc color'' for more information.');
                end
            else
                disp('The input needs to be a valid color input.');
                disp('A valid color input is a vector of length 3 of values between 0 and 1.');
            end
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
addpath tensorIO_matlab

clear

clf

tensors = readTensors('outTSDF.tensor');

tsdf = tensors(1).value;

%% visualization
%{

for i=1:size(tsdf,3)
    imagesc((tsdf(end:-1:1,end:-1:1,i))',[0 1]); axis equal; axis tight; colorbar
    title(i);
    pause(0.1);
end

%return;

for i=1:size(tsdf,2)
    imagesc((reshape(tsdf(:,i,:), [size(tsdf,1) size(tsdf,3)])),[0 1]); axis equal; axis tight; colorbar
    pause(0.1);
    title(i);
end

for i=1:size(tsdf,1)
    imagesc((reshape(tsdf(i,end:-1:1,:), [size(tsdf,2) size(tsdf,3)])),[0 1]); axis equal; axis tight; colorbar
    title(i);
    pause(0.1);
end
%}

%% meshing
disp('isosurfacing...');
tic;
fv = isosurface(tsdf,0);
toc;
% visualizaiton
%{
figure(3);
p = patch(fv);
p.FaceColor = 'red';
p.EdgeColor = 'none';
daspect([1,1,1])
view(3); axis tight
camlight
lighting gouraud
%}

unit = 0.0005;

fv.vertices = fv.vertices * unit;

fv.vertices(:,1) = fv.vertices(:,1) - mean(fv.vertices(:,1));
fv.vertices(:,2) = fv.vertices(:,2) - mean(fv.vertices(:,2));
fv.vertices(:,3) = fv.vertices(:,3) - mean(fv.vertices(:,3));

mesh2off('isosurface.off', fv.faces,fv.vertices);

author: jianxiongxiao

repo: ProfXkit


text:
function F = coeffs2diskfun(X)
%COEFFS2DISKFUN   Convert a matrix of Chebyshev-Fourier coefficients to a
%                 diskfun.
%
%   F = coeffs2diskfun( X ) returns a diskfun object F that has a
%   Chebyshev-Fourier matrix of coefficients X.  This is useful for
%   computing quantities on the disk with the function F.
%
% See also DISKFUN/COEFFS2

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% To get the correct doubled-up grid, we can only accept
% certain coeff matrix sizes:

% Get size:
[m, n] = size( X );

% If n is odd, then make it even.
if ( mod(n, 2) == 1 )
    X = [ zeros(m, 1) X ];
end

%If m is even, we need it to be odd:
if (mod(m, 2) == 0)
    m=m+1;
    X = chebtech2.alias(X, m);
end

% Convert to values on the grid.
vals = trigtech.coeffs2vals(chebtech2.coeffs2vals(X).').';

% Assume that the function is real-valued.
vals = real( vals );

% Restrict to the region of interest.
vals = vals(floor(m/2)+1:m, :);

% Finally, make a diskfun object out of the values.
F = diskfun( vals );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c = camGraphics(csize,euler)

%CAMGRAPHICS Create a camera graphics 3D object.
%   CAMGRAPHICS creates a solid representing a camera. The objective is
%   heading towards the camera optical axis and is located at camera
%   position. Initial configuration is: position at origin and optical
%   axis aligned with world's Z axis.
%
%   The result is a structure with fields:
%       .vert0  the vertices in camera frame.
%       .vert   the vertices in world frame.
%       .faces  the definition of faces.
%
%   Fields .vert and .faces are used to draw the object via the PATCH
%   command. Further object repositionning is accomplished with DRAWOBJECT.
%
%   CAMGRAPHICS(SIZE) allows for choosing the camera size. Default is 0.1.
%
%   CAMGRAPHICS(SIZE,EULER) admits a different orientation to the default
%   one. To make the camera look in the X-axis direction, use EULER =
%   [-pi/2 0 -pi/2].
%
%   See also PATCH, SET, DRAWOBJECT.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargin < 2
    euler = zeros(3,1);
    if nargin < 1
        csize = 0.1;
    end
end

d = .6; % objective diameter
f = .5; % objective length
h = .8; % body height
w = 1;  % body width
t = .3; % body thickness

% objective vertices
a         = (0:pi/4:2*pi-pi/4)';
na        = length(a);
circle    = d/2*[cos(a),sin(a)];
objective = [circle,zeros(na,1);circle,f*ones(na,1)];

% body vertices
rect      = [1 1;-1 1;-1 -1;1 -1]*diag([w h])/2;
body      = [rect,zeros(4,1);rect,-t*ones(4,1)];

% rotation matrix
R = e2R(euler);

% vertices in camera frame
vert0     = [body;objective]*R';

% graphics structure - vertices and faces
c.vert0   = csize*vert0;
c.vert    = c.vert0;
c.faces   = [ ...
    01 02 03 04
    05 06 07 08
    01 02 06 05
    02 03 07 06
    03 04 08 07
    04 01 05 08
    09 10 18 17
    10 11 19 18
    11 12 20 19
    12 13 21 20
    13 14 22 21
    14 15 23 22
    15 16 24 23
    16 09 17 24];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [M1] = spm_eeg_inv_rigidreg(data1, data2)
% Computes homogenious transformation matrix based on two sets
% of points from two coordinate systems
%
% FORMAT [M1] = spm_eeg_inv_rigidreg(data1, data2)
% Input:
% data1      - locations of the first set of points corresponding to the
%            3D surface to register onto
% data2      - locations of the second set of points corresponding to the
%            second 3D surface to be registered
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_inv_rigidreg.m 5219 2013-01-29 17:07:07Z spm $

M       = spm_detrend(data1');
S       = spm_detrend(data2');
[U,A,V] = svd(S'*M);
R1      = V*U';
if det(R1) < 0
    B      = eye(3);
    B(3,3) = det(V*U');
    R1     = V*B*U';
end
t1 = mean(data1,2) - R1*mean(data2,2);
M1 = [R1 t1; 0 0 0 1];

author: spm


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
classdef ProbabilityDensity2 < handle
    properties(Constant)
        DEFAULT_BAND_WIDTH=.014;
        DEFAULT_GRID_SIZE=256;
        N_MIN=5000;
    end

    properties(SetAccess=private)
        deltas;
        onScale;
        xyData;
        xgrid;
        ygrid;
        wmat;
        fmat;
        fmatVector;
        ye;
        xm;
        ym;
        h;
        M;
        N;
        D;
        mins;
        maxs;
        dataBinIdxs;
    end

    methods
        function this=ProbabilityDensity2(xyData, mins, maxs, gridSize, bandWidth)
            assert(nargin>0);
            [this.N, this.D]=size(xyData);
            assert(this.D==2); %2D only
            if nargin<5
                bandWidth=this.DEFAULT_BAND_WIDTH;
                if nargin<4
                    gridSize=this.DEFAULT_GRID_SIZE;
                    if nargin<3
                        maxs=[];
                        if nargin<2
                            mins=[];
                        end
                    end
                end
            end
            needsScaling=~isempty(maxs) || ~isempty(mins);
            if isempty(maxs)
                maxs=max(xyData);
            end
            if isempty(mins)
                mins=min(xyData);
            end
            if needsScaling
                this.onScale=MatBasics.FindOnScale(xyData, mins, maxs);
                cntEdge=size(xyData, 1)-sum(this.onScale);
                if cntEdge>0
                    xyData=xyData(this.onScale, :);
                    this.N=this.N-cntEdge;
                end

            end
            this.mins=mins;
            this.maxs=maxs;
            this.M =gridSize;
            this.deltas=1/(this.M-1)*(maxs-mins);
            this.xyData=xyData;
            this.h=zeros(1, this.D);
            for i =1:this.D
                if this.N<ProbabilityDensity2.N_MIN
                    this.h(i)=(this.N/ProbabilityDensity2.N_MIN)^(-1/6)*1.7...
                        *bandWidth*(maxs(i)-mins(i));
                else
                    this.h(i)=bandWidth*(maxs(i)-mins(i));
                end
            end
            this.computeWeight;
            this.computeDensity;
        end

        function drawJetColors(this, ax, numberOfJetColors)
            if nargin<3
                numberOfJetColors=128;
            end
            this.drawColors(ax, numberOfJetColors);
        end

        function drawContourOutliers(this, ax, hContours)
            levelList=get(hContours,'LevelList');
            outDensity=find(this.fmatVector<=levelList(1));
            outEvents=ismember(this.dataBinIdxs, outDensity);
            if isempty(this.onScale)
                outsideLastContour=outEvents;
            else
                outsideLastContour=false(length(this.onScale),1);
                outsideLastContour(this.onScale)=outEvents;
            end
            plot(ax, this.xyData(outsideLastContour, 1), ...
                this.xyData(outsideLastContour,2), ...
                '.', 'MarkerSize', 4, ...
                'MarkerEdgeColor', 'black',...
                'LineStyle', 'none');
        end

        function H=drawContours(this, ax, percent, color, lineWidth)
            if nargin<5
                lineWidth=1;
                if nargin<4
                    color=[.5 .5 .6];
                    if nargin<3
                        percent=10;
                    end
                end
            end
            numLevels=floor(100/percent);
            levels=this.computeLevels(numLevels);
            [~,H]=contour(ax, this.xm, this.ym, this.fmat, levels, ...
                'k', 'color', color, 'LineStyle', '-', 'LineWidth', lineWidth);
        end

        function drawColors(this, ax, numberOfJetColors, subsetOfData, ...
            colorRangeStart, colorRangeEnd)
            wasHeld=ishold(ax);
            if ~wasHeld
                hold(ax);
            end
            if nargin<5
                isJetColor=true;
                if nargin<4
                    subsetOfData=[];
                    if  nargin<3
                        numberOfJetColors=64;
                    end
                end
            else
                isJetColor=false;
            end
            if isempty(subsetOfData)
                data=this.xyData;
            else
                if size(subsetOfData, 1)==length(this.onScale)
                    data=subsetOfData(this.onScale, :);
                else
                    onScale_=MatBasics.FindOnScale(subsetOfData, ...
                        this.mins, this.maxs);
                    data=subsetOfData(onScale_, :);
                end
            end
            z=reshape(1:this.M^2,this.M,this.M);
            eb=interp2(this.xgrid, this.ygrid, z',...
                data(:,1),data(:,2),'nearest');  %this associates each data point with its nearest grid point
            [x1,~,x3]=unique(eb);
            if isJetColor
                colors=jet(numberOfJetColors);
            else
                if nargin<6
                    colors=colorRangeStart;
                else
                    a1=mean(colorRangeEnd);
                    a2=mean(colorRangeStart);
                    if a1<a2
                        m1=max(colorRangeEnd);
                        if m1>.75
                            f=.75/m1;
                            colorRangeEnd=colorRangeEnd*f;
                        end
                    end
                    nColors=32;
                    colors=zeros(nColors,3);
                    colors(1,:)=colorRangeStart;
                    colors(nColors,:)=colorRangeEnd;
                    gap=zeros(1,3);
                    for i=1:3
                        gap(i)=colorRangeEnd(i)-colorRangeStart(i);
                    end
                    for i=2:nColors-1
                        for j=1:3
                            colors(i,j)=colors(1,j)+(i/nColors*gap(1,j));
                        end
                    end
                end
            end
            nColors=length(colors);
            try
                levels=this.computeLevels(nColors);
            catch ex
            end
            if size(data,1)<10
                color=colors(1, :);
                plot(ax, data(:,1), data(:,2), 'd',...
                    'markersize', 2, 'MarkerEdgeColor',...
                    color, 'LineStyle', 'none');
                if ~wasHeld
                    hold(ax);
                end
                return;
            end
            try
                colormap(colors);
            catch ex
                disp('huh');
            end
            densities=this.fmatVector(x1);
            lookup=bsearch(levels,densities);
            eventColors=lookup(x3);
            usedColors=unique(eventColors);
            N2=length(usedColors);
            sz=size(data,1);
            if sz<10000
                marker='d';
                ms=2;
            else
                marker='.';
                ms=2;
            end
            for i=1:N2
                colorIdx=usedColors(i);
                li=eventColors==colorIdx;
                plot(ax, data(li,1), data(li,2), marker,...
                    'markersize', ms, 'MarkerEdgeColor',...
                    colors(colorIdx, :), 'LineStyle', 'none');
            end
            if ~wasHeld
                hold(ax);
            end
        end
    end

    methods(Access=private)

        function computeWeight(this)
            this.ye=zeros(2, this.M);
            pointLL=zeros(this.N,2);  %this will be the "lower left" gridpoint to each data point
            for ii = 1:2
                this.ye(ii,:) = linspace(this.mins(ii), this.maxs(ii), this.M);
                pointLL(:,ii)=floor((this.xyData(:,ii)-this.mins(ii))./this.deltas(ii)) + 1;
            end
            pointLL(pointLL==this.M)=this.M-1;  %this avoids going over grid boundary
            %% assign each data point to its closest grid point
            [this.xgrid, this.ygrid]=meshgrid(this.ye(1,:),this.ye(2,:));
            z=reshape(1:this.M^2, this.M, this.M);
            this.dataBinIdxs=interp2(this.xgrid, this.ygrid,z',...
                this.xyData(:,1),this.xyData(:,2),'nearest');  %this associates each data point with its nearest grid point

            %% compute w
            Deltmat=repmat(this.deltas, this.N,1);
            shape=this.M*ones(1,2);
            wmat_=zeros(this.M, this.M);
            for ii=0:1  %number of neighboring gridpoints in 2 dimensions
                for j=0:1
                    pointm=pointLL+repmat([j ii],this.N,1);  %indices of ith neighboring gridpoints
                    pointy=zeros(this.N,2);
                    for k=1:2
                        pointy(:,k)=this.ye(k,pointm(:,k));  %y-values of ith neighboring gridpoints
                    end
                    W=prod(1-(abs(this.xyData-pointy)./Deltmat),2);  %contribution to w from ith neighboring gridpoint from each datapoint
                    wmat_=wmat_+accumarray(pointm,W,shape);  %sums contributions for ith gridpoint over data points and adds to wmat
                end
            end
            this.wmat=wmat_;

        end

        function computeDensity(this)
            Z_=zeros(1, this.D);
            Zin=cell(1, this.D);
            for i =1:this.D
                Z_(i)=min(floor(4*this.h(i)/this.deltas(i)), this.M-1);
                Zin{i}=-Z_(i):Z_(i);
            end
            phi = @(x) 1/sqrt(2*pi)*exp(-x.^2./2);
            [L_{1},L_{2}]=meshgrid(Zin{1},Zin{2});
            Phix=phi(L_{1}*this.deltas(1)./this.h(1))./this.h(1);
            Phiy=phi(L_{2}*this.deltas(2)./this.h(2))./this.h(2);
            Phimat = (Phix.*Phiy)';
            fMat = 1/this.N*conv2(this.wmat,Phimat,'same');
            this.fmatVector=reshape(fMat,[1, this.M^2]);
            this.fmat=fMat';
            [this.xm, this.ym]=meshgrid(this.ye(1,:),this.ye(2,:));
        end

        function levels=computeLevels(this, numberOfLevels)
            T=sort(reshape(this.fmat, 1, this.M^this.D));
            CT=cumsum(T);
            NT=CT/CT(end);
            levels=zeros(1, numberOfLevels);
            for level=1:numberOfLevels
                idx=bsearch(NT, level/numberOfLevels);
                levels(level)=T(idx);
            end
        end

    end

    methods(Static)
        function Draw(ax, data, doContours, doJetColors, reset, ...
                stretchPerc, contourPerc)
            if nargin<7
                contourPerc=10;
                if nargin<6
                    stretchPerc=0;
                    if nargin<5
                        reset=true;
                        if nargin<4
                            doJetColors=true;
                            if nargin<3
                                doContours=true;
                            end
                        end
                    end
                end
            end
            if reset
                cla(ax, 'reset');
            end
            try
                if stretchPerc>0
                    pb=ProbabilityDensity2(MatBasics.StretchXyData(...
                        data, stretchPerc));
                else
                    pb=ProbabilityDensity2(data);
                end
                wasHeld=ishold(ax);
                if ~wasHeld
                    hold(ax, 'on');
                end
                if doJetColors
                    pb.drawJetColors(ax);
                end
                if doContours
                    hContour=pb.drawContours(ax, contourPerc);
                    if ~doJetColors
                        pb.drawContourOutliers(ax, hContour);
                    end
                end
                if ~wasHeld
                    hold(ax, 'off');
                end
            catch ex
                ex.getReport
                plot(ax, data(:,1), data(:,2));
            end
        end

        function [fncMotion, javaLegend, btns, btnLbls, plots, labelHs]=DrawLabeled(ax, ...
                data, lbls, lblMap, doContours, reset, priorFcn, doubleClicker, ...
                xMargin, yMargin, doJavaLegend, oldJavaBtns, southComponent, ...
                selectedCallback, stretchPerc)
            fncMotion=[];
            javaLegend=[];
            btns=[];
            btnLbls=[];
            plots=[];
            if nargin<15
                stretchPerc=0;
                if nargin<14
                    selectedCallback=[];
                    if nargin<13
                        southComponent=[];
                        if nargin<12
                            oldJavaBtns=[];
                            if nargin<11
                                doJavaLegend=false;
                                if nargin<10
                                    yMargin=-0.007;
                                    if nargin<9
                                        xMargin=-0.007;
                                        if nargin<8
                                            doubleClicker=[];
                                            if nargin<7
                                                priorFcn=[];
                                                if nargin<6
                                                    reset=true;
                                                    if nargin<5
                                                        doContours=10;
                                                        if nargin<4
                                                            lblMap=[];
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
            [addTrainingHtml, sup1, sup2, trStart, trEnd]=...
                LabelBasics.AddTrainingHtml(lblMap, doJavaLegend);
            [R, C]=size(data);
            assert(C==2, 'data''s 2nd dimension must be size 2!');
            [RL, CL]=size(lbls);
            if RL==1 && CL==R
            else
                assert(CL==1, 'Labels'' 2nd dimension must be size 1!');
                assert(R==RL, 'Labels'' 1st dimension size !- data''s 1st!');
            end
            if stretchPerc>0
                pb=ProbabilityDensity2(...
                    MatBasics.StretchXyData(data, stretchPerc));
            else
                pb=ProbabilityDensity2(data);
            end
            if reset
                cla(ax, 'reset');
            end
            try
                if doContours
                    ms=1;
                else
                    ms=2;
                end
                wasHeld=ishold(ax);
                if ~wasHeld
                    hold(ax, 'on');
                end
                u=unique(lbls);
                N=length(u);
                labelHs=zeros(1,N);
                if isempty(lblMap)
                    for i=1:N
                        l=lbls==u(i);
                        clr=Gui.HslColor(i, N);
                        labelHs(i)=plot(ax, ...
                            data(l,1), data(l,2), '.', 'MarkerSize', ms,...
                            'MarkerEdgeColor', clr, 'LineStyle', 'none');
                    end
                else
                    names={};
                    labelIdxs=[];
                    for i=1:N
                        [key, keyColor, keyTraining]=LabelBasics.Keys(u(i));
                        l=lbls==u(i);
                        colorString=char(lblMap.get(keyColor));
                        if ~isempty(colorString)
                            clr=str2num(colorString);
                            if any(clr>1)
                                clr=clr/256;
                            end
                        else
                            if u(i)==0
                                clr=[.9 .9 .9261];
                            else
                                clr=Gui.HslColor(i, N);
                            end
                            colorString=num2str(floor(clr*256));
                            lblMap.put(keyColor, colorString);
                        end
                        labelHs(i)=plot(ax, ...
                            data(l,1), data(l,2), '.', 'MarkerSize', ms,...
                            'MarkerEdgeColor', clr, 'LineStyle', 'none');
                        name=char(lblMap.get(java.lang.String(key)));
                        if addTrainingHtml
                            name=[name trStart lblMap.get(keyTraining) trEnd];
                        end
                        if ~isempty(name)
                            if doJavaLegend
                                if contains(name, '^{')
                                    name=strrep(name, '^{', sup1);
                                    name=strrep(name, '}', sup2);
                                end
                            end
                            names{end+1}=name;
                            labelIdxs(end+1)=i;
                        else
                            if u(i)==0
                                names{end+1}='unsupervised';
                            else
                                if u(i)<0
                                    names{end+1}=['Subset #' ...
                                        num2str(0-u(i)) ];
                                else
                                    names{end+1}=['Subset #' char(key)];
                                end
                            end
                            labelIdxs(end+1)=i;
                            %disp([id ' NOT found' ]);
                        end
                    end
                    if ~isempty(names)
                        [javaLegend, fncMotion, btns, sortI, plots]=...
                            Plots.Legend(labelHs,...
                             names, labelIdxs, xMargin, yMargin, true,...
                            [], priorFcn, doubleClicker, doJavaLegend, ...
                            oldJavaBtns, southComponent, selectedCallback);
                        btnLbls=u(sortI);
                    end
                end
                if (islogical(doContours) && doContours)
                    hContour=pb.drawContours(ax, 10, [0 0 0]);
                    pb.drawContourOutliers(ax, hContour);
                elseif ~isempty(doContours) && doContours>0
                    hContour=pb.drawContours(ax, doContours, [0 0 0]);
                    pb.drawContourOutliers(ax, hContour);
                end
                if ~wasHeld
                    hold(ax, 'off');
                end
            catch ex
                ex.getReport
            end
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
clearvars;
close all;
clc;
rng default;
e1 = [1 0 0]';
e2 = [0 1 0]';
e3 = [0 0 1]';

basis1 = [e1 e2];
basis2 = [e2 e3];

bases = {basis1, basis2};
cluster_sizes = [5 5];
points = spx.data.synthetic.subspaces.uniform_points_on_subspaces(bases, cluster_sizes);
X = points.X

S = size(X, 2);

C = zeros(S, S);
start_time = tic;
fprintf('Processing %d signals\n', S);
for s=1:S
    fprintf('.');
    if (mod(s, 50) == 0)
        fprintf('\n');
    end
    x = X(:, s);
    cvx_begin
    % storage for  l1 solver
    variable z(S, 1);
    minimize norm(z, 1)
    subject to
    x == X*z;
    z(s) == 0;
    cvx_end
    C(:, s)  = z;
end
elapsed_time  = toc(start_time);
fprintf('\n');
spr_stats = spx.cluster.subspace.subspace_preservation_stats(C, cluster_sizes);

labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes)

spr_flags = zeros(1, S);
spr_errors = zeros(1, S);

C = abs(C);

c1 = C(:, 1);

k = labels(1);

non_zero_indices = (c1 >= 1e-3);

non_zero_labels = labels(non_zero_indices);

spr_flags(1) = all(non_zero_labels == k)

w = labels == k;

c1k = c1(w);

spr_errors(1) = 1 - sum(c1k) / sum (c1);

author: indigits

repo: sparse-plex


text:
function v = tfocs_dot( x, y )

% TFOCS_DOT    Dot products.

if isempty( x ) || isempty( y ),
	v = 0;
else
	v = sum( cellfun( @tfocs_dot, x.value_, y.value_ ) );
end

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

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [ p, seed ] = cylinder_sample_3d ( p1, p2, r, n, seed )

%*****************************************************************************80
%
%% CYLINDER_SAMPLE_3D samples a cylinder in 3D.
%
%  Discussion:
%
%    We are sampling the interior of a right finite cylinder in 3D.
%
%    The interior of a (right) (finite) cylinder in 3D is defined by an axis,
%    which is the line segment from point P1 to P2, and a radius R.  The points
%    on or inside the cylinder are:
%    * points whose distance from the line through P1 and P2 is less than
%      or equal to R, and whose nearest point on the line through P1 and P2
%      lies (nonstrictly) between P1 and P2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3,1), P2(3,1), the first and last points
%    on the axis line of the cylinder.
%
%    Input, real R, the radius of the cylinder.
%
%    Input, integer N, the number of sample points to compute.
%
%    Input, integer SEED, the random number seed.
%
%    Input, real P(3,N), the sample points.
%
%    Output, integer SEED, the random number seed.
%

%
%  Compute the axis vector.
%
  axis(1:3,1) = p2(1:3,1) - p1(1:3,1);
  axis_length = r8vec_norm ( 3, axis );
  axis(1:3,1) = axis(1:3,1) / axis_length;
%
%  Compute vectors V2 and V3 that form an orthogonal triple with AXIS.
%
  [ v2, v3 ] = plane_normal_basis_3d ( p1, axis );
%
%  Assemble the randomized information.
%
  [ radius(1,1:n), seed ] = r8vec_uniform_01 ( n, seed );
  radius(1,1:n) = r * sqrt ( radius(1,1:n) );

  [ theta(1,1:n), seed ] = r8vec_uniform_01 ( n, seed );
  theta(1,1:n) = 2.0 * pi * theta(1,1:n);

  [ z(1,1:n), seed ] = r8vec_uniform_01 ( n, seed );
  z(1,1:n) = axis_length * z(1,1:n);

  for i = 1 : 3

    p(i,1:n) = p1(i,1) ...
              + axis(i,1) * z(1,1:n) ...
              + v2(i,1) * radius(1,1:n) .* cos ( theta(1,1:n) ) ...
              + v3(i,1) * radius(1,1:n) .* sin ( theta(1,1:n) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = isPointInEllipse(point, ellipse, varargin)
%ISPOINTINELLIPSE Check if a point is located inside a given ellipse
%
%   B = isPointInEllipse(POINT, ELLIPSE)
%   Returns true if point is located inside the given ellipse.
%
%   B = isPointInEllipse(POINT, ELLIPSE, TOL)
%   Specifies the tolerance value
%
%   Example:
%   isPointInEllipse([1 0], [0 0 2 1 0])
%   ans =
%       1
%   isPointInEllipse([0 0], [0 0 2 1 0])
%   ans =
%       1
%   isPointInEllipse([1 1], [0 0 2 1 0])
%   ans =
%       0
%   isPointInEllipse([1 1], [0 0 2 1 30])
%   ans =
%       1
%
%   See also:
%   ellipses2d, isPointInCircle
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 11/03/2011
%

%   HISTORY

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% compute ellipse to unit circle transform
rot = createRotation(-deg2rad(ellipse(5)));
sca = createScaling(1./ellipse(3:4));
trans = sca * rot;

% transform points to unit circle basis
pTrans = bsxfun(@minus, point, ellipse(:,1:2));
pTrans = transformPoint(pTrans, trans);

% test if distance to origin smaller than 1
b = sqrt(sum(power(pTrans, 2), 2)) - 1 <= tol;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [index, sequence] = sequence_find(sequences, name)
% sequence_find Find a sequence by its name
%
% Find a sequence by its name in a cell array of sequence structures. Returns index of sequence and its structure.
% If a sequence is not found the function returns an empty matrix.
%
% Input:
% - sequences: Cell array of sequence structures.
% - name: A string containing sequence namer.
%
% Output:
% - index: Index of the sequence in the cell array or empty matrix if not found.
% - sequence: The sequence structure.

index = find(cellfun(@(t) strcmp(t.name, name), sequences, 'UniformOutput', true), 1);

if isempty(index)
    sequence = [];
else
    sequence = sequences{index};
end;

author: votchallenge

repo: toolkit-legacy


text:
function [Fitness, extreme] = CalFitness(PopObj)
% Calculate the fitness of each solution
% F(x) = f(x) + d(x) Algorithm 3

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Li Li

    [N,M]  = size(PopObj); %N  No. of solutions, M No. of Objectives

   %% q
    %Zmax = max(PopObj,[],1); % Identify the max point
    Zmin = min(PopObj,[],1); % Identify the ideal point
    % Identify the extreme points
    W = zeros(M) + 1e-6;
    W(logical(eye(M))) = 1;
    ASF = zeros(N,M);
    for i = 1 : M
        ASF(:,i) = max((PopObj-repmat(Zmin,N,1))./repmat(W(i,:),N,1),[],2);
    end
    [~,extreme] = min(ASF,[],1); %'extreme' is the extreme solutions
    % Calculate the intercepts
    Hyperplane = PopObj(extreme,:)\ones(M,1); % linear equation X=A\B???A*X=B???hyperplane is the solution of linear equation PopObj(extreme,:)x Hyperplane = ones(M,1)
    a = (1./Hyperplane)';
    if any(isnan(a))
        a = max(PopObj,[],1);
    end
    % Normalization
    PopObj = (PopObj-repmat(Zmin,N,1))./repmat(a-Zmin,N,1);

    %---------------w_nad
    w_nad = a; % set a as the nadir point

    Zmin = min(PopObj,[],1); % Identify the ideal point
    for i = 1 : M
        ASF_q(:,i) = max((PopObj-repmat(Zmin,N,1))./repmat(w_nad,N,1),[],2);
    end
    [~,key_point] = min(ASF_q,[],1);
    %q
    q = sqrt(sum(PopObj(key_point).^2,2)) * sqrt(M);
    if q <= 1
        fx = sum(PopObj-repmat(Zmin,N,1),2);
    else
        fx = max(PopObj-repmat(Zmin,N,1),[],2);
    end

    %% Calculate the DMD between each two solutions
    Distance = inf(N);
    for i = 1 : N
        for j = [1:i-1,i+1:N]
            denominator = 1 + min(PopObj(i,:),PopObj(j,:));
            numerator = max(PopObj(i,:),PopObj(j,:)) - min(PopObj(i,:),PopObj(j,:));
            dis = numerator./ denominator;
            Distance(i,j) = sum(dis);
        end
    end

    %% Calculate D(i)
    Distance = sort(Distance,2);
    k = floor(sqrt(N)); %k-NN
    d = 1./(Distance(:,k)+2); % dimensionality margin distance

    %% Calculate the fitnesses
    Fitness = fx+d;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%  Computing an numerical Jordan decomposition
%                      X*J*inv(X)
%  of a given matrix A within a distance threshold theta, formulated under
%  a 'three-strikes principle'.  For details, see
%     Z. Zeng and T.Y. Li,  A numerical algorithm for computing the Jordan
%        Canonical Form, preprint, 2007
%
%  In a nutshell, let A be a matrix whose entries are given approximately
%  with error magnitude being small.  The exact JCF of A will be degraded.
%  However, it is possible to recover the underlying Jordan structure and
%  approximate X and J by NumericalJoranForm
%
%  Syntax:  There are several choices for either LHS or RHS
%        >>   LHS      =      RHS
%  ------------------- | --------------------------------------
%   [J,X]              |   RegularizedJCF(A)
%   [J,X,e,s,t]        |   RegularizedJCF(A,theta)
%                      |   RegularizedJCF(A,theta,tau,gap)
%
% Input:    A -- matrix whose AJCF is to be computed
%       theta -- (optional) distance threshold
%         tau -- (optional) deflation threshold, simple eigenvalues
%                  whose geometric condition numbers above tau will be deflated
%         gap -- (optional) singular value gap used in rank revealing
%
% Output:   J -- The numerical Jordan Canonical Form
%           X -- The principle vector matrix
%           e -- the list of distinct eigenvalues
%           s -- The Jordan block sizes of the eigenvalues
%           t -- the residuals (backward errors) and the staircase condition
%                  numbers of the eigenvalues

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP