text:
%% Bingham Distribution
%
%% Theory
%
% The Bingham distribution has the density function
%
% $$ f(g;K,U) = _1\!F_1 \left(\frac{1}{2},2,K \right)^{-1} \exp
% \left\{ g^T UKU  g \right\},\qquad g\in S^3, $$
%
% where $U$ is an $4 \times 4$ orthogonal matrix with unit quaternions
% $u_{1,..,4}\in S^3$ in the columns and $K$ is a $4 \times 4$ diagonal matrix
% with the entries $k_1,..,k_4$ describing the shape of the distribution.
% $_1F_1(\cdot,\cdot,\cdot)$ is the hypergeometric function with matrix
% argument normalizing the density.
%
% The shape parameters $k_1 \ge k_2 \ge k_3 \ge k_4$ give
%
% * a _bipolar_   distribution, if $k_1 + k_4 > k_2 + k_3$,
% * a _circular_  distribution, if $k_1 + k_4 = k_2 + k_3$,
% * a _spherical_ distribution, if $k_1 + k_4 < k_2 + k_3$,
% * a _uniform_  distribution, if $k_1 = k_2 = k_3 = k_4$,
%
%%
% The general setup of the Bingham distribution in MTEX is done as follows

cs = crystalSymmetry('1');

kappa = [100 90 80 0];   % shape parameters
U     = eye(4);          % orthogonal matrix

odf = BinghamODF(kappa,U,cs)

%%
%

h = [Miller(0,0,1,cs) Miller(1,0,0,cs) Miller(1,1,1,cs)];
plotPDF(odf,h,'antipodal','silent');

%%
%
plot(odf,'sections',6)

%% The bipolar case and unimodal distribution
% First, we define some unimodal odf

odf_spherical = unimodalODF(orientation.rand(cs),'halfwidth',20*degree)

%%
%

plotPDF(odf_spherical,h,'antipodal','silent')

%%
% Next, we simulate individual orientations from this odf, in a scattered
% axis/angle plot in which the simulated data looks like a sphere

ori_spherical = discreteSample(odf_spherical,1000);
close all
scatter(ori_spherical)

%%
% From this simulated EBSD data, we can estimate the parameters of the
% Bingham distribution,

odf_est = calcBinghamODF(ori_spherical)

plotPDF(odf_est,h,'antipodal','silent')

%% TODO
%
% where |U| is the orthogonal matrix of eigenvectors of the orientation
% tensor and |kappa| the shape parameters associated with the |U|.
%
% next, we test the different cases of the distribution on rejection

%T_spherical = bingham_test(ori_spherical,'spherical','approximated');
%T_oblate    = bingham_test(ori_spherical,'prolate',  'approximated');
%T_prolate   = bingham_test(ori_spherical,'oblate',   'approximated');

%t = [T_spherical T_oblate T_prolate]

%%
% The spherical test case failed to reject for some level of
% significance, hence we would dismiss the hypothesis prolate and oblate.

%df_spherical = BinghamODF(kappa,U,crystalSymmetry,specimenSymmetry)

%%
%

%plotPDF(odf_spherical,h,'antipodal','silent')

%% Prolate case and fibre distribution
% The prolate case correspondes to a fibre.

odf_prolate = fibreODF(Miller(0,0,1,crystalSymmetry('1')),zvector,...
  'halfwidth',20*degree)

%%
%

plotPDF(odf_prolate,h,'upper','silent')

%%
% As before, we generate some random orientations from a model odf. The
% shape in an axis/angle scatter plot reminds of a cigar

ori_prolate = discreteSample(odf_prolate,1000);
close all
scatter(ori_prolate)

%%
% We estimate the parameters of the Bingham distribution

calcBinghamODF(ori_prolate)

%%
% and test on the three cases

%T_spherical = bingham_test(ori_prolate,'spherical','approximated');
%T_oblate    = bingham_test(ori_prolate,'prolate',  'approximated');
%T_prolate   = bingham_test(ori_prolate,'oblate',   'approximated');

%t = [T_spherical T_oblate T_prolate]

%%
% The test clearly rejects the spherical and prolate case, but not the
% prolate. We construct the Bingham distribution from the parameters, it
% might show some skewness

odf_prolate = BinghamODF(kappa,U,crystalSymmetry,specimenSymmetry)

%%
%

plotPDF(odf_prolate,h,'antipodal','silent')

%% Oblate case
% The oblate case of the Bingham distribution has no direct counterpart in
% terms of texture components, thus we can construct it straightforward

odf_oblate = BinghamODF([50 50 50 0],eye(4),crystalSymmetry,specimenSymmetry)

%%
%

plotPDF(odf_oblate,h,'antipodal','silent')

  %%
% The oblate cases in axis/angle space remind on a disk

ori_oblate = discreteSample(odf_oblate,1000);
close all
scatter(ori_oblate)

%%
% We estimate the parameters again

calcBinghamODF(ori_oblate)

%%
% and do the tests

%T_spherical = bingham_test(ori_oblate,'spherical','approximated');
%T_oblate    = bingham_test(ori_oblate,'prolate',  'approximated');
%T_prolate   = bingham_test(ori_oblate,'oblate',   'approximated');

%t = [T_spherical T_oblate T_prolate]

%%
% the spherical and oblate case are clearly rejected, the prolate case
% failed to reject for some level of significance

odf_oblate = BinghamODF(kappa, U,crystalSymmetry,specimenSymmetry)

%%
%

plotPDF(odf_oblate,h,'antipodal','silent')

%%
% *Bingham unimodal ODF*

% a modal orientation
cs = crystalSymmetry('-3m');
mod = orientation.byEuler(45*degree,0*degree,0*degree,cs);

% the corresponding Bingham ODF
odf = BinghamODF(20,mod)

plot(odf,'sections',6,'silent','contourf','sigma')

%%
% *Bingham fibre ODF*

odf = BinghamODF([-10,-10,10,10],quaternion(eye(4)),cs)

plot(odf,'sections',6,'silent','sigma')

%%
% *Bingham spherical ODF*

odf = BinghamODF([-10,10,10,10],quaternion(eye(4)),cs)

plot(odf,'sections',6,'silent','sigma');

author: mtex-toolbox

repo: mtex

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


text:
function V=computeV(geo,angles,alphablocks,orig_index,varargin)
if nargin>=6
    [gpuids] = parse_inputs(varargin{1:length(varargin)});
elseif nargin==4
    gpuids = GpuIds();
else
    error('Wrong amount of inputs, 4 or 6 expected');
end

V=zeros(geo.nVoxel(1),geo.nVoxel(2) ,length(alphablocks),'single');
%geo.offDetector(1) = geo.offDetector(1) + (geo.DSD(1) / geo.DSO(1)) * geo.COR(1);
geo=checkGeo(geo,angles);

if ~isfield(geo,'mode')||~strcmp(geo.mode,'parallel')
    for ii=1:length(alphablocks)
        auxang=alphablocks{ii};
        auxindex=orig_index{ii};
        auxgeo = geo;
        % expand the detector to avoiding zeros in backprojection
        maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2);
        auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO);
        auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector;
        % subset of projection angles
        auxgeo.DSD = geo.DSD(auxindex);
        auxgeo.DSO = geo.DSO(auxindex);
        auxgeo.offOrigin = geo.offOrigin(:,auxindex);
        auxgeo.offDetector = geo.offDetector(:,auxindex);
        auxgeo.rotDetector = geo.rotDetector(:,auxindex);
        auxgeo.COR = geo.COR(auxindex);
        %auxgeo=geo;
        V(:,:,ii) = mean(Atb(ones(geo.nDetector(2),geo.nDetector(1),size(auxang,2),'single'),auxgeo,auxang,'gpuids',gpuids),3)+0.000001;
    end
    V(V==0.0)=Inf;
else
    for ii=1:length(alphablocks)
        V(:,:,ii)=ones(geo.nVoxel(1:2).','single')*length(alphablocks{ii});
    end
end
end

function [gpuids]=parse_inputs(varargin)
    %fprintf('parse_inputs0(varargin (%d))\n', length(varargin));
    if isempty(varargin)
        gpuids = GpuIds();
    else
        % create input parser
        p=inputParser;
        % add optional parameters
        addParameter(p,'gpuids', GpuIds());
        %execute
        parse(p,varargin{:});
        %extract
        gpuids=p.Results.gpuids;
    end
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
classdef TestComposeRT
    %TestComposeRT

    properties (Constant)
        fields = {'rvec3', 'tvec3', 'dr3dr1', 'dr3dt1', ...
            'dr3dr2', 'dr3dt2', 'dt3dr1', 'dt3dt1', 'dt3dr2', 'dt3dt2'};
    end

    methods (Static)
        function test_1
            rvec1 = zeros(3,1);
            tvec1 = zeros(3,1);
            rvec2 = zeros(3,1);
            tvec2 = zeros(3,1);
            S = cv.composeRT(rvec1, tvec1, rvec2, tvec2);
            validateattributes(S, {'struct'}, {'scalar'});
            assert(all(ismember(TestComposeRT.fields, fieldnames(S))));
            cellfun(@(f) validateattributes(S.(f), {'double'}, ...
                {'vector', 'numel',3}), TestComposeRT.fields(1:2));
            cellfun(@(f) validateattributes(S.(f), {'double'}, ...
                {'2d', 'size',[3 3]}), TestComposeRT.fields(3:end));
        end

        function test_2
            rvec1 = rand(3,1, 'single');
            tvec1 = rand(3,1, 'single');
            rvec2 = rand(3,1, 'single');
            tvec2 = rand(3,1, 'single');
            S = cv.composeRT(rvec1, tvec1, rvec2, tvec2);
            validateattributes(S, {'struct'}, {'scalar'});
            assert(all(ismember(TestComposeRT.fields, fieldnames(S))));
            cellfun(@(f) validateattributes(S.(f), {'single'}, ...
                {'vector', 'numel',3}), TestComposeRT.fields(1:2));
            cellfun(@(f) validateattributes(S.(f), {'single'}, ...
                {'2d', 'size',[3 3]}), TestComposeRT.fields(3:end));
        end

        function test_error_argnum
            try
                cv.composeRT();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Add a landmark to the UKF.
% We have to compute the uncertainty of the landmark given the current state
% (and its uncertainty) of the newly observed landmark. To this end, we also
% employ the unscented transform to propagate Q (sensor noise) through the
% current state

function [mu, sigma, map] = add_landmark_to_map(mu, sigma, z, map, Q);

% For computing sigma
global scale;

landmarkId = z.id;

%add landmark to the map
map = [map; landmarkId];
% TODO: Initialize its pose according to the measurement and add it to mu

% Append the measurement to the state vector
mu = [mu; z.range(); z.bearing()];
% Initialize its uncertainty and add it to sigma
sigma = blkdiag(sigma, Q);

% Transform from [range, bearing] to the x/y location of the landmark
% This operation intializes the uncertainty in the position of the landmark
% Sample sigma points
sig_pnts_new = compute_sigma_points(mu, sigma);
% Normalize!
sig_pnts_new(3,:) = normalize_angle(sig_pnts_new(3,:));
% Compute the xy location of the new landmark according to each sigma point
newX = sig_pnts_new(1,:) + sig_pnts_new(end-1,:).*cos(sig_pnts_new(3,:) + sig_pnts_new(end,:));
newY = sig_pnts_new(2,:) + sig_pnts_new(end-1,:).*sin(sig_pnts_new(3,:) + sig_pnts_new(end,:));
% The last 2 components of the sigma points can now be replaced by the xy pose of the landmark
sig_pnts_new(end-1,:) = newX;
sig_pnts_new(end,:) = newY;

% Recover mu and sigma
n = length(mu);
lambda = scale - n;
w0 = lambda/scale;
wm = [w0, repmat(1/(2*scale),1,2*n)];
wc = wm;
% Theta should be recovered by summing up the sines and cosines
cosines = sum(cos(sig_pnts_new(3,:)).*wm);
sines = sum(sin(sig_pnts_new(3,:)).*wm);
% recompute the angle and normalize it
mu_theta = normalize_angle(atan2(sines, cosines));
mu = sum(sig_pnts_new .* repmat(wm, rows(sig_pnts_new), 1), 2);
mu(3) = mu_theta;

diff = sig_pnts_new - repmat(mu,1,size(sig_pnts_new,2));
% Normalize!
diff(3,:) = normalize_angle(diff(3,:));
sigma = (repmat(wc, rows(diff), 1) .* diff) * diff';

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym numel (@var{x})
%% Return number of elements in symbolic array.
%%
%% Example:
%% @example
%% @group
%% syms x
%% A = [1 2 x; x 3 4];
%% numel(A)
%%   @result{} 6
%% @end group
%% @end example
%%
%% @seealso{@@sym/length, @@sym/size}
%% @end defmethod

function n = numel(x)

  %disp('numel call')
  d = size(x);
  n = prod(d);

end

%!test
%! a = sym([1 2 3]);
%! assert(numel(a) == 3);

%!test
%! % 2D array
%! a = sym([1 2 3; 4 5 6]);
%! assert(numel(a) == 6);

%!test
%! % empty
%! a = sym([]);
%! assert(numel(a) == 0);

repo: octsympy


text:
function out = CO_FirstMin(y,minWhat,extraParam,minNotMax)
% CO_FirstMin  Time of first minimum in a given self-correlation function
%
%---INPUTS:
% y, the input time series
% minWhat, the type of correlation to minimize: either 'ac' for autocorrelation,
%           or 'mi' for automutual information. By default, 'mi' specifies the
%           'gaussian' method from the Information Dynamics Toolkit. Other
%           options can also be implemented as 'mi-kernel', 'mi-kraskov1',
%           'mi-kraskov2' (all from Information Dynamics Toolkit implementations),
%           or 'mi-hist' (histogram-based method).
% extraParam, an additional parameter required from the minWhat (e.g., Kraskov)
% minNotMax, return the max instead of min.
%
% Note that selecting 'ac' is unusual operation: standard operations are the
% first zero-crossing of the autocorrelation (as in CO_FirstCrossing), or the first
% minimum of the mutual information function ('mi').

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check inputs:
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(minWhat)
    % Mutual information using gaussian method from Information Dynamics Toolkit:
    minWhat = 'mi-gaussian';
end
if nargin < 3
    extraParam = [];
end
if nargin < 4
    minNotMax = true;
end

N = length(y); % Time-series length

% ------------------------------------------------------------------------------
% Define the autocorrelation function
% ------------------------------------------------------------------------------
switch minWhat
case {'mi','mi-gaussian'} % default method (using Information Dynamics Toolkit)
    corrfn = @(x) IN_AutoMutualInfo(y,x,'gaussian');
case 'mi-kernel' % (using Information Dynamics Toolkit)
    corrfn = @(x) IN_AutoMutualInfo(y,x,'kernel');
case 'mi-kraskov1' % (using Information Dynamics Toolkit)
    corrfn = @(x) IN_AutoMutualInfo(y,x,'kraskov1');
case 'mi-kraskov2' % (using Information Dynamics Toolkit)
    % extraParam is the number of nearest neighbors:
    corrfn = @(x) IN_AutoMutualInfo(y,x,'kraskov2',extraParam);
case 'mi-hist'
    % Automutual information implemented in super-naive box counting as in BF_MutualInformation
    corrfn = @(x) BF_MutualInformation(y(1:end-x), y(1+x:end), 'range', 'range',extraParam);
case {'ac','corr'}
    % Autocorrelation implemented as CO_AutoCorr
    corrfn = @(x) CO_AutoCorr(y,x,'Fourier');
otherwise
    error('Unknown correlation type specified: ''%s''',minWhat);
end

% ------------------------------------------------------------------------------
% Search for a minimum
% ------------------------------------------------------------------------------
% (Incrementally through time lags until a minimum is found)

autoCorr = zeros(N-1,1); % pre-allocate maximum length autocorrelation vector
if minNotMax
    %---------------------------------------------------------------------------
    % FIRST LOCAL MINUMUM
    %---------------------------------------------------------------------------
    for i = 1:N-1
        % Calculate the auto-correlation at this lag:
        autoCorr(i) = corrfn(i);

        % Hit a NaN before got to a minimum -- there is no minimum
        if isnan(autoCorr(i))
            warning('No minimum in %s [[time series too short to find it?]]',minWhat)
            out = NaN;
            return
        end

        % We're at a minimum:
        if i==2 && (autoCorr(2) > autoCorr(1))
            % already increases at lag of 2 from lag of 1: a minimum (since ac(0) is maximal)
            out = 1;
            return
        elseif (i > 2) && (autoCorr(i-2) > autoCorr(i-1)) && (autoCorr(i-1) < autoCorr(i)); % minimum at previous i
            out = i-1; % I found the first minimum!
            return
        end
    end
else
    %---------------------------------------------------------------------------
    % FIRST LOCAL MAXIMUM
    %---------------------------------------------------------------------------
    for i = 1:N-1
        % Calculate the auto-correlation at this lag:
        autoCorr(i) = corrfn(i);

        % Hit a NaN before got to a max -- there is no max
        if isnan(autoCorr(i))
            warning('No maximum in %s [[time series too short to find it?]]',minWhat)
            out = NaN;
            return
        end

        % We're at a local maximum:
        if (i > 2) && (autoCorr(i-2) < autoCorr(i-1)) && (autoCorr(i-1) > autoCorr(i)); % minimum at previous i
            out = i-1; % I found the first maximum!
            return
        end
    end
end

% Still decreasing -- no minimum was found after searching all across the time series:
out = N;

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% -  Unary minus.
%    -A negates the elements of A.
%
%    B = UMINUS(A) is called for the syntax '-A' when A is an object.
%
%    Reference page in Doc Center
%       doc uminus
%
%    Other functions named uminus
%
%       calendarDuration/uminus    fints/uminus       sym/uminus
%       codistributed/uminus       gpuArray/uminus    ts/uminus
%       duration/uminus
%

author: jmaih


text:
%LSPB  Linear segment with parabolic blend
%
% [S,SD,SDD] = LSPB(S0, SF, M) is a scalar trajectory (Mx1) that varies
% smoothly from S0 to SF in M steps using a constant velocity segment and
% parabolic blends (a trapezoidal velocity profile).  Velocity and
% acceleration can be optionally returned as SD (Mx1) and SDD (Mx1)
% respectively.
%
% [S,SD,SDD] = LSPB(S0, SF, M, V) as above but specifies the velocity of
% the linear segment which is normally computed automatically.
%
% [S,SD,SDD] = LSPB(S0, SF, T) as above but specifies the trajectory in
% terms of the length of the time vector T (Mx1).
%
% [S,SD,SDD] = LSPB(S0, SF, T, V) as above but specifies the velocity of
% the linear segment which is normally computed automatically and a time
% vector.
%
% LSPB(S0, SF, M, V) as above but plots S, SD and SDD versus time in a single
% figure.
%
% Notes::
% - If M is given
%   - Velocity is in units of distance per trajectory step, not per second.
%   - Acceleration is in units of distance per trajectory step squared, not
%     per second squared.
% - If T is given then results are scaled to units of time.
% - The time vector T is assumed to be monotonically increasing, and time
%   scaling is based on the first and last element.
% - For some values of V no solution is possible and an error is flagged.
%
% References::
% - Robotics, Vision & Control, Chap 3,
%   P. Corke, Springer 2011.
%
% See also TPOLY, JTRAJ.

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

%TODO
% add a 'dt' option, to convert to everything to units of seconds

function [s,sd,sdd] = lspb(q0, q1, t, V)

    t0 = t;
    if isscalar(t)
        t = (0:t-1)';
    else
        t = t(:);
    end
    plotsargs = {'Markersize', 16};

    tf = max(t(:));

    if nargin < 4
        % if velocity not specified, compute it
        V = (q1-q0)/tf * 1.5;
    else
        V = abs(V) * sign(q1-q0);
        if abs(V) < abs(q1-q0)/tf
            error('V too small');
        elseif abs(V) > 2*abs(q1-q0)/tf
            error('V too big');
        end
    end

    if q0 == q1
        s = ones(size(t)) * q0;
        sd = zeros(size(t));
        sdd = zeros(size(t));
        return
    end

    tb = (q0 - q1 + V*tf)/V;
    a = V/tb;

    p = zeros(length(t), 1);
    pd = p;
    pdd = p;

    for i = 1:length(t)
        tt = t(i);

        if tt <= tb
            % initial blend
            p(i) = q0 + a/2*tt^2;
            pd(i) = a*tt;
            pdd(i) = a;
        elseif tt <= (tf-tb)
            % linear motion
            p(i) = (q1+q0-V*tf)/2 + V*tt;
            pd(i) = V;
            pdd(i) = 0;
        else
            % final blend
            p(i) = q1 - a/2*tf^2 + a*tf*tt - a/2*tt^2;
            pd(i) = a*tf - a*tt;
            pdd(i) = -a;
        end
    end

    switch nargout
        case 0
            if isscalar(t0)
                % for scalar time steps, axis is labeled 1 .. M
                xt = t+1;
            else
                % for vector time steps, axis is labeled by vector M
                xt = t;
            end

            clf
            subplot(311)
            % highlight the accel, coast, decel phases with different
            % colored markers
            hold on
            %plot(xt, p);
            k = t<= tb;
            plot(xt(k), p(k), 'r.-', plotsargs{:});
            k = (t>=tb) & (t<= (tf-tb));
            plot(xt(k), p(k), 'b.-', plotsargs{:});
            k = t>= (tf-tb);
            plot(xt(k), p(k), 'g.-', plotsargs{:});
            grid; ylabel('$s$', 'FontSize', 16, 'Interpreter','latex');

            hold off

            subplot(312)
            plot(xt, pd, '.-', plotsargs{:});
            grid;
            if isscalar(t0)
                ylabel('$ds/dk$', 'FontSize', 16, 'Interpreter','latex');
            else
                ylabel('$ds/dt$', 'FontSize', 16, 'Interpreter','latex');
            end

            subplot(313)
            plot(xt, pdd, '.-', plotsargs{:});
            grid;
            if isscalar(t0)
                ylabel('$ds^2/dk^2$', 'FontSize', 16, 'Interpreter','latex');
            else
                ylabel('$ds^2/dt^2$', 'FontSize', 16, 'Interpreter','latex');
            end

            if ~isscalar(t0)
                xlabel('t (seconds)')
            else
                xlabel('k (step)');
                for c=findobj(gcf, 'Type', 'axes')
                    set(c, 'XLim', [1 t0]);
                end
            end
            shg
        case 1
            s = p;
        case 2
            s = p;
            sd = pd;
        case 3
            s = p;
            sd = pd;
            sdd = pdd;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [Q,c,f,x,info] = quaddecomp(p,z)
%QUADDECOMP Internal function to decompose quadratic expression

[n,m]=size(p);
info = 0;

% Is it a scalar
if (n*m==1)
    % Involved in polynomial expression
    [mt,variabletype] = yalmip('monomtable');
    x_lin = getvariables(p);
    x_var = find(any(mt(x_lin,:),1));
    if nargin==2
        x_var = union(x_var,depends(z));
    end
    x = recover(x_var);
    if all(variabletype(x_lin) ==0)% is(p,'linear')
        n = length(x);
        Q = spalloc(n,n,0);
        fc = getbase(p);
        f = fc(1);
        if nargin==2
            vars = getvariables(p);
            c = zeros(length(x),1);
            for i = 1:length(vars)
                c(find(vars(i)==x_var)) = fc(1+i);
            end
        else
            c = fc(2:end);c=c(:);
        end
        return
    end
    variabletype = variabletype(x_lin);
    if all(variabletype<=2)

        base = getbase(p);
        if nnz(base(1))==0
            f = 0;
            base = base(2:end);
        else
            f = base(1);
            base = base(2:end);
        end
        mt = mt(x_lin,x_var);
        quads   = find (variabletype == 2);
        bilins  = find (variabletype == 1);
        linears  = find (variabletype == 0);
        [varsC,aux1,aux2] = find(mt(linears,:)');
        [varsQ,aux1,aux2] = find(mt(quads,:)');
        [varsB,aux1,aux2] = find(mt(bilins,:)');
        if isempty(varsQ)
            varsQ = [];
        end
        if isempty(varsB)
            varsB = [];
        end
        if isempty(varsC)
            varsC = [];
        end
        c = sparse(varsC,1,base(linears),length(x),1);
        ii = [varsQ ; varsB(1:2:end) ; varsB(2:2:end)];
        jj = [varsQ ; varsB(2:2:end) ; varsB(1:2:end)];
        kk = [base(quads)  base(bilins)/2  base(bilins)/2];
        Q = sparse(ii,jj,kk,length(x),length(x));
    else
        if nargout==5
            info = 1;
            Q = [];
            c = [];
            f = [];
            x = [];
        else
            error('Function is not quadratic');
        end
    end

else
    if nargout==5
        info = 1;
        Q = [];
        c = [];
        f = [];
        x = [];
    else
        error('Function is not scalar');
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout=load_gaimc_graph(graphname)
% LOAD_GAIMC_GRAPH Loads a graph from the gaimc library
%
% load_gaimc_graph is a helper function to load a graph provided with the
% library regardless of the current working directory.
%
% If it's called without any output arguments, it functions just like a
% load command executed on the .mat file with the graph.  If it's called
% with an output arguemnt, it functions just like a load command with
% output arguments.  It's somewhat complicated to explain because this is
% just a convinence function to make the examples work for any path, and
% not just from the gaimc root directory.
%
% Example:
%   % equivalent to load('graphs/airports.mat') run from the gaimc directory
%   load_gaimc_graph('airports')
%   % equivalent to P=load('graphs/kt-7-2.mat') run from the gaimc directory
%   P=load_gaimc_graph('kt-7-2.mat')
%   % so you don't have to put the path in for examples!

% David F. Gleich
% Copyright, Stanford University, 2008-2009

% History
% 2009-04-27: Initial coding

path=fileparts(mfilename('fullpath'));
if nargout==0
    evalin('caller',['load(''' fullfile(path,'graphs',graphname) ''');']);
else
    P = load(fullfile(path,'graphs',graphname));
    varargout{1} = P;
end

author: luanfujun

repo: deep-photo-styletransfer


text:
function [F] = registerjacobian(F)
%REGISTERJACOBIAN Register jacobians to all constraints

F = flatten(F);
Counter = length(F.LMIid);
for i = 1:Counter
    switch F.clauses{i}.type
        case {1,2,3}
            Fi = F.clauses{i}.data;
            F.clauses{i}.data = registerjacobian(Fi);
        otherwise
    end
end


text:
function ima = spm_dilate(varargin)
% Perform a dilation on an image (2D or 3D)
% It uses either the supplied kernel or a standard 6-connectivity kernel.
% FORMAT: ima = spm_dilate(ima)
% or
% FORMAT: ima = spm_dilate(ima,kernel)
%
% Input:
% ima    : 2 or 3D image
% kernel : (Optional) voxel values in ima are replaced by the
%          maximum value in a neighbourhood defined by kernel.
%          The "standard" dilation operation (in 2D) is realised
%          using the kernel
%          0 1 0
%          1 1 1
%          0 1 0
%
% Output:
% ima    : Dilated image.
%
% The functionality of this routine has been modelled on the function
% imdilate from the MATLAB Image processing toolbox. It doesn't (yet)
% have a support function such as strel to help the user to define
% kernels (you have to do it yourself if you want anything above
% 6-connectivty) and it doesnt do the clever structuring element
% decomposition that strel does (and imdilate uses). That should
% in principle mean that spm_dilate is slower than imdilate, but
% at least for small (typical) kernels it is actually more than
% twice as fast.
% The actual job is done by spm_dilate_erode.c that serves both
% spm_dilate.m and spm_erode.m
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: spm_dilate.m 4310 2011-04-18 16:07:35Z guillaume $

if exist('spm_dilate_erode','file')~=3
   error('spm_dilate_erode.c not compiled - see Makefile');
end

if nargin > 1
   kernel = varargin{2};
else
   if length(size(varargin{1})) == 2
      kernel = [0 1 0; 1 1 1; 0 1 0];
   elseif length(size(varargin{1})) == 3
      kernel = cat(3,[0 0 0; 0 1 0; 0 0 0],[0 1 0; 1 1 1; 0 1 0],[0 0 0; 0 1 0; 0 0 0]);
   else
      error('Input ima must be 2- or 3-dimensional');
   end
end

ima = spm_dilate_erode(varargin{1},kernel,'dilate');

author: spm


text:
% gravity_demo.m

% Copyright 2003-2010 The MathWorks, Inc.

clear, close all, clc
get_video_data, disp('press any key to select region'), pause
select_region, disp('press any key to remove background'), pause
remove_background, disp('press any key to segment balls'), pause
segment_by_threshold, disp('press any key to suppress noise'), pause
suppress_noise, disp('press any key to locate ball positions'), pause
locate_ball_positions, disp('press any key to fit circle'), pause
fit_circle, disp('press any key to transform XY points to polar'), pause
%transform_to_polar, disp('press any key to fit damped sinusoid'), pause
%fit_damped_sinusoid, disp('press any key to calibrate pixel size'), pause
transform_to_polar, disp('press any key to model system'), pause
edit model_system
model_system, disp('press any key to calibrate pixel size'), pause
calibrate_resolution, disp('press any key to calculate gravity'), pause
edit calculate_gravity
calculate_gravity

author: Sable

repo: mcbench-benchmarks


text:
function varargout=resgram(f,varargin)
%RESGRAM  Reassigned spectrogram plot
%   Usage: resgram(f,op1,op2, ... );
%          resgram(f,fs,op1,op2, ... );
%
%   `resgram(f)` plots a reassigned spectrogram of *f*.
%
%   `resgram(f,fs)` does the same for a signal with sampling rate *fs*
%   (sampled with *fs* samples per second).
%
%   Because reassigned spectrograms can have an extreme dynamical range,
%   consider always using the `'dynrange'` or `'clim'` options (see below)
%   in conjunction with the `'db'` option (on by default). An example:::
%
%     resgram(greasy,16000,'dynrange',70);
%
%   This will produce a reassigned spectrogram of the |greasy| signal
%   without drowning the interesting features in noise.
%
%   `C=resgram(f, ... )` returns the image to be displayed as a matrix. Use this
%   in conjunction with `imwrite` etc. These coefficients are **only** intended to
%   be used by post-processing image tools. Reassignment should be done
%   using the |gabreassign| function instead.
%
%   `resgram` accepts the following additional arguments:
%
%
%     'dynrange',r  Limit the dynamical range to *r* by using a colormap in
%                   the interval $[chigh-r,chigh]$, where *chigh* is the highest
%                   value in the plot. The default value of `[]` means to not
%                   limit the dynamical range.
%
%     'db'         Apply `20*log10` to the coefficients. This makes it possible to
%                  see very weak phenomena, but it might show too much noise. A
%                  logarithmic scale is more adapted to perception of sound.
%                  This is the default.
%
%     'sharp',alpha
%                  Set the sharpness of the plot. If $alpha=0$ the regular
%                  spectrogram is obtained. $alpha=1$ means full
%                  reassignment. Anything in between will produce a partially
%                  sharpened picture. Default is $alpha=1$.
%
%     'lin'        Show the energy of the coefficients on a linear scale.
%
%     'tfr',v      Set the ratio of frequency resolution to time resolution.
%                  A value $v=1$ is the default. Setting $v>1$ will give better
%                  frequency resolution at the expense of a worse time
%                  resolution. A value of $0<v<1$ will do the opposite.
%
%     'wlen',s     Window length. Specifies the length of the window
%                  measured in samples. See help of |pgauss| on the exact
%                  details of the window length.
%
%     'posfreq'    Display only the positive frequencies. This is the
%                  default for real-valued signals.
%
%     'nf'         Display negative frequencies, with the zero-frequency
%                  centered in the middle. For real signals, this will just
%                  mirror the upper half plane. This is standard for complex
%                  signals.
%
%     'tc'         Time centering. Move the beginning of the signal to the
%                  middle of the plot. This is useful for visualizing the
%                  window functions of the toolbox.
%
%     'image'      Use `imagesc` to display the spectrogram. This is the
%                  default.
%
%     'clim',clim  Use a colormap ranging from $clim(1)$ to $clim(2)$. These
%                  values are passed to `imagesc`. See the help on `imagesc`.
%
%     'thr',r      Keep only the largest fraction *r* of the coefficients, and
%                  set the rest to zero.
%
%     'fmax',y     Display *y* as the highest frequency. Default value of `[]`
%                  means to use the Nyquist frequency.
%
%     'xres',xres  Approximate number of pixels along x-axis / time.
%                  The default value is 800
%
%     'yres',yres  Approximate number of pixels along y-axis / frequency
%                  The default value is 600
%
%     'contour'    Do a contour plot to display the spectrogram.
%
%     'surf'       Do a surf plot to display the spectrogram.
%
%     'mesh'       Do a mesh plot to display the spectrogram.
%
%     'colorbar'   Display the colorbar. This is the default.
%
%     'nocolorbar'  Do not display the colorbar.
%
%   In addition to these parameters, `sgram` accepts any of the flags from
%   |setnorm|. The window used to calculate the spectrogram will be
%   normalized as specified.
%
%   See also: sgram

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: NA
%   REFERENCE: NA

% BUG: Setting the sharpness different from 1 produces a black line in the
% middle of the plot.

if nargin<1
  error('Too few input arguments.');
end;

if sum(size(f)>1)>1
  error('Input must be a vector.');
end;

definput.import={'ltfattranslate','setnorm','tfplot'};

% Define initial value for flags and key/value pairs.
definput.flags.wlen={'nowlen','wlen'};
definput.flags.thr={'nothr','thr'};
definput.flags.tc={'notc','tc'};
definput.flags.plottype={'image','contour','mesh','pcolor'};

% Override the setting from tfplot, because SGRAM does not support the
% 'dbsq' setting (it does not make sense).
definput.flags.log={'db','lin'};

if isreal(f)
  definput.flags.posfreq={'posfreq','nf'};
else
  definput.flags.posfreq={'nf','posfreq'};
end;

definput.keyvals.sharp=1;
definput.keyvals.tfr=1;
definput.keyvals.wlen=0;
definput.keyvals.thr=0;
definput.keyvals.fmax=[];
definput.keyvals.xres=800;
definput.keyvals.yres=600;

[flags,kv,fs]=ltfatarghelper({'fs','dynrange'},definput,varargin);

if (kv.sharp<0 || kv.sharp >1)
  error(['RESGRAM: Sharpness parameter must be between (including) ' ...
	 '0 and 1']);
end;

% Downsample
if ~isempty(kv.fmax)
  if ~isempty(kv.fs)
    resamp=kv.fmax*2/fs;
  else
    resamp=kv.fmax*2/length(f);
  end;

  f=fftresample(f,round(length(f)*resamp));
  kv.fs=2*kv.fmax;
end;

Ls=length(f);

if flags.do_posfreq
   kv.yres=2*kv.yres;
end;

try
  [a,M,L,N,Ndisp]=gabimagepars(Ls,kv.xres,kv.yres);
catch
  err=lasterror;
  if strcmp(err.identifier,'LTFAT:noframe')
    error(sprintf(['The signal is too long. RESGRAM cannot visualize all the details.\n' ...
                   'Try a shorter signal or increase the image resolution by calling:\n\n' ...
                   '  sgram(...,''xres'',xres,''yres'',yres);\n\n' ...
                   'for larger values of xres and yres.\n'...
                   'The current values are:\n  xres=%i\n  yres=%i'],kv.xres,kv.yres));
  else
    rethrow(err);
  end;
end;

% Set an explicit window length, if this was specified.
if flags.do_wlen
  kv.tfr=kv.wlen^2/L;
end;

g={'gauss',kv.tfr,flags.norm};

[tgrad,fgrad,c]=gabphasegrad('dgt',f,g,a,M);
coef=gabreassign(abs(c).^2,kv.sharp*tgrad,kv.sharp*fgrad,a);

% Cut away zero-extension.
coef=coef(:,1:Ndisp);

if flags.do_thr
  % keep only the largest coefficients.
  coef=largestr(coef,kv.thr);
end

if flags.do_posfreq
  coef=coef(1:floor(M/2)+1,:);
  plotdgtreal(coef,a,M,'argimport',flags,kv);
else
  plotdgt(coef,a,'argimport',flags,kv);
end;

if nargout>0
  varargout={coef};
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function f = changeMap(f, newdom)
%CHANGEMAP   Map the domain of a DELTAFUN via a linear change of variable.
%   G = MAP(F,NEWDOM), where the DELTAFUN F has a domain [a, b], returns a
%   DELTAFUN G defined on [c, d], where c = NEWDOM(1), d = NEWDOM(2), such that
%       G(x) = F(a*(d - x)/(d - c) + b*(x - c)/(d - c)) for all x in [c, d].

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

% Store the old mapping:
oldMapping = f.funPart.mapping;

% Map the funPart:
f.funPart = changeMap(f.funPart, newdom);

% Grab the new mapping:
newMapping = f.funPart.mapping;

% Map the deltaLocs:
f.deltaLoc = newMapping.For(oldMapping.Inv(f.deltaLoc));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function v = compute_crossp_vf_vf(v1,v2)

% compute_crossp_vf_vf - compute the crossproduct of 2 3D vector fields.
%
%   v = compute_crossp_vf_vf(v1,v2);
%
%   Copyright (c) 2004 Gabriel Peyre

if size(v1,3)~=3
    error('Works only for 3D vector fields.');
end

v = v1;

v(:,:,1) =  v1(:,:,2).*v2(:,:,3)-v1(:,:,3).*v2(:,:,2);
v(:,:,2) = -v1(:,:,1).*v2(:,:,3)+v1(:,:,3).*v2(:,:,1);
v(:,:,3) =  v1(:,:,1).*v2(:,:,2)-v1(:,:,2).*v2(:,:,1);

author: gpeyre

repo: matlab-toolboxes


text:
classdef SubMesher < handle

    properties (Access = public)
       subMesh
       localMesh
    end

    properties (Access = private)
        mesh
        lastNode
    end

    methods (Access = public)

        function obj = SubMesher(cParams)
            obj.init(cParams);
            obj.computeLocalMesh();
            obj.computeSubMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh        = cParams.mesh;
            obj.lastNode    = cParams.lastNode;
        end

        function computeLocalMesh(obj)
            s.connec = [1 2 5;
                2 3 5;
                3 4 5;
                1 5 4];
            s.coord =  obj.computeCoordLocal();
            obj.localMesh = Mesh(s);
        end

        function xIsoCoord = computeCoordLocal(obj)
            xIsoQuad = obj.computeXisoQuad();
            xIsoCent = [0 0];
            xIsoCoord = [xIsoQuad;xIsoCent];
        end

        function computeSubMesh(obj)
            s.connec = obj.computeSubMeshConnec();
            s.coord  = obj.computeSubMeshCoord();
            obj.subMesh = Mesh(s);
        end

        function x = computeXisoQuad(obj)
            int = Interpolation.create(obj.mesh,'LINEAR');
            x = int.pos_nodes;
        end

        function newConnec = computeSubMeshConnec(obj)
            connecE = obj.computeExtendedConnectivities();
            nnodeN   = size(obj.localMesh.connec,2);
            nsubCells = size(obj.localMesh.connec,1);
            nelemB = obj.mesh.nelem;
            nelemN = nsubCells*nelemB;
            newConnec = zeros(nelemN,nnodeN);
            for isubcell = 1:nsubCells
                element = obj.computeGlobalElements(isubcell,nsubCells,nelemN);
                for inode = 1:nnodeN
                    node = obj.localMesh.connec(isubcell,inode);
                    newConnec(element,inode) = connecE(:,node);
                end
            end
        end

        function connecE = computeExtendedConnectivities(obj)
            oldConnec = obj.mesh.connec;
            nnode     = obj.mesh.nnodeElem;
            nnodeE    = nnode + 1;
            newNodes = obj.computeNewNodes();
            connecE(:,1:nnodeE) = [oldConnec newNodes];
        end

        function nNodes = computeNewNodes(obj)
            nelemOld = obj.mesh.nelem;
            nNodes(:,1) = (obj.lastNode+1):(obj.lastNode + nelemOld);
        end

        function newCoord = computeSubMeshCoord(obj)
            xC = obj.mesh.computeBaricenter();
            newCoord = [obj.mesh.coord; xC' ];
        end

    end

    methods (Access = private, Static)

        function gElement = computeGlobalElements(isub,nsubCells,nelem)
            gElement = isub:nsubCells:nelem ;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function imdata = mcmcComputeImageData(im, imsegs)

[imh, imw] = size(imsegs.segimage);
imdata.yim = 1-repmat([(0:imh-1)/(imh-1)]', 1, imw);
imdata.xim = repmat([(0:imw-1)/(imw-1)], imh, 1);

[gx, gy] = gradient(im);
imdata.gradim = sqrt(sum(gx.^2 + gy.^2, 3));

minEdgeLen = sqrt(imh^2+imw^2)*0.02;
[vpdata.lines, vpdata.spinfo] = ...
    APPgetLargeConnectedEdges(rgb2gray(im), minEdgeLen, imsegs);
[vpdata.v, vpdata.vars, vpdata.p, vpdata.hpos] = ...
    APPestimateVp(vpdata.lines, [imh imw], 0);

imdata.vpdata = vpdata;

% get pixels in each superpixe
stats = regionprops(imsegs.segimage, 'PixelIdxList');
imdata.pixlist = {stats.PixelIdxList};

% get boundary pixels for each superpixel
[tx, ty] = gradient(double(imsegs.segimage));
segedgeim = (tx~=0) | (ty ~=0);
stats = regionprops(double(imsegs.segimage).*segedgeim, 'PixelIdxList');
imdata.bndpixlist = {stats.PixelIdxList};
imdata.bndnpix = zeros(size(imdata.bndpixlist));
for s = 1:imsegs.nseg
    try
        imdata.bndnpix(s) = numel(imdata.bndpixlist{s});
    catch
        imdata.bndnpix(s) = 0;
    end
end

% get a random subset of pixels and boundary pixels to speed computation of
% some features
for s = 1:imsegs.nseg
    npix = imsegs.npixels(s);
    nbndpix = imdata.bndnpix(s);
    imdata.nsmpix(s) = min(npix, 1000);
    imdata.nsmbndpix(s) = floor(nbndpix/4);
    if npix>1000
        rind = randperm(npix);
        imdata.smpixlist{s} = imdata.pixlist{s}(rind(1:1000));
    else
        imdata.smpixlist(s) = imdata.pixlist(s);
    end
    rind = randperm(nbndpix);
    try
        imdata.smbndpixlist{s} = imdata.bndpixlist{s}(rind(1:imdata.nsmbndpix(s)));
    catch
        imdata.smbndpixlist{s} = [];
    end
end

author: Cloud-CV

repo: object-proposals


text:
function idx = gui_select_bracketed(listhandle)
% make sure that the user does not select [...] entries (except if not otherwise possible)

contents = cellstr(get(listhandle,'String'));
idx = get(listhandle,'Value');
for k=idx-1:idx-1+length(contents)
    idx = mod(k,length(contents))+1;
    set(listhandle,'Value',idx)
    sel = contents{idx};
    if ~isempty(sel) && sel(1) ~= '['
        break; end
end

repo: REST


text:
function hh = uitabdemo(option)

% UITABDEMO  This script demostrates the use of UITABPANEL.

% Author: Shiying Zhao (zhao@arch.umsl.edu)
% Version: 1.0
% First created: May 02, 2006
% Last modified: May 20, 2006

if nargin<1, option = 1; end
switch option
  case 1
    style1 = 'popup';
    style2 = 'lefttop';
    style3 = 'rightbottom';
  case 2
    style1 = 'popup';
    style2 = 'righttop';
    style3 = 'leftbottom';
  otherwise
    style1 = 'popup';
    style2 = 'centertop';
    style3 = 'centerbottom';
end

width = 0;
hfig = figure(...
  'Name','uitabpanel demo',...
  'NumberTitle','off', ...
  'Menubar','none',...
  'Toolbar','none');

%--------------------------------------------------------------------------
% Creating the tabs
%--------------------------------------------------------------------------
htab(1) = uitabpanel(...
  'Parent',hfig,...
  'Style',style1,...
  'Units','normalized',...
  'Position',[0,0,0.3,1],...
  'FrameBackgroundColor',[0.4314,0.5882,0.8431],...
  'FrameBorderType','etchedin',...
  'Title',{'This is popup tab','Click me!','See also','About'},...
  'PanelHeights',[8,30,10,10],...
  'HorizontalAlignment','left',...
  'FontWeight','bold',...
  'TitleBackgroundColor',[0.9255 0.9490 0.9765],...
  'TitleForegroundColor',[0.1294,0.3647,0.8510],...
  'PanelBackgroundColor',[0.7725,0.8431,0.9608],...
  'PanelBorderType','line',...
  'CreateFcn',@CreateTab1);

htab(2) = uitabpanel(...
  'Parent',hfig,...
  'TabPosition',style2,...
  'Position',[0.3,0.5,0.7,0.5],...
  'Margins',{[0,-1,1,0],'pixels'},...
  'PanelBorderType','line',...
  'Title',{'Introduction','Syntax','Features','Todo'},...
  'CreateFcn',@CreateTab2);

htab(3) = uitabpanel(...
  'Parent',hfig,...
  'TitlePosition',style3,...
  'Position',[0.3,0,0.7,0.5],...
  'Margins',{[0,1,1,0],'pixels'},...
  'Title',{'The Earth','The Globe','More Styles'},...
  'CreateFcn',@CreateTab3,...
  'SelectedItem',3);

if nargout, hh = htab; end

%--------------------------------------------------------------------------
  function CreateTab1(htab,evdt,hpanel,hstatus)
    uicontrol(...
      'Parent',hpanel(1),...
      'Units','normalized',...
      'Position',[.1,.1,.9,.8],...
      'BackgroundColor',get(hpanel(3),'BackgroundColor'),...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      'Popup panel is another   '
      'type of UITABPANEL       '
      'with the ''Style''         '
      'property set to ''Popup''. '
      'Click on its title bar   '
      'toggle open/close of the '
      'popup panel.             ']);

    uicontrol(...
      'Parent',hpanel(2),...
      'Units','normalized',...
      'Position',[.1,.02,.8,.95],...
      'BackgroundColor',get(hpanel(3),'BackgroundColor'),...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      'Here is a long panel. '
      'Using the vertical    '
      'slider to scroll up   '
      'and down for the rest '
      'of content in this    '
      'popup panel.          '
      '                      '
      'All popup panels have '
      'the same vertical     '
      'layout.               '
      '                      '
      'The length of each    '
      'panel in the group can'
      'be specified using the'
      'the property          '
      '''PanelHeights'' at     '
      'creation time, which  '
      'should be a numeric   '
      'array with the same   '
      'length of ''Title'' in  '
      'the unit of           '
      '''Characters''.         '
      '                      '
      'Interactive change of '
      'the length of a panel '
      'in the group can only '
      'be done when it is    '
      'deactivated.          ']);

    uicontrol(...
      'Parent',hpanel(3),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.82],...
      'BackgroundColor',get(hpanel(3),'BackgroundColor'),...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      'tabpanel, by Dirk Tenne '
      'Tab panel example, by   '
      '  Bill York             '
      'TabPanel Constructor    '
      '  v1.3, by Elmar Tarajan'
      'Tab Panel (Yet another  '
      '  one), by Laine Berhane'
      '  Kahsay                '
      'in File Exchange Site.  ']);

    uicontrol(...
      'Parent',hpanel(4),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.76],...
      'BackgroundColor',get(hpanel(3),'BackgroundColor'),...
      'Style','text',...
      'String',[...
      '     UITABPANEL     '
      '     version 1.0    '
      '                    '
      '     Created by     '
      '    Shiying Zhao    '
      '(zhao@arch.umsl.edu)'
      '                    '
      '    May 20, 2006    ']);

    uicontrol(...
      'Parent',hstatus,...
      'Units','normalized',...
      'Position',[0,0,1,.3],...
      'BackgroundColor',get(hstatus,'BackgroundColor'),...
      'ForegroundColor',[0.6,1,1],...
      'HitTest','off',...
      'Style','text',...
      'String',[...
      ' Copyright(c) 2006  '
      '   Shiying Zhao     '
      'All rights reserved.']);

    set(htab,'ResizeFcn',@TabResize1);
  end

%--------------------------------------------------------------------------
  function CreateTab2(htab,evdt,hpanel,hstatus)
    uicontrol(...
      'Parent',hpanel(1),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.8],...
      'BackgroundColor',get(hpanel(1),'BackgroundColor'),...
      'ForegroundColor','w',...
      'HorizontalAlignment','left',...
      'FontSize',12,...
      'Style','text',...
      'String',[ ...
      'UITABPANEL creates a group of tabbed     '
      'panels with a consistent look and feel to'
      'the builtin MATLAB UI objects.           '
      '                                         '
      'A special type of UITABPANEL is the      '
      'popup tabpanel as shown on the left-hand '
      'side of this figure.                     ']);

    uicontrol(...
      'Parent',hpanel(2),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.8],...
      'BackgroundColor','w',...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      '                                                                '
      '   The syntax of using UITABPANEL is the same as                '
      '   UIPANEL except a few property changes. For a complete        '
      '   description, see the online help of UITABPANEL and the       '
      '   MATLAB script of this demo. Examples:                        '
      '                                                                '
      '     htab = uitabpanel( ''Title'',{''A'',''B'',''C'',''D''} )             '
      '     htab = uitabpanel( ''Title'',{''A'',''B''},''Style'',''righttop'' )  '
      '     htab = uitabpanel( ''Title'',{''A'',''B'',''C''},''Position'',my_pos,'
      '                ''CreateFcn'',@my_fun1,''ResizeFcn'',@my_fun2)      ']);

    uicontrol(...
      'Parent',hpanel(3),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.8],...
      'BackgroundColor','w',...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      '                                                             '
      '   The following is a brief summary of some advanced features'
      '   of the current implementation of UITABPANEL:              '
      '   1. UITABPANEL can be resizable if ''ResizeFcn'' is set.     '
      '   2. A callback routine ''SelectionChangeFcn'' can be defined '
      '      through ''set/getappdata'', which will be executed when  '
      '      the selected tab changes.                              '
      '   3. An empty UICONTAINER in each UITABPANEL is             '
      '      designed mainly for a functional purpose, but is can be'
      '      used by users for other purposes in a number of ways.  ']);

    uicontrol(...
      'Parent',hpanel(4),...
      'Units','normalized',...
      'Position',[.1,.1,.8,.8],...
      'BackgroundColor','w',...
      'HorizontalAlignment','left',...
      'Style','text',...
      'String',[...
      '                                                            '
      '   Although the current implementation of UITABPANEL works  '
      '   well for small number of tabs, many improvements and new '
      '   features are still desirable.                            '
      '                                                            '
      '   An immediate task is to add more styles, such as vertical'
      '   tabs and multi-row tabs. However, the most desirable     '
      '   improvement to me is to rewrite UITABPANEL as a          '
      '   MATLAB class, so that it will behave in the same way as  '
      '   the builtin MATLAB user interface (UI) objects.          ']);

    uicontrol(...
      'Parent',hstatus,...
      'Units','normalized',...
      'Position',[0,0.1,1,.68],...
      'BackgroundColor',get(hstatus,'BackgroundColor'),...
      'ForegroundColor','k',...
      'HitTest','off',...
      'Style','text',...
      'String',['This is ',style2,' tab']);
    setappdata(htab,'SelectionChangeFcn',@TabStatus2);

    set(htab,'ResizeFcn',{@TabResize2,[0.5,0.5]});
  end

%--------------------------------------------------------------------------
  function CreateTab3(htab,evdt,hpanel,hstatus)
    uicontrol(...
      'Parent',hpanel(3),...
      'Units','normalized',...
      'Position',[.3,.1,.6,.6],...
      'BackgroundColor',get(hpanel(1),'BackgroundColor'),...
      'ForegroundColor','w',...
      'HorizontalAlignment','left',...
      'FontSize',12,...
      'Style','text',...
      'String',[...
      'Run this demo with the optional '
      'argument 1, 2, and 3 to view all'
      'implemented styles.             ']);

    load('topo.mat','topo','topomap1');
    axes('Parent',hpanel(1));
    PlotTheEarth(topo,topomap1);

    axes('Parent',hpanel(2),'Position',[0,.1,.8,.8]);
    axis square off
    uicontrol(...
      'Parent',hpanel(2),...
      'Units','normalized',...
      'Position',[0.7,.4,.2,.1],...
      'Style','pushbutton',...
      'Callback',{@PlotTheGlobel,topo},...
      'String','Create it!');
    uicontrol(...
      'Parent',hpanel(2),...
      'Units','normalized',...
      'Position',[0.7,.2,.2,.1],...
      'Style','pushbutton',...
      'Callback','cla',...
      'String','Remove it!');

    uicontrol(...
      'Parent',hstatus,...
      'Units','normalized',...
      'Position',[0,0.1,1,.68],...
      'BackgroundColor',get(hstatus,'BackgroundColor'),...
      'ForegroundColor','k',...
      'HitTest','off',...
      'Style','text',...
      'String',['This is ',style3,' tab']);

    set(htab,'ResizeFcn',{@TabResize2,[0,0.5]});
  end

%--------------------------------------------------------------------------
  function TabResize1(hobj,evdt)
    figpos = get(hfig,'Position');
    tabpos = get(hobj,'Position');
    tabpos(4) = figpos(4);
    set(hobj,'Position',tabpos);

    width = tabpos(3)/figpos(3);
  end

%--------------------------------------------------------------------------
  function TabResize2(hobj,evdt,ysiz)
    figpos = get(hfig,'Position');
    tabpos = get(hobj,'Position');
    tabpos([1,3]) = [width,1-width]*figpos(3)+[1,0];
    tabpos([2,4]) = ysiz*figpos(4)+[1,0];
    set(hobj,'Position',tabpos);
  end

%--------------------------------------------------------------------------
  function TabStatus2(hobj,evdt)
    set(get(evdt.Status,'Children'),'String',['page ',num2str(evdt.NewValue)]);
  end

%--------------------------------------------------------------------------
  function PlotTheEarth(topo,topomap1)
    contour(0:359,-89:90,topo,[0 0],'b')
    axis equal
    box on
    set(gca,'XLim',[0 360],'YLim',[-90 90], ...
      'XTick',[0 60 120 180 240 300 360], ...
      'YTick',[-90 -60 -30 0 30 60 90]);
    hold on
    image([0 360],[-90 90],topo,'CDataMapping', 'scaled');
    colormap(topomap1);
  end

%--------------------------------------------------------------------------
  function PlotTheGlobel(hobj,evdt,topo)
    [x,y,z] = sphere(50);
    props.AmbientStrength = 0.1;
    props.DiffuseStrength = 1;
    props.SpecularColorReflectance = .5;
    props.SpecularExponent = 20;
    props.SpecularStrength = 1;
    props.FaceColor= 'texture';
    props.EdgeColor = 'none';
    props.FaceLighting = 'phong';
    props.Cdata = topo;
    surface(x,y,z,props);
    light('position',[-1 0 1]);
    light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]);
  end

%--------------------------------------------------------------------------

end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function obj = addblobs(obj, cl, varargin)
% This is a method for fmridisplay objects that adds blobs to one or more montages and other surface plot(s).
%
% :Usage:
% ::
%
%     obj = addblobs(obj, cl, varargin)
%
% See addthreshblobs and multi_threshold methods for a multiple thresholds version
% render_blobs does most of the hard work.
%
% :Inputs:
%
%   **obj:**
%        an fmridisplay object
%
%   **cl:**
%        a region object. If you're using an fmri_data object pass in region(fmri_data_obj)
%
% :Optional inputs:
%
% There are many optional inputs that control display features of blobs.
% These are determined by render_blobs
%
%   **COLOR:**
%
%   **'color':**
%        followed by color vector, e.g., [0 1 1], for solid-color blobs
%
%   **'maxcolor':**
%        followed by color vector for max color range, e.g., [0 1 1]
%
%   **'mincolor':**
%        followed by color vector for min color range, e.g., [0 0 1]
%
%   **'onecolor':**
%        force solid-color blobs
%
%   **'splitcolor':**
%        Positive and negative values are mapped to different
%        colormaps. Default is +=hot, -=cool colors. Followed
%        optionally by cell array with
%        vectors of 4 colors defining max/min for +/- range,
%        {minnegcolor maxnegcolor minposcolor maxposcolor}
%        e.g., {[0 0 1] [.3 0 .8] [.8 .3 0] [1 1 0]}
%
%   **'indexmap':**
%       Interpret data as indices into colormap. Followed by an n x 3
%       vector of colormap values. Each element of your input data should
%       specify a row of this colormap.
%
%   **'OUTLINING:**
%
%   **''outline'
%
%   **'linewidth':**
%        followed by width value, e.g., 1
%
%        Note: add 'no_surface' to stop updating the existing surface blobs
%
%
%   **'COLOR RANGE:**
%
%   **'cmaprange':**
%        followed by range of values, e.g., [0 40], [-3 3]. Used in
%        color and transparency setting under some circumstances.
%
%   **'TRANSPARENCY:**
%        {'trans', 'transparent','scaledtransparency', 'constanttrans', [val], 'transvalue', [val]}
%
%   **'trans':**
%        Transparent blobs; with no other input, transparency = 0.75 (1 is opaque, 0 is transparent/invisible)
%
%   **'scaledtransparency':**
%        Transparency is a function of voxel value, lower values are more transparent
%
%   **'transvalue':**
%        Followed by width value, e.g., 1. also 'constanttrans'
%
%
% :Other Options:
%
%   **'smooth':**
%        Smooth blobs
%
%   **'contour':**
%
%   **'interp':**
%      Options to pass to interp2 when rendering blobs on standard
%      brain. See help interp2 for valid interpolation methods
%
%   **'no_surface':**
%        Do not add blobs to surface handles, if they exist
%
%   **'noverbose'** : turn off verbose text reporting
%
%   **'depth'** : followed by number (mm) changes search depth (default=2)
%
% CONTROL OF WHICH MONTAGE
%
%   **'wh_montages':**
%        followed by vector of montage numbers as they appear in
%        the list of registered montages in the fmridisplay object
%
% CONTROL OF WHICH SURFACE
%
%   **'wh_surfaces':**
%        followed by vector of surface numbers as they appear in
%                   the list of registered surfaces in the fmridisplay object
%
% :Examples:
% ::
%
%    obj = addblobs(obj, cl, 'color', [0 1 1]);
%    obj = addblobs(obj, cl, 'color', [0 0 1], 'outline');
%    obj = addblobs(obj, cl, 'color', [0 1 0], 'outline', 'linewidth', 1, 'smooth');
%    obj = addblobs(obj, cl, 'color', [1 0 0], 'smooth', 'cmaprange', [0 40]);
%    obj = addblobs(obj, cl, ... 'wh_montages', 1);
%    obj = addblobs(obj, cl, 'splitcolor', ... 'cmaprange', ... 'trans');
%
% Add only to montage 2 in vector of montages in obj.montage
% ::
%
%    obj = addblobs(obj, cl, 'which_montages', 2);
%
% Map values to the colormap red->yellow.
%
% This uses the default percentile-based mapping so that 20% of voxels
% will have the low color and 20% will have the high color, and the rest will be in between:
% ::
%
%    obj = addblobs(obj, cl, 'maxcolor', [1 1 0], 'mincolor', [1 0 0]);
%
% Same, but now Map a specific range of values in image ([0 to .05])
% ::
%
%    obj = addblobs(obj, cl, 'maxcolor', [1 1 0], 'mincolor', [1 0 0], 'cmaprange', [0 .05]);
%
% Separate positive and negative activations and map to a split colormap;
%
% See render_blobs
% ::
%
%    o2 = addblobs(o2, cl, 'splitcolor', {[0 0 1] [.3 0 .8] [.8 .3 0] [1 1 0]}, 'wh_montages', 1);
%
% It is possible to transparency-map values in a statistic image so you
% can show 'unthresholded' statistic values.  e.g.:
% ::
%
%    o2 = addblobs(o2, cl, 'splitcolor', {[0 0 1] [0 1 1] [1 .5 0] [1 1 0]}, 'cmaprange', [-2 2], 'trans', 'scaledtransparency');
%
% ..
%    Copyright Tor Wager, 2011
% ..

% ..
%    Add the volume to activation maps in fmridisplay object
% ..

% Check and convert to region
% -------------------------------------------------------------------
if isstruct(cl) && ~isa(cl, 'region'), cl = cluster2region(cl); end

if isa(cl, 'image_vector'), cl = region(cl); end

if ~isa(cl, 'region')
    error('cl input must be a region object. Try region() constructor method.');
end

if isempty(cl), return, end

% Map to volume space and interpolate
% -------------------------------------------------------------------

[~, mask] = clusters2mask2011(cl); % turn clusters into mask and volume info

if sum(mask(:)) == 0, warning('No voxels in cl! Empty/zero cl.Z field? Bad cl? No results?'); return, end

XYZ = cat(2, cl(:).XYZ);

dim = max(XYZ, [], 2);

V = struct('mat', cl(1).M, 'dim', dim');

SPACE = map_to_world_space(V);

obj.activation_maps{end + 1} = struct('mapdata', mask, 'V', V, 'SPACE', SPACE, 'blobhandles', [], 'cmaprange', [], ...
    'mincolor', [0 0 1], 'maxcolor', [1 0 0], 'color', []);

% Montage selection and options
% -------------------------------------------------------------------

% select which montages; default = all
wh_montage = 1:length(obj.montage);

whm = strcmp(varargin, 'wh_montages') | strcmp(varargin, 'wh_montage') | strcmp(varargin, 'which_montages') | strcmp(varargin, 'which montages');
if any(whm)
    whm = find(whm);
    wh_montage = varargin{whm(1) + 1};
end

% select which surfaces; default = all
wh_surface = 1:length(obj.surface);
addsurfaceblobs = 1;

whs = strcmp(varargin, 'wh_surfaces') | strcmp(varargin, 'wh_surface') | strcmp(varargin, 'which_surfaces') | strcmp(varargin, 'which surfaces');
if any(whs)
    whs = find(whs);
    wh_surface = varargin{whs(1) + 1};
end

% Default color values
% -------------------------------------------------------------------

% default values
dosplitcolor = 1;
doonecolor = 0;
domaxcolor = 0;
domincolor = 0;

maxposcolor = [1 1 0]; % max pos, most extreme values
minposcolor = [1 .4 .5]; % [.8 .3 0]; % min pos
maxnegcolor = [0 .8 .8]; % [.3 .6 .9]; % max neg
minnegcolor = [0 0 1]; % min neg, most extreme values
indexmap = [];

% More purply orange and deeper blue:  {[0 0 1] [.3 0 .8] [.8 .3 0] [1 1 0]}
% more straight orange to yellow, cyan: {[0 0 1] [0 1 1] [1 .5 0] [1 1 0]}

%add_splitcolor_to_varargin = 0; % internal control, not input option
% pos_colormap = colormap_tor([1 0 .5], [1 1 0], [.9 .6 .1]);  %reddish-purple to orange to yellow
% neg_colormap = colormap_tor([0 0 1], [0 1 1], [.5 0 1]);  % cyan to purple to dark blue

depth = 2;

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

            case {'color', 'onecolor', 'solid'} % many are passed into render_blobs - don't need to do anything
                dosplitcolor = 0;
                 doonecolor = 1;

            case 'maxcolor'
                dosplitcolor = 0;
                 domaxcolor = 1;
                 maxcolor = varargin{i + 1};

            case 'mincolor'
                dosplitcolor = 0;
                 domincolor = 1;
                 mincolor = varargin{i + 1};

            case 'indexmap'
                for exclusive = {'color', 'onecolor', 'solid', 'maxcolor','mincolor','splitcolor'}
                    assert(~contains(exclusive{1},varargin(cellfun(@ischar,varargin))),...
                        sprintf('Cannot evaluate addblobs() with both ''indexmap'' and ''%s'' arguments. These are mutually exclusive',...
                            exclusive{1}));
                end
                dosplitcolor = 0;

                try
                    indexmap = varargin{i+1};
                catch
                    error('''indexmap'' argument must be followed by an n x 3 matrix of colormap values');
                end

            case 'splitcolor'

                if length(varargin) > i && iscell(varargin{i + 1})
                    % we have entered colors

                    splitcolors = varargin{i + 1};
                    if ~iscell(splitcolors) || length(splitcolors) ~= 4
                        error('Enter splitcolor followed by a 4-element cell vector of 4 colors\n{Min_neg} {Max_neg} {Min_pos} {Max_pos}');
                    end
                    maxposcolor = splitcolors{4}; % max pos
                    minposcolor = splitcolors{3}; % min pos
                    maxnegcolor = splitcolors{2}; % max neg
                    minnegcolor = splitcolors{1}; % min neg

                else
                    % use defaults - but add the default colors to varargin
                    % because these are passed on to render_blobs
                    % this will be done below

                end

             case 'depth'
                depth = varargin{i + 1};

            case 'no_surface'
                addsurfaceblobs = 0;

            otherwise %suppress warning because other options passed on.  warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% Add relevant colors and args to varargin, because we may have passed in
% keywords without following args, intending to use defaults specified above
% These are passed to render_blobs and also used to update legend registry
% in fmridisplay obj

if dosplitcolor

    mysplitcolors = {minnegcolor maxnegcolor minposcolor maxposcolor};

    wh = strcmp(varargin, 'splitcolor');
    varargin(wh) = [];  % remove, add to end
    varargin{end + 1} = 'splitcolor';
    varargin{end + 1} = mysplitcolors;

elseif doonecolor
    % do not edit var args, just pass through

else  % add maxcolor/mincolor to varargin

    if domaxcolor

        wh = strcmp(varargin, 'maxcolor');
        varargin(wh) = [];  % remove, add to end
        varargin{end + 1} = 'maxcolor';
        varargin{end + 1} = maxcolor;

    end

    if domincolor

        wh = strcmp(varargin, 'mincolor');
        varargin(wh) = [];  % remove, add to end
        varargin{end + 1} = 'mincolor';
        varargin{end + 1} = mincolor;

    end

end

% Resampling whole map seems to be too slow: do this in render slice...
% fprintf('Resampling map data to underlay space.');
%
% [resampled_dat, SPACEto] = resample_space(mask, V, obj.SPACE);
%
% fprintf(' Done.\n');
%
% obj.activation_maps{end + 1} = struct('mapdata', resampled_dat, 'V_original', V, 'blobhandles', []);

wh_to_display = length(obj.activation_maps);

% Find valid handles and render blobs on them
% -------------------------------------------------------------------------

currentmap = obj.activation_maps{wh_to_display};

if isempty(currentmap) || ~isfield(currentmap, 'mapdata') || isempty(currentmap.mapdata)

    error('Map to display is not a valid activation map on fmridisplay object.');

end

% Montages
% -------------------------------------------------------------------------

for i = wh_montage

    if length(obj.montage) < i
        error('Requested montage does not exist! Check input montage indices.');
    end

    % render, and return color ranges (splitcolor will change these)
    % Note: mincolor, maxcolor used in setting fields in fmridisplay for legend

    [blobhan, cmaprange, mincolor, maxcolor] = render_blobs(currentmap, obj.montage{i}, obj.SPACE, varargin{:});

    % Register blob handles and colormap-defining range in object
    if ~isempty(blobhan)
        obj.activation_maps{wh_to_display}.blobhandles = [obj.activation_maps{wh_to_display}.blobhandles; blobhan];

        obj.activation_maps{wh_to_display}.cmaprange = cmaprange;

        % update color info, for legend (or use defaults..set in
        % render_blobs)_
        for Arg = {'color'}

            whmax = find(strcmp(varargin, Arg{1}));
            if ~isempty(whmax)
                obj.activation_maps{wh_to_display}.(Arg{1}) = [obj.activation_maps{wh_to_display}.(Arg{1}) varargin{whmax(1) + 1}];
            end

        end

        for Arg = {'mincolor' 'maxcolor'}

            str = (['obj.activation_maps{wh_to_display}.(Arg{1}) = ' Arg{1} ';']);
            eval(str);

        end

    end
end

% Surfaces
% -------------------------------------------------------------------------

if addsurfaceblobs
    for i = wh_surface

        if length(obj.surface) < i
            warning('Requested surface does not exist! Check input surface indices.');
            continue
        end

        % Set color maps for + / - values
        if dosplitcolor || doonecolor || domaxcolor || domincolor
            pos_colormap = colormap_tor(minposcolor, maxposcolor);
            neg_colormap = colormap_tor(minnegcolor, maxnegcolor);
        end

        % OLD method: cluster_surf
        % cluster_surf(cl, depth, 'heatmap', 'colormaps', pos_colormap, neg_colormap, obj.surface{i}.object_handle);

        img = region2imagevec(cl);
        if dosplitcolor
            render_on_surface(img, obj.surface{i}.object_handle, 'pos_colormap', pos_colormap, 'neg_colormap', neg_colormap, varargin{:});
        else
            if isempty(indexmap)
                render_on_surface(img, obj.surface{i}.object_handle, varargin{:});
            else
                render_on_surface(img, obj.surface{i}.object_handle, varargin{:},'colormap', indexmap, 'indexmap');
            end
        end

    end
end

end  % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = hlp_getresult(idx,f,varargin)
% Returns the Result-Idx's output of the function, given the supplied arguments.
% Results... = hlp_getresult(Result-Idx, Function, Arguments...)
%
% In:
%   Result-Idx : index of the result (output) of the given function application; can also
%                be a vector of indices, then re-emitting those outputs as outputs.
%                note: if this is passed as a cell array, the outputs will be wrapped into a cell
%                      array.
%
%   Function : function to apply to the Arguments
%
%   Arguments... : list of arguments to the function
%
% Out:
%   Results... : one or more outputs, selected according to Result-Idx from the outputs of
%                Function(Arguments...), and optionally wrapped into a cell array.
%
% Examples:
%   % get the second output of the sort() function, when applied to some data
%   hlp_getresult(2,@sort,[1 4 2 5 2 1 0 6])
%
% See also:
%   hlp_wrapresults
%
%				          		 Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-03-28

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

if isnumeric(idx)
    [tmp{1:max(idx)}] = f(varargin{:});
    varargout = tmp(idx);
elseif iscell(idx)
    [tmp{1:max([idx{:}])}] = f(varargin{:});
    varargout = {tmp([idx{:}])};
else
    error('Unsupported index format (must be either numeric or cell array).');
end


text:
function a = asin(a)
%ASIN         Gradient inverse sine asin(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = 1 ./ sqrt( 1 - sqr(full(a.x(:))) );
  a.x = asin(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function calcAllFeatureSets_HN(pathSet,training,pathMINE,param,setSize,alpha,delta,nBoot,seed,batchNum)
% -------------------------------------------------------------------------
% function calcAllFeatureSets_HN(pathExperiment,trainingStruct,pathMINE,setSize,alpha,delta,nBoot,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set reduction for a given feature set type
% and for all experiments with different degrees of freedom. See ref. [1,2]
% for more details.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathSet: Full path to the Feature Set folder of the corresponding experiment.
%             --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/FSET'
% 2. training: Structure defining all parameters for the given experiment
%              to perform. See masterScript_HN.m for more details.
% 3. pathMINE: Full path to the MINE.jar executable. The executable can be
%              downloaded at: <http://www.exploredata.net/>.
%              --> Ex: '/myProject/radiomics/MultivariableModeling/MINE'
% 4. param: Cell of two strings, defining 1) The feature set type/name; and
%           2) the outcome to model
% 4. setSize: Size of the output feature set (typically set to 25 in ref. [1]).
%             --> Ex: 25
% 5. alpha: Numerical values specifying the coefficient of the first part of
%           the Gain equation, as defined in ref. [1,2].
%           --> Ex: 0.5
% 6. delta: Numerical values specifying the coefficient of the second part
%           of the Gain equation, as defined in ref. [1,2] (third part is set
%           to 0 in this function).
%           --> Ex: 0.5
% 7. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 8. seed: Numerical number to use as seed for bootstrapping experiment
%          --> Ex: 54288
%
% See <https://github.com/mvallieres/radiomics/tree/master/STUDIES/HNstudy/Functions>
% to find computeAllNonTextureFeatures_HN.m, organizeSeparateTextures_HN.m
% and organizeFusedTextures_HN.m. See masterScript_HN.m for a complete
% example of how to utilize the current function.
% -------------------------------------------------------------------------
% OUTPUTS: Feature sets are saved in a folder named 'FSET' in the
% corresponding folder of "pathExperiment"
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off
if nargin < 10
    batchNum = 1;
end

fSetName = param{1};
outcomeName = param{2};
outcome = training.(outcomeName).outcome;
nonText = training.(outcomeName).nonText;
% if strcmp(outcomeName,'Death'), nonText = rmfield(nonText,'PercentInactive'); end % Creates problem for that outcome, this needs to be solved
textCellsName = fieldnames(training.(outcomeName).text.(fSetName)); textCellsName(1:4) = []; % Removing 'param', 'baseline', 'freedom', 'paramName'
nTextCell = numel(textCellsName);
textCells = cell(1,nTextCell);
for i = 1:nTextCell
    textCells{i} = training.(outcomeName).text.(fSetName).(textCellsName{i});
end
paramAll = training.(outcomeName).text.(fSetName).param;
freedomMat = training.(outcomeName).text.(fSetName).freedom;
baseline = training.(outcomeName).text.(fSetName).baseline;

if strcmp(fSetName,'CT')
    nonText = rmfield(nonText,{'SUVmax','SUVpeak','SUVmean','aucCSH','TLG','PercentInactive','gETU'}); % These features were calculated from PET.
end

tic
fprintf(['\n--> COMPUTING "',fSetName,'" FEATURE SET FOR "',outcomeName,'" OUTCOME ... '])
[fSet] = featureSetReduction_HN(pathMINE,fSetName,outcome,setSize,nonText,textCells,textCellsName,paramAll,freedomMat,baseline,alpha,delta,nBoot,seed,batchNum);
cd(pathSet), save(['FSET_',fSetName,'_',outcomeName],'fSet')
fprintf('DONE!\n'), toc

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
classdef PlaneStressTransformer < handle

    properties (Access = protected)
        psTensor
    end

    methods (Access = public,Static)

        function psTensor = transform(Tensor)
            factory       = PlaneStressTransformerFactory();
            planeStresser = factory.create(Tensor);
            psTensor      = planeStresser.getTensor();
        end

    end

    methods (Access = private)

        function t = getTensor(obj)
            t = obj.psTensor;
        end

    end

    methods (Access = protected,Abstract)
        createPlaneStressTensor(obj)
    end

end


text:
function [x, y] = ijcvConfidenceImage2contours(confim, maxy)

[imh, imw] = size(confim);

confim = imfilter(confim, fspecial('gaussian', 25, 7), 'same');

if ~exist('maxy')
    maxy = imh-15;
end

x = [15:imw-15];

cval = [0.25 0.5 0.75];

for c = 1:numel(cval)

    for i = 1:numel(x)
        try
            y{c}(i) = max(find(confim(1:maxy, x(i))>cval(c)));
        catch
            y{c}(i) = 1;
        end
    end

end


text:
% Run this in order to construct the results table
file_out = fopen('results/300VW_66.txt', 'w');

load('results/300VW_OpenFace.mat');
load('../../matlab_version/experiments_300VW/results/cfss_errors.mat');
load('../../matlab_version/experiments_300VW/results/cfan_errors.mat');
load('../../matlab_version/experiments_300VW/results/drmf_errors.mat');
load('../../matlab_version/experiments_300VW/results/iccr_errors.mat');
load('../../matlab_version/experiments_300VW/results/pocr_errors.mat');

fprintf(file_out, 'Errors with outline (66 points)\n');
fprintf(file_out, '------------------------------\n');
fprintf(file_out, 'Method\tcat 1\tcat 2\tcat3\n');
fprintf(file_out, 'CFSS\t%.2f\t%.2f\t%.2f\n', median(cfss_error_66_cat_1)*100, median(cfss_error_66_cat_2)*100, median(cfss_error_66_cat_3)*100);
fprintf(file_out, 'ICCR\t%.2f\t%.2f\t%.2f\n', median(iccr_error_66_cat_1)*100, median(iccr_error_66_cat_2)*100, median(iccr_error_66_cat_3)*100);
fprintf(file_out, 'CFAN\t%.2f\t%.2f\t%.2f\n', median(cfan_error_66_cat_1)*100, median(cfan_error_66_cat_2)*100, median(cfan_error_66_cat_3)*100);
fprintf(file_out, '------------------------------\n');
fprintf(file_out, 'OpenFace CLNF\t%.2f\t%.2f\t%.2f\n', median(clnf_error_66_cat_1)*100, median(clnf_error_66_cat_2)*100, median(clnf_error_66_cat_3)*100);
fprintf(file_out, 'OpenFace CE-CLM\t%.2f\t%.2f\t%.2f\n', median(ceclm_error_66_cat_1)*100, median(ceclm_error_66_cat_2)*100, median(ceclm_error_66_cat_3)*100);
fclose(file_out);

file_out = fopen('results/300VW_49.txt', 'w');

fprintf(file_out, 'Errors without outline (49 points)\n');
fprintf(file_out, '------------------------------\n');
fprintf(file_out, 'Method\tcat 1\tcat 2\tcat3\n');
fprintf(file_out, 'CFSS\t%.2f\t%.2f\t%.2f\n', median(cfss_error_49_cat_1)*100, median(cfss_error_49_cat_2)*100, median(cfss_error_49_cat_3)*100);
fprintf(file_out, 'ICCR\t%.2f\t%.2f\t%.2f\n', median(iccr_error_49_cat_1)*100, median(iccr_error_49_cat_2)*100, median(iccr_error_49_cat_3)*100);
fprintf(file_out, 'DRMF\t%.2f\t%.2f\t%.2f\n', median(drmf_error_49_cat_1)*100, median(drmf_error_49_cat_2)*100, median(drmf_error_49_cat_3)*100);
fprintf(file_out, 'PO-CR\t%.2f\t%.2f\t%.2f\n', median(pocr_error_49_cat_1)*100, median(pocr_error_49_cat_2)*100, median(pocr_error_49_cat_3)*100);
fprintf(file_out, 'CFAN\t%.2f\t%.2f\t%.2f\n', median(cfan_error_49_cat_1)*100, median(cfan_error_49_cat_2)*100, median(cfan_error_49_cat_3)*100);
fprintf(file_out, '------------------------------\n');
fprintf(file_out, 'OpenFace CLNF\t%.2f\t%.2f\t%.2f\n', median(clnf_error_49_cat_1)*100, median(clnf_error_49_cat_2)*100, median(clnf_error_49_cat_3)*100);
fprintf(file_out, 'OpenFace CE-CLM\t%.2f\t%.2f\t%.2f\n', median(ceclm_error_49_cat_1)*100, median(ceclm_error_49_cat_2)*100, median(ceclm_error_49_cat_3)*100);
fclose(file_out);

author: TadasBaltrusaitis

repo: OpenFace


text:
%% Image Alignment using ECC algorithm
%
% This sample demonstrates the use of the ECC image alignment algorithm. When
% one image is given, the template image is artificially formed by a random
% warp given the motion type. Otherwise supply both images. If input warp
% matrix is not specified, the identity transformation is used to initialize
% the algorithm.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/image_alignment.cpp>
%

%% Options
INPUT_IMAGE = fullfile(mexopencv.root(), 'test', 'fruits.jpg');
LOAD_TEMPLATE = false;
LOAD_MAT = false;
MOTION_TYPE = 'Affine';

%%
% type of motion
MOTION_TYPE = validatestring(MOTION_TYPE, ...
    {'Translation', 'Euclidean', 'Affine', 'Homography'});

%%
% termination criteria (ECC iterations and convergnce epsilon)
CRIT = struct('type','Count+EPS', 'maxCount',50, 'epsilon',1e-4);
if CRIT.maxCount > 200
    warning('Too many iterations');
end

%% Input image
% load input image
imgInput = cv.imread(INPUT_IMAGE, 'Grayscale',true);

%% Template image
% initialize or load template image
M0 = [];
if LOAD_TEMPLATE
    % load an existing template image
    fmts = imformats();
    filtspec = strjoin(strcat('*.', [fmts.ext]), ';');
    [fn,fp] = uigetfile(filtspec, 'Select an image');
    if fp==0, error('No file selected'); end
    imgTemplate = cv.imread(fullfile(fp,fn), 'Grayscale',true);
else
    % create a template image by applying a random warp to input image
    imgInput = cv.resize(imgInput, [216 216]);
    opts = {'DSize',[200 200], 'Interpolation','Linear', 'WarpInverse',true};
    switch MOTION_TYPE
        case 'Translation'
            M0 = [1 0 rand*10+10;
                  0 1 rand*10+10];
            imgTemplate = cv.warpAffine(imgInput, M0, opts{:});
        case 'Euclidean'
            theta = pi/30 + pi*randi([-2 2])/180;
            M0 = [cos(theta) -sin(theta) rand*10+10;
                  sin(theta)  cos(theta) rand*10+10];
            imgTemplate = cv.warpAffine(imgInput, M0, opts{:});
        case 'Affine'
            M0 = [1-(rand*0.1-0.05)    rand*0.06-0.03 rand*10+10;
                     rand*0.06-0.03 1-(rand*0.1-0.05) rand*10+10];
            imgTemplate = cv.warpAffine(imgInput, M0, opts{:});
        case 'Homography'
            M0 = [1-(rand*0.1-0.05)    rand*0.06-0.03 rand*10+10;
                     rand*0.06-0.03 1-(rand*0.1-0.05) rand*10+10;
                     rand*2e-4+2e-4    rand*2e-4+2e-4 1];
            imgTemplate = cv.warpPerspective(imgInput, M0, opts{:});
    end
    display(M0)  % ground truth
end
sz = size(imgTemplate);

%% Input warp matrix
% initialize or load warp matrix
M = [];
if LOAD_MAT
    % load from a MAT-file
    uiopen('load');
    assert(~isempty(M), 'Failed to load warp matrix M');
    if strcmp(MOTION_TYPE, 'Homography')
        assert(isequal(size(M), [3 3]));
    else
        assert(isequal(size(M), [2 3]));
    end
else
    % identity matrix
    if strcmp(MOTION_TYPE, 'Homography')
        M = eye(3,3);
    else
        M = eye(2,3);
    end
    warning(['Performance Warning: Identity warp ideally assumes images ' ...
        'of similar size. If the deformation is strong, the identity ' ...
        'warp may not be a good initialization, and estimation may fail.']);
end
display(M)

%% Estimate transformation
fprintf('Estimating "%s" transformation...\n', MOTION_TYPE);
tic
M = cv.findTransformECC(imgTemplate, imgInput, ...
    'InputWarp',M, 'MotionType',MOTION_TYPE, 'Criteria',CRIT);
toc
display(M)

%%
% compare against ground truth
if ~isempty(M0)
    err = norm(M - M0)
end

%% Apply estimated transformation
% warped image
opts = {'DSize',[sz(2) sz(1)], 'Interpolation','Linear', 'WarpInverse',true};
if strcmp(MOTION_TYPE, 'Homography')
    imgWarped = cv.warpPerspective(imgInput, M, opts{:});
else
    imgWarped = cv.warpAffine(imgInput, M, opts{:});
end

%%
% compare against template image
if ~mexopencv.isOctave() && mexopencv.require('vision')
    imgError = imfuse(imgTemplate, imgWarped, 'diff');
else
    imgError = abs(double(imgTemplate) - double(imgWarped));
    imgError = imgError / max(imgError(:));
end

%%
% compute region boundary in input image corresponding to template image
pts = [1 1; sz(2) 1; sz(2) sz(1); 1 sz(1)]; % corners to warp (TL/TR/BR/BL)
pts(:,3) = 1;  % homogeneous coordinates
if strcmp(MOTION_TYPE, 'Homography')
    pts = pts * M.';
else
    pts = pts * [M; 0 0 1].';
end
pts = bsxfun(@rdivide, pts(:,1:2), pts(:,3));

%%
% show results
subplot(221), imshow(imgInput), title('image')
line(pts([1:end 1],1), pts([1:end 1],2), 'Color','m', 'LineWidth',2)
axis tight
subplot(222), imshow(imgTemplate), title('template')
subplot(223), imshow(imgWarped), title('warped image')
subplot(224), imshow(imgError), title('error')

author: kyamagu

repo: mexopencv


text:
function A=vech2Mat(theVec,isSymmetric,offDiagScalFact)
%%VECH2MAT Given the main diagonal and lower-triangular elements of a
%          matrix stacked column-wise, obtain the original
%          symmetric or lower-triangular matrix. This function is
%          essentially the inverse of the vech function when applied to a
%          symmetric or lower-triangular matrix with the proper
%          parameterization.
%
%INPUTS: theVect An nX1 or 1Xn vector holding the diagonal and
%                lower-triangular components of a symmetric matrix stacked
%                column-wise.
%    isSymmetric An optional boolean parameter specifying whether the
%                matrix being recreated is symmetric. If false, a lower-
%                triangular matrix is formed. The default if this parameter
%                is omitted or an empty matrix is passed is true.
%   offDiagScalFact An optional scalar factor by which the off-diagonal
%                elements are scaled. If this parameter is omitted or an
%                empty matrix is passed, then the default value of 1 is
%                used.
%
%OUTPUTS: A The matrix impled by theVec and isSymmetric. This is a dXd
%           matrix, where n=d*(d+1)/2.
%
%The opposite of this function is vech.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(offDiagScalFact))
    offDiagScalFact=1;
end

if(nargin<2||isempty(isSymmetric))
   isSymmetric=true;
end

n=length(theVec);
theVec=theVec(:);

d=(1/2)*(-1+sqrt(1+8*n));
A=zeros(d,d);

curStart=1;
if(isSymmetric)
    for curCol=1:d
        %The diagonal element
        A(curCol,curCol)=theVec(curStart);
        span=(curStart+1):(curStart+d-curCol);

        %The lower-triangular part
        A((curCol+1):end,curCol)=theVec(span)*offDiagScalFact;
        %The upper triangular part.
        A(curCol,(curCol+1):end)=theVec(span)'*offDiagScalFact;
        curStart=curStart+d-curCol+1;
    end
else
    for curCol=1:d
        %The diagonal element
        A(curCol,curCol)=theVec(curStart);
        span=(curStart+1):(curStart+d-curCol);

        %The lower-triangular part
        A((curCol+1):end,curCol)=theVec(span)*offDiagScalFact;

        %The upper-triangular part is just zero.
        curStart=curStart+d-curCol+1;
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Hist] = makeHistogram(Bin,XSize,YSize,NrX,NrY,NrBins)
%  This function classifies the greylevels present in the array image into
%  a greylevel histogram. The pLookupTable specifies the relationship
%  between the greyvalue of the pixel (typically between 0 and 4095) and
%  the corresponding bin in the histogram (usually containing only 128 bins).

Hist=zeros(NrX,NrY,NrBins);

for i=1:NrX
    for j=1:NrY
        bin=Bin(1+(i-1)*XSize:i*XSize,1+(j-1)*YSize:j*YSize);
        for i1=1:XSize
            for j1=1:YSize
                Hist(i,j,bin(i1,j1)) = Hist(i,j,bin(i1,j1)) + 1;
            end
        end
    end
end


text:
% MKLTRNLS  Solve a NLS using MKLTRNLS
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_mkltrnls() INSTEAD!
%
% mktrnls uses the Intel MKL Trust Region NLS Solver.
%
%   [x,fval,exitflag,iter] = mktrnls(fun,grad,x0,ydata,lb,ub,opts)
%
%   Input arguments:
%       fun - nonlinear fitting function handle
%       grad - gradient of nonlinear fitting function handle (optional)
%       x0 - initial solution guess
%       ydata - fitting data
%       lb - decision variable lower bounds (optional)
%       ub - decision variable upper bounds (required if lb present)
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       iter - number of iterations taken by the solver
%       feval - number of function evaluations used by the solver
%
%   Option Fields (all optional):
%       display - solver display level [0,1,2]
%       maxiter - maximum solver iterations
%       maxtime - maximum solver execution time
%       tolafun - absolute function tolerance
%       tolrfun - relative function tolerance
%       iterfun - Iteration Callback Function, stop = iterfun(iter,fval,x)
%
%   Return Status:
%       1 - optimal
%       0 - maximum iterations exceeded
%      -1 - could not converge / tolerance too small
%      -2 - singular or other error
%      -3 - unknown error
%
%
%
%   Copyright (C) 2012 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus