text:
 function [b, d, f, z, mm] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%function [b, d, f, z, mm] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%Function that computes the complex-valued SPSP RF pulse waveform,
%using parameters specified in structure rfp.
%Inputs:
%kp: k-space trajectory parameter structure, generated by kparameterSPSP.m
%rfp:RF waveform design parameter structure, generated by rfparameterSPSP.m
%gz: z gradient waveform vector
%kz,kf: SPSP k-space (kf-kz space) trajectory vectors
%Output:
%b:complex-value SPSP RF pulse waveform in units of Gauss (g).
%
%Chun-yu Yip, 4/1/2009

%Physical parameters
gam = 26751;                       %rad/s/g; gyromagnetic ratio
gambar = gam/2/pi;                 %Hz/g

%Create grid in SPSP space
dresz = rfp.dfovz/rfp.ddimz;       %cm
dresf = rfp.dfovf/rfp.ddimf;       %Hz
z = [-rfp.dfovz/2:dresz:rfp.dfovz/2-dresz];
f = [-rfp.dfovf/2:dresf:rfp.dfovf/2-dresf];
[Z,F] = ndgrid(z,f);

printm('Creating desired pattern for SPSP pulse...');

%phase pattern
TD = rfp.TE-kp.pointtime*kp.npnts/2; %pulse end to acquisition of DC sample
dphs = gam * rfp.alpha * TD * F .* (Z-rfp.zshift); %Eq. 5

%magnitude pattern
if streq(rfp.profileshape,'rect')

    dmag = zeros(size(Z));
    dmag(find((Z<(rfp.slthickz/2)) & (Z>(-rfp.slthickz/2)))) = 1;

elseif streq(rfp.profileshape,'gaussian')

    a = -4*log(2)/rfp.slthickz/rfp.slthickz;
    dmag = exp(a*(Z.^2));

else

    printm('Slice profile shape unknown.');

end

if rfp.smoothprofile

    kernel = gaussian(z,0,rfp.kernelstdv)';
    dmag=convn_fft(dmag,kernel);     %Smoothing implemented in Fourier domain

end;

dmag = dmag/max(max(dmag));
dmag = sin(rfp.flipangle/180 * pi) * dmag;
d = exp(1i*dphs).*dmag;

%Define matrices and vectors for conjugate gradient
b0 = zeros(length(kz),1);            %Initiation of CG
mask = ones(size(d));                %ROI in f-z space
%W = spdiag(mask); % weighting based on ROI definition
W = diag_sp(mask); % jf version

%total energy penalty
beta = rfp.beta;
% C = beta*speye(length(kz)); % original version
C = beta; % jf version

M = length(Z(:));
%precon =speye(length(kz));           %CG Preconditioner
precon = 1; % jf version
Niter = rfp.Niter;                   %Number of CG iterations

%Define matlab object "A" to represent the gigantic system matrix.
%"A" is of type "Fatrix" from Jeff Fessler's toolbox.
printm('Forming spatial-spectral fast object...');
A=spsp_Af(z,f,kz,kf,kp.pointtime,rfp.dfovf,rfp.dfovz);

%%%%%%%%%%%Conjugate Gradient%%%%%%%%%%%%%%%%%%%
printm('Running preconditioned conjugate gradient...');
[b_all, info] = qpwls_pcg(b0, A, W, d(:), zeros(M,1), sqrt(C), precon, Niter);
b = b_all(:,Niter);

if rfp.hamming
	scale = 2;          %This parameter controls how wide the window is.
	hamm=(0.54-0.46.*cos(pi*((kz/scale/max(kz))+1)));
	b= b.*hamm;         %Intended to smooth the excited slice profile.
end

m = A*b;            %Small-tip-angle approximation of excitation pattern
mm = reshape(m,size(d));

if 0 % no, superceded by ir_mri_rf_spsp_plot.m
	printm('Displaying SPSP pulse design results...');
	t = kp.pointtime*[0:1:length(gz)-1].';
	figure;
	subplot(2,1,1);
	plot(t*1000,gz);
	grid;
	x1=xlabel('Time (ms)');
	y1=ylabel('g/cm');
	t1=title('z gradient');
	set(t1,'Fontsize',14);
	set(x1,'Fontsize',12);
	set(y1,'Fontsize',12);

	subplot(2,1,2);
	hold on;
	plot(t*1000,real(b),'b');
	plot(t*1000,imag(b),'r');
	grid;
	x2=xlabel('Time (ms)');
	y2=ylabel('g');
	t2=title('SPSP pulse');
	l2=legend('Real','Imaginary');
	set(t2,'Fontsize',14);
	set(x2,'Fontsize',12);
	set(y2,'Fontsize',12);
	set(l2,'Fontsize',12);

	%Display desired pattern of the SPSP recovery pulse.
	figure;
	subplot(2,1,1);
	imagesc(f,z,abs(d));
	colormap default;colorbar;
	x1 = xlabel('Frequency offset (Hz)');
	y1 = ylabel('z (cm)');
	t1 = title('Desired SPSP pattern (magnitude)');

	set(x1,'Fontsize',12);
	set(y1,'Fontsize',12);
	set(t1,'Fontsize',12);

	subplot(2,1,2);
	imagesc(f,z,angle(d).*abs(d),[-pi,pi]);colorbar;
	x2 =xlabel('Frequency offset (Hz)');
	y2=ylabel('z (cm)');
	t2=title('Desired SPSP pattern (phase)');
	set(x2,'Fontsize',12);
	set(y2,'Fontsize',12);
	set(t2,'Fontsize',12);

	figure;
	subplot(2,1,1);
	imagesc(f,z,abs(mm));
	xlabel('frequency (Hz)');
	ylabel('z (cm)');
	title('Resulting pattern (linear regime prediction)');

	subplot(2,1,2);
	imagesc(f,z,angle(mm));
	xlabel('frequency (Hz)');
	ylabel('z (cm)');
	title('Resulting phase pattern (linear regime prediction)');
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%Camera  Camera superclass
%
% An abstract superclass for Toolbox camera classes.
%
% Methods::
%
% plot           plot projection of world point to image plane
% hold           control figure hold for image plane window
% ishold         test figure hold for image plane
% clf            clear image plane
% figure         figure holding the image plane
% mesh           draw shape represented as a mesh
% point          draw homogeneous points on image plane
% homline        draw homogeneous lines on image plane
% lineseg        draw line segment defined by points
% plot_camera    draw camera in world view
%-
% rpy            set camera attitude
% move           clone Camera after motion
% centre         get world coordinate of camera centre
%-
% delete         object destructor
% char           convert camera parameters to string
% display        display camera parameters
%-
% Properties (read/write)::
% npix    image dimensions (2x1)
% pp      principal point (2x1)
% rho     pixel dimensions (2x1) in metres
% T       camera pose as homogeneous transformation
%
% Properties (read only)::
% nu    number of pixels in u-direction
% nv    number of pixels in v-direction
% u0    principal point u-coordinate
% v0    principal point v-coordinate
%
% Notes::
%  - Camera is a reference object.
%  - Camera objects can be used in vectors and arrays
%  - This is an abstract class and must be subclassed and a project()
%    method defined.
%  - The object can create a window to display the Camera image plane, this
%    window is protected and can only be accessed by the plot methods of
%    this object.
%  - The project method is implemented by the concrete subclass.
%
% See also CentralCamera, SphericalCamera, FishEyeCamera, CatadiptricCamera.

% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

% TODO:
%   make a parent imaging class and subclass perspective, fisheye, panocam
%   test for points in front of camera and set to NaN if not
%   test for points off the image plane and set to NaN if not
%     make clipping/test flags
classdef Camera < handle

    properties
        name    % camera name
        type
        rho     % pixel dimensions 1x2
        pp      % principal point 1x2
        npix    % number of pixel 1x2
        T       % camera pose SE3 object
        noise   % pixel noise 1x2
        image
    end

    properties (SetAccess = protected)
        limits
        perspective
        h_imageplane     % handle for image plane
        h_3dview % handle for camera 3D view
        P           % world points (last plotted)
        holdon
        color
    end

    properties (Dependent = true, SetAccess = protected)
        u0
        v0
        nu
        nv
    end

    methods (Abstract)
        p = project(c, P, varargin);
    end

    methods

        function c = Camera(varargin)
            %Camera.Camera Create camera object
            %
            % Constructor for abstact Camera class, used by all subclasses.
            %
            % C = Camera(OPTIONS) creates a default (abstract) camera with null parameters.
            %
            % Options::
            % 'name',N          Name of camera
            % 'image',IM        Load image IM to image plane
            % 'resolution',N    Image plane resolution: NxN or N=[W H]
            % 'sensor',S        Image sensor size in metres (2x1) [metres]
            % 'centre',P        Principal point (2x1)
            % 'pixel',S         Pixel size: SxS or S=[W H]
            % 'noise',SIGMA     Standard deviation of additive Gaussian noise added to
            %                   returned image projections
            % 'pose',T          Pose of the camera as a homogeneous transformation
            % 'color',C         Color of image plane background (default [1 1 0.8])
            %
            % Notes::
            % - Normally the class plots points and lines into a set of axes that represent
            %   the image plane.  The 'image' option paints the specified image onto the
            %   image plane and allows points and lines to be overlaid.
            %
            % See also CentralCamera, FisheyeCamera, CatadioptricCamera, SphericalCamera.

            % default values
            c.type = '**abstract**';
            c.T = eye(4,4);
            c.pp = [];
            c.limits = [-1 1 -1 1];
            c.perspective = false;
            c.h_imageplane = [];
            c.h_3dview = [];
            c.holdon = false;
            c.color = [1 1 0.8];

            if nargin == 1 && isa(varargin{1}, 'Camera')
                return;
            else
                opt.name = 'noname';
                opt.image = [];
                opt.resolution = [];
                opt.centre = [];
                opt.sensor = [];
                opt.pixel = [1 1];
                opt.noise = [];
                opt.pose = [];
                opt.color = [];
                c.pp = [0 0];

                [opt,args] = tb_optparse(opt, varargin);

                c.name = opt.name;
                if ~isempty(opt.image)
                    c.image = opt.image;
                    c.npix = [size(c.image,2) size(c.image,1)];
                end
                if ~isempty(opt.resolution)
                    if length(opt.resolution) == 1
                        c.npix = [opt.resolution opt.resolution];
                    elseif length(opt.resolution) == 2
                        c.npix = opt.resolution;
                    else
                        error('resolution must be a 1- or 2-vector');
                    end
                end

                c.pp = opt.centre;
                c.rho = opt.pixel;
                if ~isempty(opt.color)
                    c.color = opt.color;
                end
                if ~isempty(opt.noise)
                    if length(opt.noise) == 1
                        c.noise = [opt.noise opt.noise];
                    elseif length(opt.noise) == 2
                        c.noise = opt.noise;
                    else
                        error('noise must be a 1- or 2-vector');
                    end
                end
                if ~isempty(opt.pose)
                    c.T = SE3(opt.pose);
                end
                if ~isempty(opt.sensor)
                    c.rho = opt.sensor ./ c.npix;
                end
            end

            if length(c.rho) == 1
                c.rho = ones(1,2) * c.rho;
            end
            if isempty(c.pp)
                fprintf('principal point not specified, setting it to centre of image plane\n');
                c.pp = c.npix / 2;
            end
        end

        function delete(c)
            %Camera.delete Camera object destructor
            %
            % C.delete() destroys all figures associated with the Camera object and
            % removes the object.
            %disp('delete camera object');
            if ~isempty(c.h_imageplane) && isgraphics(c.h_imageplane)
                delete(get(c.h_imageplane, 'Parent'));
            end
            if ~isempty(c.h_3dview) && isgraphics(c.h_3dview)
                delete(get(c.h_3dview, 'Parent'));
            end
        end

        function display(c)
            %Camera.display Display value
            %
            % C.display() displays a compact human-readable representation of the camera
            % parameters.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a Camera object and the command has no trailing
            %   semicolon.
            %
            % See also Camera.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            if loose
                disp(' ');
            end
            disp(char(c))
            if loose
                disp(' ');
            end
        end

        function s = char(c, s)
            %Camera.char Convert to string
            %
            % S = C.char() is a compact string representation of the camera parameters.

            s = '';
            if ~isempty(c.rho)
                s = strvcat(s, sprintf('  pixel size:     (%.4g, %.4g)', c.rho(1), c.rho(2)));
            end
            if ~isempty(c.pp)
                s = strvcat(s, sprintf('  principal pt:   (%.4g, %.4g)', c.u0, c.v0));
            end
            if ~isempty(c.npix)
                s = strvcat(s, sprintf('  number pixels:  %d x %d', c.nu, c.nv));
            end
            if ~isempty(c.noise)
                s = strvcat(s, sprintf('  noise:          %.4g,%.4g pix', c.noise));
            end
            s = strvcat(s,     sprintf('  pose:           %s', print(c.T, 'camera')));
            % = strvcat(s, [repmat('      ', 4,1) num2str(c.T)]);
        end

        function rpy(c, roll, pitch, yaw)
            %Camera.rpy Set camera attitude
            %
            % C.rpy(R, P, Y) sets the camera attitude to the specified roll-pitch-yaw angles.
            %
            % C.rpy(RPY) as above but RPY=[R,P,Y].
            if nargin == 2,
                pitch = roll(2);
                yaw = roll(3);
                roll = roll(1);
            end
            c.T = SE3.Ry(yaw) * SE3.Rx(pitch) * SE3.Rz(roll);
        end

        function c = centre(c)
            %Camera.centre Get camera position
            %
            % P = C.centre() is the 3-dimensional position of the camera centre (3x1).

            c = c.T.t;
        end

        function ishold = hold(c, flag)
            %Camera.hold Control hold on image plane graphics
            %
            % C.hold() sets "hold on" for the camera's image plane.
            %
            % C.hold(H) hold mode is set on if H is true (or > 0), and off if
            % H is false (or 0).

            % cam.ishold(); set hold
            % cam.ishold(flag); set hold
            if nargin < 2
                flag = true;
            end
            c.holdon = flag;
            if flag
                set(c.h_imageplane, 'NextPlot', 'add');
            else
                set(c.h_imageplane, 'NextPlot', 'replacechildren');
            end
            if nargout > 0
                % i = cam.ishold(); test hold condition
                ishold = c.holdon;

            end
        end

        function v = ishold(c)
            %Camera.ishold Return image plane hold status
            %
            % H = C.ishold() returns true (1) if the camera's image plane is in hold mode,
            % otherwise false (0).
            v = c.holdon();
        end

        function clf(c, flag)
            %Camera.clf Clear the image plane
            %
            % C.clf() removes all graphics from the camera's image plane.
            h = c.h_imageplane;
            if ~isempty(h) && isgraphics(h)
                % remove all children of the figure
                children = get(h, 'Children');
                for child=children
                    delete(child)
                end

                % if the camera is displaying an image
                if ~isempty(c.image)
                    c.figure
                    idisp(c.image, 'nogui');
                    hold on
                end
            end
        end

        function h = figure(c)
            %Camera.figure Return figure handle
            %
            % H = C.figure() is the handle of the figure that contains the camera's
            % image plane graphics.
            fig  = get(c.h_imageplane, 'Parent');
            figure( fig );
            if nargout > 0
                h = fig;
            end
        end

        % Return the graphics handle for this camera's image plane
        % and create the graphics if it doesnt exist
        %
        function h = plot_create(c, hin)

            if ~isempty(c.image)
                % if this camera is created from an image, then display that image

                if isempty(c.h_imageplane) || ~isgraphics(c.h_imageplane)
                    idisp(c.image, 'nogui');
                    set(gcf, 'name', sprintf('%s(%s) - image plane', class(c), c.name));
                    set(gcf, 'MenuBar', 'none');
                    hold on
                    h = gca;
                    %title(h, ['CentralCamera image plane:' c.name);
                    c.h_imageplane = h;
                    set(gcf, 'HandleVisibility', 'off');
                    set(h, 'HandleVisibility', 'off');
                else
                    h = c.h_imageplane;
                end
                return;
            end

            if isgraphics(c.h_imageplane)
                % it already exists, just return the handle
                h = c.h_imageplane;
                return;
            end
            c.h_imageplane = [];  % handle is invalid

            if (nargin == 2) && isgraphics(hin)
                % draw camera in an existing axes
                disp('draw image plane in existing axes')
                h = hin;
                h.HandleVisibility = 'Off';
            else
                disp('creating new figure for camera')

                figure
                h = axes
                fig = get(h, 'Parent');
                disp('make axes');
                axis square
                set(fig, 'MenuBar', 'none');
                set(fig, 'Tag', 'camera');
                set(h, 'Color', c.color);
                set(fig, 'HandleVisibility', 'off');
                set(fig, 'name', sprintf('%s(%s) - image plane', class(c), c.name));
            end
            % create an axis for camera view
            set(h, 'XLim', c.limits(1:2), 'YLim', c.limits(3:4), ...
                'DataAspectRatio', [1 1 1], ...
                'Xgrid', 'on', 'Ygrid', 'on', ...
                'Ydir' , 'reverse', ...
                'NextPlot', 'add', ...
                'Tag', c.name ...
                );
            c.h_imageplane = h;       % keep this around
            c.labelaxes();

            title(h, c.name);
            %figure( fig );   % raise the camera view
            set(h, 'NextPlot', 'replacechildren');
        end

        function labelaxes(c)
            h = c.h_imageplane;

            if ~isempty(c.npix)
                xlabel(h, 'u (pixels)');
                ylabel(h, 'v (pixels)');
            else
                xlabel(h, 'x (m)');
                ylabel(h, 'y (m)');
            end
        end

        function v =  plot(c, points, varargin)
            %Camera.plot Plot points on image plane
            %
            % C.plot(P, OPTIONS) projects world points P (3xN) to the image plane and plots them.  If P is 2xN
            % the points are assumed to be image plane coordinates and are plotted directly.
            %
            % UV = C.plot(P) as above but returns the image plane coordinates UV (2xN).
            %
            % - If P has 3 dimensions (3xNxS) then it is considered a sequence of point sets and is
            %   displayed as an animation.
            %
            % C.plot(L, OPTIONS) projects the world lines represented by the
            % array of Plucker objects (1xN) to the image plane and plots them.
            %
            % LI = C.plot(L, OPTIONS) as above but returns an array (3xN) of
            % image plane lines in homogeneous form.
            %
            % Options::
            % 'objpose',T     Transform all points by the homogeneous transformation T before
            %                  projecting them to the camera image plane.
            % 'pose',T         Set the camera pose to the homogeneous transformation T before
            %                  projecting points to the camera image plane.  Overrides the current
            %                  camera pose C.T.
            % 'fps',N          Number of frames per second for point sequence display
            % 'sequence'       Annotate the points with their index
            % 'textcolor',C    Text color for annotation (default black)
            % 'textsize',S     Text size for annotation (default 12)
            % 'drawnow'        Execute MATLAB drawnow function
            %
            % Additional options are considered MATLAB line or marker style parameters and are passed
            % directly to plot.
            %
            % See also Camera.mesh, Camera.hold, Camera.clf, Plucker.

            opt.fps = 5;
            opt.sequence = false;
            opt.textcolor = 'k';
            opt.textsize = 12;
            opt.drawnow = false;

            [opt,arglist,ls] = tb_optparse(opt, varargin);

            % get handle for this camera image plane
            h = c.plot_create();

            if isa(points, 'Plucker')
                % plot lines

                % project 3D world lines using the class project() method
                uv = c.project(points, arglist{:});
                c.hold(true);
                for line=uv
                    c.homline(line);
                end
                c.hold(false)
            else
                % plot points
                nr = numrows(points);

                if nr == 3
                    % project 3D world points using the class project() method
                    uv = c.project(points, arglist{:});
                else
                    uv = points;
                end

                if isempty(ls) || ~any(cellfun(@(x) ischar(x)&&contains(x, 'Marker'), arglist))
                    % set default style if none given
                    ls = {'Marker', 'o', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k', 'LineStyle', 'none'};
                else
                    ls = {ls};
                end

                for i=1:size(uv,3)
                    % for every frame in the animation sequence
                    plot(uv(1,:,i), uv(2,:,i), ls{:}, 'Parent', h);
                    if opt.sequence
                        for j=1:size(uv,2)
                            text(uv(1,j,i), uv(2,j,i), sprintf('  %d', j), ...
                                'HorizontalAlignment', 'left', ...
                                'VerticalAlignment', 'middle', ...
                                'FontUnits', 'pixels', ...
                                'FontSize', opt.textsize, ...
                                'Color', opt.textcolor, ...
                                'Parent', h);
                        end
                    end

                    if size(uv,3) > 1
                        pause(1/opt.fps);
                    end
                end
            end

            if opt.drawnow
                drawnow
            end

            if nargout > 0,
                v = uv;
            end
        end % plot

        function mesh(c, X, Y, Z, varargin)
            %Camera.mesh Plot mesh object on image plane
            %
            % C.mesh(X, Y, Z, OPTIONS) projects a 3D shape defined by the matrices X, Y, Z
            % to the image plane and plots them.  The matrices X, Y, Z are of the same size
            % and the corresponding elements of the matrices define 3D points.
            %
            % Options::
            % 'objpose',T   Transform all points by the homogeneous transformation T before
            %               projecting them to the camera image plane.
            % 'pose',T      Set the camera pose to the homogeneous transformation T before
            %               projecting points to the camera image plane.  Temporarily overrides
            %               the current camera pose C.T.
            %
            % Additional arguments are passed to plot as line style parameters.
            %
            % See also MESH, CYLINDER, SPHERE, MKCUBE, Camera.plot, Camera.hold, Camera.clf.

            % check that mesh matrices conform
            assert( all(size(X) == size(Y)) && all(size(X) == size(Z)), 'matrices must be the same size');

            opt.objpose = [];
            opt.pose = [];

            [opt,arglist,ls] = tb_optparse(opt, varargin);
            if isempty(opt.pose)
                opt.pose = c.T;
            end

            if isempty(ls)
                ls = {'k'};
            end

            % get handle for this camera image plane
            h = c.plot_create();

            % draw 3D line segments
            nsteps = 21;

            c.clf
            holdon = c.hold(1);
            s = linspace(0, 1, nsteps);

            for i=1:numrows(X)-1
                for j=1:numcols(X)-1
                    P0 = [X(i,j), Y(i,j), Z(i,j)]';
                    P1 = [X(i+1,j), Y(i+1,j), Z(i+1,j)]';
                    P2 = [X(i,j+1), Y(i,j+1), Z(i,j+1)]';

                    if c.perspective
                        % straight world lines are straight on the image plane
                        uv = c.project([P0 P1], 'setopt', opt);
                    else
                        % straight world lines are not straight, plot them piecewise
                        P = bsxfun(@times, (1-s), P0) + bsxfun(@times, s, P1);
                        uv = c.project(P, 'setopt', opt);
                    end
                    plot(uv(1,:)', uv(2,:)', ls{:}, arglist{:}, 'Parent', c.h_imageplane);

                    if c.perspective
                        % straight world lines are straight on the image plane
                        uv = c.project([P0 P2], 'setopt', opt);
                    else
                        % straight world lines are not straight, plot them piecewise
                        P = bsxfun(@times, (1-s), P0) + bsxfun(@times, s, P2);
                        uv = c.project(P, 'setopt', opt);
                    end
                    plot(uv(1,:)', uv(2,:)', ls{:}, arglist{:}, 'Parent', c.h_imageplane);
                end
            end

            for j=1:numcols(X)-1
                P0 = [X(end,j), Y(end,j), Z(end,j)]';
                P1 = [X(end,j+1), Y(end,j+1), Z(end,j+1)]';

                if c.perspective
                    % straight world lines are straight on the image plane
                    uv = c.project([P0 P1], 'setopt', opt);
                else
                    % straight world lines are not straight, plot them piecewise
                    P = bsxfun(@times, (1-s), P0) + bsxfun(@times, s, P1);
                    uv = c.project(P, 'setopt', opt);
                end
                plot(uv(1,:)', uv(2,:)', ls{:}, arglist{:}, 'Parent', c.h_imageplane);
            end
            c.hold(holdon); % turn hold off if it was initially off

        end % mesh

        function h =  point(c, p, varargin)
            %Camera.point Plot homogeneous points on image plane
            %
            % C.point(P) plots points on the camera image plane which are defined by columns
            % of P (3xN) considered as points in homogeneous form.

            % get handle for this camera image plane
            h = c.create

            uv = e2h(p);
            h = plot(uv(1,:), uv(2,:), varargin{:});
        end % point

        function h =  homline(c, lines, varargin)
            %Camera.homline Plot homogeneous lines on image plane
            %
            % C.homline(L) plots lines on the camera image plane which are defined by columns
            % of L (3xN) considered as lines in homogeneous form: a.u + b.v + c = 0.

            % get handle for this camera image plane
            h = c.plot_create;
            xlim = get(h, 'XLim');
            ylim = get(h, 'YLim');

            if numel(lines) == 3
                lines = lines(:);
            end

            for l=lines
                if abs(l(1)/l(2)) > 1
                    % steeper than 45deg
                    x = (-l(3) - l(2)*ylim) / l(1);
                    h = plot(x, ylim, varargin{:}, 'Parent', c.h_imageplane);
                else
                    % less than 45deg
                    y = (-l(3) - l(1)*xlim) / l(2);

                    h = plot(xlim, y, varargin{:}, 'Parent', c.h_imageplane);
                end
            end
        end % line

        function h = lineseg(c, p0, p1, varargin)
            % get handle for this camera image plane
            h = c.plot_create
            c.hold(1)
            for i=1:numcols(p0)
                plot([p0(1,i) p1(1,i)], [p0(2,i) p1(2,i)], varargin{:}, 'Parent', c.h_imageplane);
            end
        end

        function newcam = move(cam, T)
            %Camera.move Instantiate displaced camera
            %
            % C2 = C.move(T) is a new camera object that is a clone of C but its pose
            % is displaced by the homogeneous transformation T with respect to the
            % current pose of C.
            newcam = CentralCamera(cam);
            newcam.T = newcam.T * SE3.convert(T);
        end

        function movedby(c, robot)
            robot.addlistener('Moved', @(src,data)cameramove_callback(src,data,c));

            function cameramove_callback(robot, event, camera)
                camera.T = robot.fkine(robot.q);
            end
        end

        % return components of principal point and image size
        function v = get.u0(c)
            v = c.pp(1);
        end

        function v = get.v0(c)
            v = c.pp(2);
        end

        function v = get.rho(c)
            v = c.rho;
        end

        function v = get.nu(c)
            v = c.npix(1);
        end

        function v = get.nv(c)
            v = c.npix(2);
        end

        function c = set.T(c, Tc)
            c.T = SE3(Tc);

            if ~isempty(c.h_3dview) && isgraphics(c.h_3dview)
                set(c.h_3dview, 'Matrix', c.T.T);
            end
        end

        function c = set.rho(c, sxy)
            if isempty(sxy)
                c.rho = sxy;
            elseif length(sxy) == 1
                c.rho = [sxy sxy];
            elseif length(sxy) == 2
                c.rho = sxy(:)';
            else
                error('need 1 or 2 scale elements');
            end
        end

        function c = set.pp(c, pp)
            if isempty(pp)
                c.pp = [];
            elseif length(pp) == 1
                c.pp = [pp pp];
            elseif length(pp) == 2
                c.pp = pp(:)';
            else
                error('need 1 or 2 pp elements');
            end
        end

        function c = set.npix(c, npix)
            if ~isempty(npix)
                if length(npix) == 1,
                    c.npix = [npix npix];
                elseif length(npix) == 2,
                    c.npix = npix(:)';
                else
                    error('need 1 or 2 npix elements');
                end
                c.limits = [0 c.npix(1) 0 c.npix(2)];
            end
        end

        function help(c)
            disp(' C.plot(P)     return image coordinates for world points  P');
            disp(' C.point(P)     return image coordinates for world points  P');
            disp(' C.line(P)     return image coordinates for world points  P');
            disp(' C.clf     return image coordinates for world points  P');
            disp(' C.hold     return image coordinates for world points  P');
            disp(' C.project(P)     return image coordinates for world points  P');
            disp(' C.project(P, Tobj)  return image coordinates for world points P ');
            disp(' C.project(P, To, Tcm)  return image coordinates for world points P ');
            disp(' transformed by T prior to projection');
        end
    end % methods

end % class

author: petercorke

repo: machinevision-toolbox-matlab


text:
N  = 24;                            % Prediction horizon (number of iterations)
Nu  = N;                            % Control horizon (number of iterations)
Q = [1 0 1 0 0];                    % State weights
R = 0;                              % du weights
Ru = 1;                             % u weights
LB = 0*ones(Nu,1);                  % Lower bound of control input
UB = 1*ones(Nu,1);                  % Upper bound of control input
LBdu = nan;                         % Lower bound of control input rate
UBdu = nan;                         % Upper bound of control input rate
LBo = 0;                            % Lower bound of output
UBo = nan;                          % Upper bound of output

author: eurika-kaiser

repo: SINDY-MPC


text:
function tests = test_ft_connectivity_corr

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivity_corr

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_rpt_chan_chan(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nrpt   = 10;
nchan  = 3;
ntime  = 100;
dimord = 'rpt_chan_chan';

input = zeros(nrpt, nchan, nchan);
for i=1:nrpt
  dat = 10 * randn(nchan, ntime);
  input(i,:,:) = cov(dat');
end

result = {};
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', false);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1:2);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1:3);

assert(all(diag(result{2})==1));
assert(all(size(result{2})==[3 3]));
assert(all(size(result{3})==[2 2]));
assert(all(size(result{4})==[1 1]));
assert(all(size(result{5})==[0 0]));

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_rpt_chan_chan_time(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nrpt   = 10;
nchan  = 3;
nseg   = 20;  % sliding window, there are 20 segments per trial
ntime  = 50; % sliding window, there are 50 samples per segment
dimord = 'rpt_chan_chan_time';

input = zeros(nrpt, nchan, nchan);
for rpt=1:nrpt
  for seg=1:nseg
    dat = 10 * randn(nchan, ntime);
    input(rpt,:,:,seg) = cov(dat');
  end
end

result = {};
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', false);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1, 'allchanindx', 1:3);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 2, 'allchanindx', 1:3);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1:2, 'allchanindx', 1:3);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_rpt_chan_chan_freq(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nrpt   = 10;
nchan  = 3;
nfreq  = 5;
dimord = 'rpt_chan_chan_freq';

input = zeros(nrpt, nchan, nchan, nfreq);
for rpt=1:nrpt
  for freq=1:nfreq
    fdat = 10 * randn(nchan, 1) + 1i * 10 * randn(nchan, 1);
    input(rpt,:,:,freq) = fdat * ctranspose(fdat); % compute the cross-spectral density
  end
end

result = {};
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', false);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1, 'allchanindx', 1:3);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 2, 'allchanindx', 1:3);
result{end+1} = ft_connectivity_corr(input, 'dimord', dimord, 'hasjack', false, 'pownorm', true, 'pchanindx', 1:2, 'allchanindx', 1:3);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ y, m, d ] = day_borrow_hebrew ( y, m, d )

%*****************************************************************************80
%
%% DAY_BORROW_HEBREW borrows days from months in a Hebrew date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, a year, month, and day
%    representing a date.  On input, D might be negative.  On output,
%    M should have decreased by one month, and D gone up by the
%    number of days in the month we "cashed in".  Y may be affected
%    if the input value of M was 1.
%
  while ( d <= 0 )

    m = m - 1;

    [ y, m ] = month_borrow_hebrew ( y, m );

    days = month_length_hebrew ( y, m );

    d = d + days;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout=sigstar(groups,stats,nosort)
    % sigstar - Add significance stars to bar charts, boxplots, line charts, etc,
    %
    % H = sigstar(groups,stats,nsort)
    %
    % Purpose
    % Add stars and lines highlighting significant differences between pairs of groups.
    % The user specifies the groups and associated p-values. The function handles much of
    % the placement and drawing of the highlighting. Stars are drawn according to:
    %   * represents p<=0.05
    %  ** represents p<=1E-2
    % *** represents p<=1E-3
    %
    %
    % Inputs
    % groups - a cell array defining the pairs of groups to compare. Groups defined
    %          either as pairs of scalars indicating locations along the X axis or as
    %          strings corresponding to X-tick labels. Groups can be a mixture of both
    %          definition types.
    % stats -  a vector of p-values the same length as groups. If empty or missing it's
    %          assumed to be a vector of 0.05s the same length as groups. Nans are treated
    %          as indicating non-significance.
    % nsort -  optional, 0 by default. If 1, then significance markers are plotted in
    %          the order found in groups. If 0, then they're sorted by the length of the
    %          bar.
    %
    % Outputs
    % H - optionally return handles for significance highlights. Each row is a different
    %     highlight bar. The first column is the line. The second column is the text (stars).
    %
    %
    % Examples
    % 1.
    % bar([5,2,1.5])
    % sigstar({[1,2], [1,3]})
    %
    % 2.
    % bar([5,2,1.5])
    % sigstar({[2,3],[1,2], [1,3]},[nan,0.05,0.05])
    %
    % 3.  **DOESN'T WORK IN 2014b**
    % R=randn(30,2);
    % R(:,1)=R(:,1)+3;
    % boxplot(R)
    % set(gca,'XTick',1:2,'XTickLabel',{'A','B'})
    % H=sigstar({{'A','B'}},0.01);
    % ylim([-3,6.5])
    % set(H,'color','r')
    %
    % 4. Note the difference in the order with which we define the groups in the
    %    following two cases.
    % x=[1,2,3,2,1];
    % subplot(1,2,1)
    % bar(x)
    % sigstar({[1,2], [2,3], [4,5]})
    % subplot(1,2,2)
    % bar(x)
    % sigstar({[2,3],[1,2], [4,5]})
    %
    % ALSO SEE: demo_sigstar
    %
    % KNOWN ISSUES:
    % 1. Algorithm for identifying whether significance bar will overlap with
    %    existing plot elements may not work in some cases (see line 277)
    % 2. Bars may not look good on exported graphics with small page sizes.
    %    Simply increasing the width and height of the graph with the
    %    PaperPosition property of the current figure should fix things.
    %
    % Rob Campbell - CSHL 2013

    %Input argument error checking

    %If the user entered just one group pair and forgot to wrap it in a cell array
    %then we'll go easy on them and wrap it here rather then generate an error
    if ~iscell(groups) & length(groups)==2
        groups={groups};
    end

    if nargin<2
        stats=repmat(0.05,1,length(groups));
    end
    if isempty(stats)
        stats=repmat(0.05,1,length(groups));
    end
    if nargin<3
        nosort=0;
    end

    %Check the inputs are of the right sort
    if ~iscell(groups)
        error('groups must be a cell array')
    end

    if ~isvector(stats)
        error('stats must be a vector')
    end

    if length(stats)~=length(groups)
        error('groups and stats must be the same length')
    end

    %Each member of the cell array groups may be one of three things:
    %1. A pair of indices.
    %2. A pair of strings (in cell array) referring to X-Tick labels
    %3. A cell array containing one index and one string
    %
    % For our function to run, we will need to convert all of these into pairs of
    % indices. Here we loop through groups and do this.

    xlocs=nan(length(groups),2); %matrix that will store the indices
    xtl=get(gca,'XTickLabel');

    for ii=1:length(groups)
        grp=groups{ii};

        if isnumeric(grp)
            xlocs(ii,:)=grp; %Just store the indices if they're the right format already

        elseif iscell(grp) %Handle string pairs or string/index pairs

            if isstr(grp{1})
                a=strmatch(grp{1},xtl);
            elseif isnumeric(grp{1})
                a=grp{1};
            end
            if isstr(grp{2})
                b=strmatch(grp{2},xtl);
            elseif isnumeric(grp{2})
                b=grp{2};
            end

            xlocs(ii,:)=[a,b];
        end

        %Ensure that the first column is always smaller number than the second
        xlocs(ii,:)=sort(xlocs(ii,:));

    end

    %If there are any NaNs we have messed up.
    if any(isnan(xlocs(:)))
        error('Some groups were not found')
    end

    %Optionally sort sig bars from shortest to longest so we plot the shorter ones first
    %in the loop below. Usually this will result in the neatest plot. If we waned to
    %optimise the order the sig bars are plotted to produce the neatest plot, then this
    %is where we'd do it. Not really worth the effort, though, as few plots are complicated
    %enough to need this and the user can define the order very easily at the command line.
    if ~nosort
        [~,ind]=sort(xlocs(:,2)-xlocs(:,1),'ascend');
        xlocs=xlocs(ind,:);groups=groups(ind);
        stats=stats(ind);
    end

    %-----------------------------------------------------
    %Add the sig bar lines and asterisks
    holdstate=ishold;
    hold on

    H=ones(length(groups),2); %The handles will be stored here

    y=ylim;
    yd=myRange(y)*0.05; %separate sig bars vertically by 5%

    for ii=1:length(groups)
        thisY=findMinY(xlocs(ii,:))+yd;
        H(ii,:)=makeSignificanceBar(xlocs(ii,:),thisY,stats(ii));
    end
    %-----------------------------------------------------

    %Now we can add the little downward ticks on the ends of each line. We are
    %being extra cautious and leaving this it to the end just in case the y limits
    %of the graph have changed as we add the highlights. The ticks are set as a
    %proportion of the y axis range and we want them all to be the same the same
    %for all bars.
    yd=myRange(ylim)*0.01; %Ticks are 1% of the y axis range
    for ii=1:length(groups)
        y=get(H(ii,1),'YData');
        y(1)=y(1)-yd;
        y(4)=y(4)-yd;
        set(H(ii,1),'YData',y)
    end

    %Be neat and return hold state to whatever it was before we started
    if ~holdstate
        hold off
    elseif holdstate
        hold on
    end

    %Optionally return the handles to the plotted significance bars (first column of H)
    %and asterisks (second column of H).
    if nargout>0
        varargout{1}=H;
    end

end %close sigstar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Internal functions

function H=makeSignificanceBar(x,y,p)
    %makeSignificanceBar produces the bar and defines how many asterisks we get for a
    %given p-value

    if p<=1E-3
        stars='***';
    elseif p<=1E-2
        stars='**';
    elseif p<=0.05
        stars='*';
    elseif isnan(p)
        stars='n.s.';
    else
        stars='';
    end

    x=repmat(x,2,1);
    y=repmat(y,4,1);

    H(1)=plot(x(:),y,'-k','LineWidth',3,'Tag','sigstar_bar');

    %Increase offset between line and text if we will print "n.s."
    %instead of a star.
    if ~isnan(p)
        offset=0.005;
    else
        offset=0.02;
    end

    starY=mean(y)+myRange(ylim)*offset;
    H(2)=text(mean(x(:)),starY,stars,...
        'HorizontalAlignment','Center',...
        'BackGroundColor','none',...
        'Tag','sigstar_stars');

    Y=ylim;
    if Y(2)<starY
        ylim([Y(1),starY+myRange(Y)*0.05])
    end

end %close makeSignificanceBar

function Y=findMinY(x)
    % The significance bar needs to be plotted a reasonable distance above all the data points
    % found over a particular range of X values. So we need to find these data and calculat the
    % the minimum y value needed to clear all the plotted data present over this given range of
    % x values.
    %
    % This version of the function is a fix from Evan Remington
    oldXLim = get(gca,'XLim');
    oldYLim = get(gca,'YLim');

    axis(gca,'tight')
    set(gca,'xlim',x) %Matlab automatically re-tightens y-axis

    yLim = get(gca,'YLim'); %Now have max y value of all elements within range.
    Y = max(yLim);

    axis(gca,'normal')
    set(gca,'XLim',oldXLim,'YLim',oldYLim)

end %close findMinY

function rng=myRange(x)
    %replacement for stats toolbox range function
    rng = max(x) - min(x);
end %close myRange

author: mvallieres

repo: radiomics


text:
% DURATION_FEATURE Calculate the log-duration of an object
%
% Usage
%    duration = DURATION_FEATURE(x, object)
%
% Input
%    x (numeric): The file data (not used).
%    object (struct): The objects contained in the data.
%
% Output
%    duration (numeric): The log-duration of the objects.
%
% See also
%    PREPARE_DATABASE

function t = duration_feature(x, object)
	t = permute(log([object.u2]-[object.u1]+1),[1 3 2]);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pr01.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function summe=pr01(name,x)
% prints a (0,1) profile of x and returns the number of nonzeros
%
function summe=pr01(name,x)

% check for row or column vector
[m,n]=size(x);row=(m==1);col=(n==1);n=max(m,n);
if ~(row|col), error('x must be a vector'); end;

text=[name,': '];summe=0;
for k=1:n,
  if x(k),
    text=[text,'1'];
    summe=summe+1;
  else
    text=[text,'0'];
  end;
end;
disp([text,'   ',num2str(summe),' nonzeros']);

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [f,J] = spm_fx_phase (phi,u,P,M)
% State equation for a phase-coupled oscillator
% FORMAT [f,J] = spm_fx_phase (phi,u,P,M)
%
% phi       state variable
% u         []
% P         model (variable) parameter structure
% M         model (fixed) parameter structure
%
% f         Flow vector, dphi/dt
% J         Jacobian, J(i,j)=df_i/dphi_j
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_fx_phase.m 2908 2009-03-20 14:54:03Z will $

% Sin terms
if isfield(P,'As')
    Nr=size(P.As,1);
    Ns=size(P.As,3);

    % negative abs
    Ahs=~(P.As==0);
    As=-abs(P.As).*Ahs;
else
    Ns=0;
end

% Cos terms
if isfield(P,'Ac')
    Nc=size(P.Ac,3);
    Nr=size(P.Ac,1);

    % Positive abs
    Ahc=~(P.Ac==0);
    Ac=abs(P.Ac).*Ahc;
else
    Nc=0;
end

% Flow vector
f=zeros(Nr,1);
for i=1:Nr,
    change=0;
    for j=1:Nr,
        phi_diff=phi(i)-phi(j);
        for n=1:Ns,
            change=change+As(i,j,n)*sin(n*phi_diff);
        end
        for n=1:Nc,
            change=change+Ac(i,j,n)*cos(n*phi_diff);
        end
    end
    exc=M.freq+P.df(i)+change;
    f(i)=2*pi*exc;
end

if nargout == 1; return, end

% Jacobian
As=-As;
for i=1:Nr,
    for j=1:Nr,
        jac=0;
        if i==j
            % Diagonal
            for jj=1:Nr,
                for n=1:Ns,
                    jac=jac-n*As(i,jj,n);
                end
            end
        else
            % Off-diagonal
            phi_diff=phi(i)-phi(j);
            for n=1:Ns,
                jac=jac+n*As(i,j,n)*cos(n*phi_diff);
            end
            for n=1:Nc,
                jac=jac+n*Ac(i,j,n)*sin(n*phi_diff);
            end
        end
        J(i,j)=jac;
    end
end
J=J*2*pi;

author: spm


text:
function DI = detectChange(obj, t1, t2)
% Use logarithm for ease of the following thresholding
DI = abs(log(double(t1)+eps) - log((double(t2) + eps)));
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
global_fp = 0;
global_fn = 0;
global_tn = 0;
global_tp = 0;
sample = 90;

for i = 1:sample
    close all;
    load(['Detection/img' num2str(i) '/img' num2str(i) '_detection.mat']);
    img = rgb2gray(im2single(imread(['Detection/img' num2str(i) '/img' num2str(i) '.bmp'])));

    % -------------------------------------------------------------------------
    % Process the ground true data
    % -------------------------------------------------------------------------
    temp = zeros(500);
    for j=1:size(detection, 1)
        temp(max(floor(detection(j, 2)), 1), max(floor(detection(j, 1)), 1)) = 1;
    end
    neg = ~imdilate(temp, strel('disk', 5, 0)) ; % draw a circle around the pixel
    pos = ~neg;

%     Plot the ground true
    figure('Name',['Image: ' num2str(i)]) ; clf ;
    subplot(1,3,1) ; imagesc(img) ; axis equal ; title('image') ;
    hold on;
    plot(detection(:, 1), detection(:, 2), 's', 'MarkerSize',10, 'Color', 'g');
    subplot(1,3,2) ; imagesc(pos) ; axis equal ; title('positive points (blob centres)') ;
    subplot(1,3,3) ; imagesc(neg) ; axis equal ; title('negative points (not a blob)') ;
    colormap gray ;

    % -------------------------------------------------------------------------
    % Detection code
    % -------------------------------------------------------------------------
    min_blob_size = 40;

    img = vl_imsmooth(img,2); % Blob pixels would melt together
    mask = img < 0.6; % Leave only the darker blob

    % Label blob
    mask=imfill(mask,'holes');
    map=bwlabel(mask);
    labels = setdiff(unique(map),0)';

    result = zeros(size(img,1),size(img,2));

    % Filter blob size
    for j=labels

        blob_pixel = (map==j);
        blob_area = sum(sum(blob_pixel));

        if blob_area > min_blob_size
            result = result | blob_pixel;
        end

    end

    % -------------------------------------------------------------------------
    % Measurement
    % -------------------------------------------------------------------------
    y = zeros(size(pos),'single') ;
    y(pos) = +1 ;
    y(neg) = -1 ;

    fp = result > 0 & y < 0 ;
    fn = result < 1 & y > 0 ;
    tn = result <= 0 & y < 0 ;
    tp = result >= 1 & y > 0 ;

    global_fp = global_fp + sum(sum(fp>0));
    global_fn = global_fn + sum(sum(fn>0));
    global_tn = global_tn + sum(sum(tn>0));
    global_tp = global_tp + sum(sum(tp>0));

end

recall = global_tp/(global_tp+global_fn);
precision = global_tp/(global_tp+global_fp);

f1_score = 2*((precision * recall)/(precision + recall));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ unaryTermsNew, pairwiseTermsNew, energyConstant ] = projectEnergyBinaryPairwise( unaryTerms, pairwiseTerms, partialLabels )
%projectEnergyBinaryPairwise assigns spesified values to the subset of variables and constructs the energy from the unlabelled ones

maskUnlabeled = partialLabels == 0;
numUnlabeled = sum(maskUnlabeled);
numNodes = size( unaryTerms, 1 );
numEdges = size( pairwiseTerms, 1);
numLabels = 2;

maskEdgesNew = maskUnlabeled( pairwiseTerms(:, 1) ) & maskUnlabeled( pairwiseTerms(:, 2) );

pairwiseTermsNew = pairwiseTerms(maskEdgesNew, :);
newIds = nan(numUnlabeled, 1);
newIds(maskUnlabeled) = 1 : numUnlabeled;
pairwiseTermsNew(:, 1) = newIds( pairwiseTermsNew(:, 1) );
pairwiseTermsNew(:, 2) = newIds( pairwiseTermsNew(:, 2) );

energyConstant = 0;
unaryTermsNew = unaryTerms(maskUnlabeled, :);
labeledUnary = unaryTerms(~maskUnlabeled, :);
goodLabels = partialLabels(~maskUnlabeled);
energyConstant = energyConstant + ...
    sum( labeledUnary( (1 : numNodes - numUnlabeled)' + (goodLabels - 1) * (numNodes - numUnlabeled) ) );

labelMap = [ 3, 4; 5, 6];

for iEdge = 1 : numEdges
    node1 = pairwiseTerms(iEdge, 1);
    node2 = pairwiseTerms(iEdge, 2);
    label1 = partialLabels(node1);
    label2 = partialLabels(node2);
    if ~maskUnlabeled(node1) && ~maskUnlabeled(node2)
        energyConstant = energyConstant + pairwiseTerms( iEdge, labelMap(label1, label2) );
    elseif maskUnlabeled(node1) && ~maskUnlabeled(node2)
        unaryTermsNew( newIds(node1), : ) = unaryTermsNew( newIds(node1), : ) + ...
            reshape( pairwiseTerms(iEdge, labelMap(:, label2)), 1, numLabels);
    elseif ~maskUnlabeled(node1) && maskUnlabeled(node2)
        unaryTermsNew( newIds(node2), : ) = unaryTermsNew( newIds(node2), : ) + ...
            reshape( pairwiseTerms(iEdge, labelMap(label1, :)), 1, numLabels);
    end
end

end

author: aosokin


text:
function checkdiff(problem, x, d)
% Checks the consistency of the cost function and directional derivatives.
%
% function checkdiff(problem)
% function checkdiff(problem, x)
% function checkdiff(problem, x, d)
%
% checkdiff performs a numerical test to check that the directional
% derivatives defined in the problem structure agree up to first order with
% the cost function at some point x, along some direction d. The test is
% based on a truncated Taylor series (see online Manopt documentation).
%
% Both x and d are optional and will be sampled at random if omitted.
%
% See also: checkgradient checkhessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Verify that the problem description is sufficient.
    if ~canGetCost(problem)
        error('It seems no cost was provided.');
    end
    if ~canGetDirectionalDerivative(problem)
        error('It seems no directional derivatives were provided.');
    end

    x_isprovided = exist('x', 'var') && ~isempty(x);
    d_isprovided = exist('d', 'var') && ~isempty(d);

    if ~x_isprovided && d_isprovided
        error('If d is provided, x must be too, since d is tangent at x.');
    end

    % If x and / or d are not specified, pick them at random.
    if ~x_isprovided
        x = problem.M.rand();
    end
    if ~d_isprovided
        d = problem.M.randvec(x);
    end

    % Compute the value f0 at f and directional derivative at x along d.
    storedb = StoreDB();
    xkey = storedb.getNewKey();
    f0 = getCost(problem, x, storedb, xkey);
    df0 = getDirectionalDerivative(problem, x, d, storedb, xkey);

    % Compute the value of f at points on the geodesic (or approximation
    % of it) originating from x, along direction d, for stepsizes in a
    % large range given by h.
    h = logspace(-8, 0, 51);
    value = zeros(size(h));
    for i = 1 : length(h)
        y = problem.M.exp(x, d, h(i));
        ykey = storedb.getNewKey();
        value(i) = getCost(problem, y, storedb, ykey);
    end

    % Compute the linear approximation of the cost function using f0 and
    % df0 at the same points.
    model = polyval([df0 f0], h);

    % Compute the approximation error
    err = abs(model - value);

    % And plot it.
    loglog(h, err);
    title(sprintf(['Directional derivative check.\nThe slope of the '...
                   'continuous line should match that of the dashed '...
                   '(reference) line\nover at least a few orders of '...
                   'magnitude for h.']));
    xlabel('h');
    ylabel('Approximation error');

    line('xdata', [1e-8 1e0], 'ydata', [1e-8 1e8], ...
         'color', 'k', 'LineStyle', '--', ...
         'YLimInclude', 'off', 'XLimInclude', 'off');

    % In a numerically reasonable neighborhood, the error should decrease
    % as the square of the stepsize, i.e., in loglog scale, the error
    % should have a slope of 2.
    window_len = 10;
    [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
    hold on;
        loglog(h(range), 10.^polyval(poly, log10(h(range))), ...
               'r-', 'LineWidth', 3);
    hold off;

    fprintf('The slope should be 2. It appears to be: %g.\n', poly(1));
    fprintf(['If it is far from 2, then directional derivatives ' ...
             'might be erroneous.\n']);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function FPDF = GaussPDF
% Return function handles to routines to calculate the area, mean, and
% second moment of a unit-variance, zero-mean, Gaussian probability
% density function.
%                 b
%   Farea(a,b) = Int p(x) dx
%                 a
%                 b
%   Fmean(a,b) = Int x p(x) dx
%                 a
%                b
%   Fvar(a,b) = Int x^2 p(x) dx
%                a
% where p(x) = 1/sqrt(2 pi) e^(-x^2/2)

FPDF = {@Garea, @Gmean, @Gvar};

% ----- ------
function v = Garea (a, b)

% Evaluate the function so as to avoid taking differences
% between nearly equal quantities (e.g. when a<0 and b<0)
if (b >= 0)
  if (a >= 0)
    v = F0(b) - F0(a);                    % Both a and b positive
  else
    v = (F0(b) + 0.5) + (F0(-a) + 0.5);   % a negative, b positive
  end
else
  if (a < 0)
    v = F0(-a) - F0(-b);                  % Both a and b negative
  else
    v = (-0.5 - F0(a)) + (-0.5 - F0(-b)); % a positive, b negative
  end
end

return

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F0(b) - F0(a)
%
% F0(x) = c Int e^(-x^2/2) dx  [c = 1/sqrt(2 pi)]
%       = -Q(x)
%       = -1/2 erfc(x/sqrt(2))

function v = F0(x)

v = -0.5*erfc(x/sqrt(2));

return

% ----- ------
function v = Gmean (a, b)

% Since the integrand x*p(x) is odd, Gmean(A,B)=Gmean(abs(A),abs(B))
v = F1(abs(b)) - F1(abs(a));

return

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F1(b) - F1(a)
%
% F1(x) = c Int x e^(-x^2/2) du,
%
% A change of variables v=x^2/2, dv=x dx, gives
%
% F1(x) = c Int e^(-v) dv
%       = -c e^(-v)
%       = -c e^(-x^2/2)

function v = F1 (x)

v = -exp(-0.5*x^2)/sqrt(2*pi);

return

% ----- ------
function v = Gvar (a, b)

% Evaluate the function so as to avoid taking differences
% between nearly equal quantities (e.g. when a<0 and b<0)
if (b >= 0)
  if (a >= 0)
    v = F2(b) - F2(a);                    % Both a and b positive
  else
    v = (F2(b) + 0.5) + (F2(-a) + 0.5);   % a negative, b positive
  end
else
  if (a < 0)
    v = F2(-a) - F2(-b);                  % Both a and b negative
  else
    v = (-0.5 - F2(a)) + (-0.5 - F2(-b)); % a positive, b negative
  end
end

return

% -----
% Evaluate the indefinite integral
% The definite integral between a and b is F2(b) - F2(a)
%
% F2(x) = c Int x^2 e^(-x^2/2) dx
%
% Using integration by parts, identify
%    f(x) = x,           f'(x) = 1,
%   g'(u) = x e^(-x2/2),  g(x) = -e^(-x^2/2)
%
%  Int f(x) g'(x) dx = f(x) g(x) - Int f'(x) g(x) dx
%
% Then
% F2(x) = -c x e^(-x^2/2) + c Int e^(-x^2/2) dx
%       = -c x e^(-x^2/2) - Q(x)

function v = F2 (x)

% The function evaluation fails for x = Inf because the term x e^(-x^2/2)
% in the expression returns NaN. For non-infinite x, the expression
% evaluates to 0 for large x.
xmax = 40;
if (x > xmax)
  v = 0;
else
  v = -x*exp(-0.5*x^2)/sqrt(2*pi) - 0.5*erfc(x/sqrt(2));
end

return

author: Sable

repo: mcbench-benchmarks


text:
function Del = EnvironmentalSelection(Population,nSub)
% The environmental selection of IM-MOEA

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

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,nSub);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front by crowding distance
    Last     = find(FrontNo==MaxFNo);
    CrowdDis = CrowdingDistance(Population(Last).objs);
    [~,rank] = sort(CrowdDis,'descend');
    Next(Last(rank(1:nSub-sum(Next)))) = true;
    % The index of deleted solutions
    Del = ~Next;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function  [pl2] = PL2GetFileIndex(filename)
% PL2GetFileIndex(filename): gets pl2 file index. File indexes are cached.
%        If cached file index for the specified filename exists, returns cached file index,
%       otherwise, calls PL2ReadFileIndex(filename)
%
% [pl2] = PL2GetFileIndex(filename)
%
% INPUT:
%   filename - pl2 file path. If empty string, will clear all cached file indexes.
%
% OUTPUT:
%   pl2 - pl2 file index object (described in detail in PL2ReadFileIndex.m)
%       keeps pl2 file indexes cached in memory.
%       To clear indexes, use this command:
%       PL2GetFileIndex('')
%       or
%       clear PL2GetFileIndex
%

% this array of pl2 indexes is persisted between calls to this function
persistent PL2Indexes;

if nargin ~= 1
    error 'expected 1 input argument';
end

if isempty(filename)
    pl2 = 0;
    PL2Indexes={};
    return
end

if exist(filename,'file') ~= 2
    error('file does not exist');
end

if isempty(PL2Indexes)
    PL2Indexes={};
    PL2Indexes{1}.Index = PL2ReadFileIndex(filename);
    PL2Indexes{1}.Name = filename;
    pl2 = PL2Indexes{1}.Index;
    return;
else
    n = size(PL2Indexes,2);
    for i=1:n
        if strcmp(PL2Indexes{i}.Name, filename) == 1
            % we found Index with specified file Name
            if isfield(PL2Indexes{i}, 'Index')
                % if the index structure has pl2 index assigned, return it
                pl2 = PL2Indexes{i}.Index;
                return;
            else
                % read pl2 index from file, assign to persistent indexes and return
                PL2Indexes{i}.Index = PL2ReadFileIndex(filename);
                pl2 = PL2Indexes{i}.Index;
                return;
            end
        end
    end
    % we did not find persisted Index.
    % read the index from file
    PL2Indexes{n+1}.Index = PL2ReadFileIndex(filename);
    PL2Indexes{n+1}.Name = filename;
    pl2 = PL2Indexes{n+1}.Index;
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef TabPanel < uix.Container & uix.mixin.Panel
    %uix.TabPanel  Tab panel
    %
    %  p = uix.TabPanel(p1,v1,p2,v2,...) constructs a tab panel and sets
    %  parameter p1 to value v1, etc.
    %
    %  A tab panel shows one of its contents and hides the others according
    %  to which tab is selected.
    %
    %  From R2014b, MATLAB provides uitabgroup and uitab as standard
    %  components.  Consider using uitabgroup and uitab for new code if
    %  these meet your requirements.
    %
    %  See also: uitabgroup, uitab, uix.CardPanel

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1480 $ $Date: 2017-02-15 16:56:13 +0100 (Wed, 15 Feb 2017) $

    properties( Access = public, Dependent, AbortSet )
        FontAngle % font angle
        FontName % font name
        FontSize % font size
        FontWeight % font weight
        FontUnits % font weight
        ForegroundColor % tab text color [RGB]
        HighlightColor % border highlight color [RGB]
        ShadowColor % border shadow color [RGB]
    end

    properties
        SelectionChangedFcn = '' % selection change callback
    end

    properties( Access = public, Dependent, AbortSet )
        TabEnables % tab enable states
        TabLocation % tab location [top|bottom]
        TabTitles % tab titles
        TabContextMenus % tab context menus
        TabWidth % tab width
    end

    properties( Access = private )
        FontAngle_ = get( 0, 'DefaultUicontrolFontAngle' ) % backing for FontAngle
        FontName_ = get( 0, 'DefaultUicontrolFontName' ) % backing for FontName
        FontSize_ = get( 0, 'DefaultUicontrolFontSize' ) % backing for FontSize
        FontWeight_ = get( 0, 'DefaultUicontrolFontWeight' ) % backing for FontWeight
        FontUnits_ = get( 0, 'DefaultUicontrolFontUnits' ) % backing for FontUnits
        ForegroundColor_ = get( 0, 'DefaultUicontrolForegroundColor' ) % backing for ForegroundColor
        HighlightColor_ = [1 1 1] % backing for HighlightColor
        ShadowColor_ = [0.7 0.7 0.7] % backing for ShadowColor
        ParentBackgroundColor = get( 0, 'DefaultUicontrolForegroundColor' ) % default parent background color
        Tabs = gobjects( [0 1] ) % tabs
        TabListeners = event.listener.empty( [0 1] ) % tab listeners
        TabLocation_ = 'top' % backing for TabPosition
        TabHeight = -1 % cache of tab height (-1 denotes stale cache)
        TabWidth_ = 50 % backing for TabWidth
        Dividers % tab dividers
        BackgroundColorListener % listener
        SelectionChangedListener % listener
        ParentListener % listener
        ParentBackgroundColorListener % listener
    end

    properties( Access = private, Constant )
        FontNames = listfonts() % all available font names
        DividerMask = uix.TabPanel.getDividerMask() % divider image data
        DividerWidth = 8 % divider width
        TabMinimumHeight = 9 % tab minimum height
        Tint = 0.85 % tint factor for unselected tabs
    end

    methods

        function obj = TabPanel( varargin )
            %uix.TabPanel  Tab panel constructor
            %
            %  p = uix.TabPanel() constructs a tab panel.
            %
            %  p = uix.TabPanel(p1,v1,p2,v2,...) sets parameter p1 to value
            %  v1, etc.

            % Create dividers
            dividers = matlab.ui.control.UIControl( 'Internal', true, ...
                'Parent', obj, 'Units', 'pixels', 'Style', 'checkbox',...
                'Tag', 'TabPanelDividers' );

            % Create listeners
            backgroundColorListener = event.proplistener( obj, ...
                findprop( obj, 'BackgroundColor' ), 'PostSet', ...
                @obj.onBackgroundColorChange );
            selectionChangedListener = event.listener( obj, ...
                'SelectionChanged', @obj.onSelectionChanged );
            parentListener = event.proplistener( obj, ...
                findprop( obj, 'Parent' ), 'PostSet', ...
                @obj.onParentChanged );

            % Store properties
            obj.Dividers = dividers;
            obj.BackgroundColorListener = backgroundColorListener;
            obj.SelectionChangedListener = selectionChangedListener;
            obj.ParentListener = parentListener;

            % Set properties
            try
                uix.set( obj, varargin{:} )
            catch e
                delete( obj )
                e.throwAsCaller()
            end

        end % constructor

    end % structors

    methods

        function value = get.FontAngle( obj )

            value = obj.FontAngle_;

        end % get.FontAngle

        function set.FontAngle( obj, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'normal','italic','oblique'} ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''FontAngle'' must be ''normal'', ''italic'' or ''oblique''.' )

            % Set
            obj.FontAngle_ = value;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.FontAngle = value;
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

        end % set.FontAngle

        function value = get.FontName( obj )

            value = obj.FontName_;

        end % get.FontName

        function set.FontName( obj, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, obj.FontNames ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''FontName'' must be a valid font name.' )

            % Set
            obj.FontName_ = value;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.FontName = value;
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

        end % set.FontName

        function value = get.FontSize( obj )

            value = obj.FontSize_;

        end % get.FontSize

        function set.FontSize( obj, value )

            % Check
            assert( isa( value, 'double' ) && isscalar( value ) && ...
                isreal( value ) && ~isinf( value ) && ...
                ~isnan( value ) && value > 0, ...
                'uix:InvalidPropertyValue', ...
                'Property ''FontSize'' must be a positive scalar.' )

            % Set
            obj.FontSize_ = value;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.FontSize = value;
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

        end % set.FontSize

        function value = get.FontWeight( obj )

            value = obj.FontWeight_;

        end % get.FontWeight

        function set.FontWeight( obj, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'normal','bold'} ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''FontWeight'' must be ''normal'' or ''bold''.' )

            % Set
            obj.FontWeight_ = value;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.FontWeight = value;
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

        end % set.FontWeight

        function value = get.FontUnits( obj )

            value = obj.FontUnits_;

        end % get.FontUnits

        function set.FontUnits( obj, value )

            % Check
            assert( ischar( value ) && ...
                any( strcmp( value, {'inches','centimeters','points','pixels'} ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''FontUnits'' must be ''inches'', ''centimeters'', ''points'' or ''pixels''.' )

            % Compute size in new units
            oldUnits = obj.FontUnits_;
            oldSize = obj.FontSize_;
            newUnits = value;
            newSize = oldSize * convert( oldUnits ) / convert( newUnits );

            % Set size and units
            obj.FontSize_ = newSize;
            obj.FontUnits_ = newUnits;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.FontUnits = newUnits;
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

            function factor = convert( units )
                %convert  Compute conversion factor to points
                %
                %  f = convert(u) computes the conversion factor from units
                %  u to points.  For example, convert('inches') since 1
                %  inch equals 72 points.

                persistent SCREEN_PIXELS_PER_INCH
                if isequal( SCREEN_PIXELS_PER_INCH, [] ) % uninitialized
                    SCREEN_PIXELS_PER_INCH = get( 0, 'ScreenPixelsPerInch' );
                end

                switch units
                    case 'inches'
                        factor = 72;
                    case 'centimeters'
                        factor = 72 / 2.54;
                    case 'points'
                        factor = 1;
                    case 'pixels'
                        factor = 72 / SCREEN_PIXELS_PER_INCH;
                end

            end % convert

        end % set.FontUnits

        function value = get.ForegroundColor( obj )

            value = obj.ForegroundColor_;

        end % get.ForegroundColor

        function set.ForegroundColor( obj, value )

            % Check
            assert( isnumeric( value ) && isequal( size( value ), [1 3] ) && ...
                all( isreal( value ) ) && all( value >= 0 ) && all( value <= 1 ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''ForegroundColor'' must be an RGB triple.' )

            % Set
            obj.ForegroundColor_ = value;

            % Update existing tabs
            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tab = tabs(ii);
                tab.ForegroundColor = value;
            end

        end % set.ForegroundColor

        function value = get.HighlightColor( obj )

            value = obj.HighlightColor_;

        end % get.HighlightColor

        function set.HighlightColor( obj, value )

            % Check
            assert( isnumeric( value ) && isequal( size( value ), [1 3] ) && ...
                all( isreal( value ) ) && all( value >= 0 ) && all( value <= 1 ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''HighlightColor'' must be an RGB triple.' )

            % Set
            obj.HighlightColor_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.HighlightColor

        function set.SelectionChangedFcn( obj, value )

            % Check
            if ischar( value ) % string
                % OK
            elseif isa( value, 'function_handle' ) && ...
                    isequal( size( value ), [1 1] ) % function handle
                % OK
            elseif iscell( value ) && ndims( value ) == 2 && ...
                    size( value, 1 ) == 1 && size( value, 2 ) > 0 && ...
                    isa( value{1}, 'function_handle' ) && ...
                    isequal( size( value{1} ), [1 1] ) %#ok<ISMAT> % cell callback
                % OK
            else
                error( 'uix:InvalidPropertyValue', ...
                    'Property ''SelectionChangedFcn'' must be a valid callback.' )
            end

            % Set
            obj.SelectionChangedFcn = value;

        end % set.SelectionChangedFcn

        function value = get.ShadowColor( obj )

            value = obj.ShadowColor_;

        end % get.ShadowColor

        function set.ShadowColor( obj, value )

            % Check
            assert( isnumeric( value ) && isequal( size( value ), [1 3] ) && ...
                all( isreal( value ) ) && all( value >= 0 ) && all( value <= 1 ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''ShadowColor'' must be an RGB triple.' )

            % Set
            obj.ShadowColor_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.ShadowColor

        function value = get.TabEnables( obj )

            value = get( obj.Tabs, {'Enable'} );
            value(strcmp( value, 'inactive' )) = {'on'};

        end % get.TabEnables

        function set.TabEnables( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Retrieve tabs
            tabs = obj.Tabs;
            tabListeners = obj.TabListeners;

            % Check
            assert( iscellstr( value ) && ...
                isequal( size( value ), size( tabs ) ) && ...
                all( strcmp( value, 'on' ) | strcmp( value, 'off' ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''TabEnables'' should be a cell array of strings ''on'' or ''off'', one per tab.' )

            % Set
            tf = strcmp( value, 'on' );
            value(tf) = {'inactive'};
            for ii = 1:numel( tabs )
                tabs(ii).Enable = value{ii};
                tabListeners(ii).Enabled = tf(ii);
            end

            % Show selected child
            obj.showSelection()

            % Mark as dirty
            obj.Dirty = true;

        end % set.TabEnables

        function value = get.TabLocation( obj )

            value = obj.TabLocation_;

        end % get.TabLocation

        function set.TabLocation( obj, value )

            % Check
            assert( ischar( value ) && ...
                any( strcmp( value, {'top','bottom'} ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''TabLocation'' should be ''top'' or ''bottom''.' )

            % Set
            obj.TabLocation_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.TabLocation

        function value = get.TabTitles( obj )

            value = get( obj.Tabs, {'String'} );

        end % get.TabTitles

        function set.TabTitles( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Retrieve tabs
            tabs = obj.Tabs;

            % Check
            assert( iscellstr( value ) && ...
                isequal( size( value ), size( tabs ) ), ...
                'uix:InvalidPropertyValue', ...
                'Property ''TabTitles'' should be a cell array of strings, one per tab.' )

            % Set
            n = numel( tabs );
            for ii = 1:n
                tabs(ii).String = value{ii};
            end

            % Mark as dirty
            obj.TabHeight = -1;
            obj.Dirty = true;

        end % set.TabTitles

        function value = get.TabContextMenus( obj )

            tabs = obj.Tabs;
            n = numel( tabs );
            value = cell( [n 1] );
            for ii = 1:n
                value{ii} = tabs(ii).UIContextMenu;
            end

        end % get.TabContextMenus

        function set.TabContextMenus( obj, value )

            tabs = obj.Tabs;
            n = numel( tabs );
            for ii = 1:n
                tabs(ii).UIContextMenu = value{ii};
            end

        end % set.TabContextMenus

        function value = get.TabWidth( obj )

            value = obj.TabWidth_;

        end % get.TabWidth

        function set.TabWidth( obj, value )

            % Check
            assert( isa( value, 'double' ) && isscalar( value ) && ...
                isreal( value ) && ~isinf( value ) && ...
                ~isnan( value ) && value ~= 0, ...
                'uix:InvalidPropertyValue', ...
                'Property ''TabWidth'' must be a non-zero scalar.' )

            % Set
            obj.TabWidth_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.TabWidth

    end % accessors

    methods( Access = protected )

        function redraw( obj )

            % Compute positions
            bounds = hgconvertunits( ancestor( obj, 'figure' ), ...
                [0 0 1 1], 'normalized', 'pixels', obj );
            w = ceil( bounds(1) + bounds(3) ) - floor( bounds(1) ); % width
            h = ceil( bounds(2) + bounds(4) ) - floor( bounds(2) ); % height
            p = obj.Padding_; % padding
            tabs = obj.Tabs;
            n = numel( tabs ); % number of tabs
            tH = obj.TabHeight; % tab height
            if tH == -1 % cache stale, refresh
                if n > 0
                    cTabExtents = get( tabs, {'Extent'} );
                    tabExtents = vertcat( cTabExtents{:} );
                    tH = max( tabExtents(:,4) );
                end
                tH = max( tH, obj.TabMinimumHeight ); % apply minimum
                tH = ceil( tH ); % round up
                obj.TabHeight = tH; % store
            end
            cH = max( [h - 2 * p - tH, 1] ); % contents height
            switch obj.TabLocation_
                case 'top'
                    cY = 1 + p; % contents y
                    tY = cY + cH + p; % tab y
                case 'bottom'
                    tY = 1; % tab y
                    cY = tY + tH + p; % contents y
            end
            cX = 1 + p; % contents x
            cW = max( [w - 2 * p, 1] ); % contents width
            tW = obj.TabWidth_; % tab width
            dW = obj.DividerWidth; % tab divider width
            if tW < 0 && n > 0 % relative
                tW = max( ( w - (n+1) * dW ) / n, 1 );
            end
            tW = ceil( tW ); % round up
            for ii = 1:n
                tabs(ii).Position = [1 + (ii-1) * tW + ii * dW, tY, tW, tH];
            end
            obj.Dividers.Position = [0 tY w+1 tH];
            contentsPosition = [cX cY cW cH];

            % Redraw tabs
            obj.redrawTabs()

            % Redraw contents
            selection = obj.Selection_;
            if selection ~= 0 && strcmp( obj.TabEnables{selection}, 'on' )
                uix.setPosition( obj.Contents_(selection), contentsPosition, 'pixels' )
            end

        end % redraw

        function addChild( obj, child )
            %addChild  Add child
            %
            %  c.addChild(d) adds the child d to the container c.

            % Create new tab
            n = numel( obj.Tabs );
            tab = matlab.ui.control.UIControl( 'Internal', true, ...
                'Parent', obj, 'Style', 'text', 'Enable', 'inactive', ...
                'Units', 'pixels', 'FontUnits', obj.FontUnits_, ...
                'FontSize', obj.FontSize_, 'FontName', obj.FontName_, ...
                'FontAngle', obj.FontAngle_, 'FontWeight', obj.FontWeight_, ...
                'ForegroundColor', obj.ForegroundColor_, ...
                'String', sprintf( 'Page %d', n + 1 ) );
            tabListener = event.listener( tab, 'ButtonDown', @obj.onTabClicked );
            obj.Tabs(n+1,:) = tab;
            obj.TabListeners(n+1,:) = tabListener;

            % Mark as dirty
            obj.TabHeight = -1;

            % Check for bug
            if verLessThan( 'MATLAB', '8.5' ) && strcmp( child.Visible, 'off' )
                obj.G1218142 = true;
            end

            % Select new content
            oldSelection = obj.Selection_;
            if numel( obj.Contents_ ) == 0
                newSelection = 1;
                obj.Selection_ = newSelection;
            else
                newSelection = oldSelection;
            end

            % Call superclass method
            addChild@uix.mixin.Container( obj, child )

            % Show selected child
            obj.showSelection()

            % Notify selection change
            if oldSelection ~= newSelection
                obj.notify( 'SelectionChanged', ...
                    uix.SelectionData( oldSelection, newSelection ) )
            end

        end % addChild

        function removeChild( obj, child )
            %removeChild  Remove child
            %
            %  c.removeChild(d) removes the child d from the container c.

            % Find index of removed child
            contents = obj.Contents_;
            index = find( contents == child );

            % Remove tab
            delete( obj.Tabs(index) )
            obj.Tabs(index,:) = [];
            obj.TabListeners(index,:) = [];

            % Call superclass method
            removeChild@uix.mixin.Panel( obj, child )

        end % removeChild

        function reorder( obj, indices )
            %reorder  Reorder contents
            %
            %  c.reorder(i) reorders the container contents using indices
            %  i, c.Contents = c.Contents(i).

            % Reorder
            obj.Tabs = obj.Tabs(indices,:);
            obj.TabListeners = obj.TabListeners(indices,:);

            % Call superclass method
            reorder@uix.mixin.Panel( obj, indices )

        end % reorder

        function reparent( obj, oldFigure, newFigure )
            %reparent  Reparent container
            %
            %  c.reparent(a,b) reparents the container c from the figure a
            %  to the figure b.

            if ~isequal( oldFigure, newFigure )
                contextMenus = obj.TabContextMenus;
                for ii = 1:numel( contextMenus )
                    contextMenu = contextMenus{ii};
                    if ~isempty( contextMenu )
                        contextMenu.Parent = newFigure;
                    end
                end
            end

            % Call superclass method
            reparent@uix.mixin.Panel( obj, oldFigure, newFigure )

        end % reparent

        function showSelection( obj )
            %showSelection  Show selected child, hide the others
            %
            %  c.showSelection() shows the selected child of the container
            %  c, and hides the others.

            % Call superclass method
            showSelection@uix.mixin.Panel( obj )

            % If not enabled, hide selected contents too
            selection = obj.Selection_;
            if selection ~= 0 && strcmp( obj.TabEnables{selection}, 'off' )
                child = obj.Contents_(selection);
                child.Visible = 'off';
                if isa( child, 'matlab.graphics.axis.Axes' )
                    child.ContentsVisible = 'off';
                end
                % As a remedy for g1100294, move off-screen too
                margin = 1000;
                if isa( child, 'matlab.graphics.axis.Axes' ) ...
                        && strcmp(child.ActivePositionProperty, 'outerposition' )
                    child.OuterPosition(1) = -child.OuterPosition(3)-margin;
                else
                    child.Position(1) = -child.Position(3)-margin;
                end
            end

        end % showSelection

    end % template methods

    methods( Access = private )

        function redrawTabs( obj )
            %redrawTabs  Redraw tabs
            %
            %  p.redrawTabs() redraws the tabs.

            % Get relevant properties
            selection = obj.Selection_;
            tabs = obj.Tabs;
            t = numel( tabs );
            dividers = obj.Dividers;

            % Handle no tabs as a special case
            if t == 0
                dividers.Visible = 'off'; % hide
                return
            end

            % Repaint tabs
            backgroundColor = obj.BackgroundColor;
            for ii = 1:t
                tab = tabs(ii);
                if ii == selection
                    tab.BackgroundColor = backgroundColor;
                else
                    tab.BackgroundColor = obj.Tint * backgroundColor;
                end
            end

            % Repaint dividers
            d = t + 1;
            dividerNames = repmat( 'F', [d 2] ); % initialize
            dividerNames(1,1) = 'E'; % end
            dividerNames(end,2) = 'E'; % end
            if selection ~= 0
                dividerNames(selection,2) = 'T'; % selected
                dividerNames(selection+1,1) = 'T'; % selected
            end
            tH = obj.TabHeight;
            assert( tH >= obj.TabMinimumHeight, 'uix:InvalidState', ...
                'Cannot redraw tabs with invalid TabHeight.' )
            tW = obj.Tabs(1).Position(3);
            dW = obj.DividerWidth;
            allCData = zeros( [tH 0 3] ); % initialize
            map = [obj.ShadowColor; obj.BackgroundColor; ...
                obj.Tint * obj.BackgroundColor; obj.HighlightColor;...
                obj.ParentBackgroundColor];
            for ii = 1:d
                % Select mask
                iMask = obj.DividerMask.( dividerNames(ii,:) );
                % Resize
                iData = repmat( iMask(5,:), [tH 1] );
                iData(1:4,:) = iMask(1:4,:);
                iData(end-3:end,:) = iMask(end-3:end,:);
                % Convert to RGB
                cData = ind2rgb( iData+1, map );
                % Orient
                switch obj.TabLocation_
                    case 'bottom'
                        cData = flipud( cData );
                end
                % Insert
                allCData(1:tH,(ii-1)*(dW+tW)+(1:dW),:) = cData; % center
                if ii > 1 % extend left under transparent uicontrol edge
                    allCData(1:tH,(ii-1)*(dW+tW),:) = cData(:,1,:);
                end
                if ii < d % extend right under transparent uicontrol edge
                    allCData(1:tH,(ii-1)*(dW+tW)+dW+1,:) = cData(:,end,:);
                end
            end
            dividers.CData = allCData; % paint
            dividers.BackgroundColor = obj.ParentBackgroundColor;
            dividers.Visible = 'on'; % show

        end % redrawTabs

    end % helper methods

    methods( Access = private )

        function onTabClicked( obj, source, ~ )

            % Update selection
            oldSelection = obj.Selection_;
            newSelection = find( source == obj.Tabs );
            if oldSelection == newSelection, return, end % abort set
            obj.Selection_ = newSelection;

            % Show selected child
            obj.showSelection()

            % Mark as dirty
            obj.Dirty = true;

            % Notify selection change
            obj.notify( 'SelectionChanged', ...
                uix.SelectionData( oldSelection, newSelection ) )

        end % onTabClicked

        function onBackgroundColorChange( obj, ~, ~ )

            % Mark as dirty
            obj.Dirty = true;

        end % onBackgroundColorChange

        function onSelectionChanged( obj, source, eventData )

            % Call callback
            callback = obj.SelectionChangedFcn;
            if ischar( callback ) && isequal( callback, '' )
                % do nothing
            elseif ischar( callback )
                feval( callback, source, eventData )
            elseif isa( callback, 'function_handle' )
                callback( source, eventData )
            elseif iscell( callback )
                feval( callback{1}, source, eventData, callback{2:end} )
            end

        end % onSelectionChanged

        function onParentChanged( obj, ~, ~ )

            % Update ParentBackgroundColor and ParentBackgroundColor
            if isprop( obj.Parent, 'BackgroundColor' )
                prop = 'BackgroundColor';
            elseif isprop( obj.Parent, 'Color' )
                prop = 'Color';
            else
                prop = [];
            end

            if ~isempty( prop )
                obj.ParentBackgroundColorListener = event.proplistener( obj.Parent, ...
                    findprop( obj.Parent, prop ), 'PostSet', ...
                    @( src, evt ) obj.updateParentBackgroundColor( prop ) );
            else
                obj.ParentBackgroundColorListener = [];
            end

            obj.updateParentBackgroundColor( prop );

        end % onParentChanged

        function updateParentBackgroundColor( obj, prop )

            if isempty( prop )
                obj.ParentBackgroundColor = obj.BackgroundColor;
            else
                obj.ParentBackgroundColor = obj.Parent.(prop);
            end

            % Mark as dirty
            obj.Dirty = true;

        end

    end % event handlers

    methods( Access = private, Static )

        function mask = getDividerMask()
            %getDividerMask  Get divider image data
            %
            %  m = uix.TabPanel.getDividerMask() returns the image masks
            %  for tab panel dividers.  Mask entries are 0 (shadow), 1
            %  (background), 2 (tint) and 3 (highlight).

            mask.EF = indexColor( uix.loadIcon( 'tab_NoEdge_NotSelected.png' ) );
            mask.ET = indexColor( uix.loadIcon( 'tab_NoEdge_Selected.png' ) );
            mask.FE = indexColor( uix.loadIcon( 'tab_NotSelected_NoEdge.png' ) );
            mask.FF = indexColor( uix.loadIcon( 'tab_NotSelected_NotSelected.png' ) );
            mask.FT = indexColor( uix.loadIcon( 'tab_NotSelected_Selected.png' ) );
            mask.TE = indexColor( uix.loadIcon( 'tab_Selected_NoEdge.png' ) );
            mask.TF = indexColor( uix.loadIcon( 'tab_Selected_NotSelected.png' ) );

            function mask = indexColor( rgbMap )
                %indexColor  Returns a map of index given an RGB map
                %
                %  mask = indexColor( rgbMap ) returns a mask of color
                %  index based on the supplied rgbMap.
                %  black  : 0
                %  red    : 1
                %  yellow : 2
                %  white  : 3
                %  blue   : 4
                mask = nan( size( rgbMap, 1 ),size( rgbMap, 2 ) );
                % Black
                colorIndex = isColor( rgbMap, [0 0 0] );
                mask(colorIndex) = 0;
                % Red
                colorIndex = isColor( rgbMap, [1 0 0] );
                mask(colorIndex) = 1;
                % Yellow
                colorIndex = isColor( rgbMap, [1 1 0] );
                mask(colorIndex) = 2;
                % White
                colorIndex = isColor( rgbMap, [1 1 1] );
                mask(colorIndex) = 3;
                % Blue
                colorIndex = isColor( rgbMap, [0 0 1] );
                mask(colorIndex) = 4;
                % Nested
                function boolMap = isColor( map, color )
                    %isColor  Return a map of boolean where map is equal to color
                    boolMap = all( bsxfun( @eq, map, permute( color, [1 3 2] ) ), 3 );
                end
            end

        end % getDividerMask

    end % static helper methods

end % classdef

author: goGPS-Project