text:
%SerialLink.plot Graphical display and animation
%
% R.plot(Q, options) displays a graphical animation of a robot based on
% the kinematic model.  A stick figure polyline joins the origins of
% the link coordinate frames. The robot is displayed at the joint angle Q (1xN), or
% if a matrix (MxN) it is animated as the robot moves along the M-point trajectory.
%
% Options::
% 'workspace', W    Size of robot 3D workspace, W = [xmn, xmx ymn ymx zmn zmx]
% 'floorlevel',L    Z-coordinate of floor (default -1)
%-
% 'delay',D         Delay betwen frames for animation (s)
% 'fps',fps         Number of frames per second for display, inverse of 'delay' option
% '[no]loop'        Loop over the trajectory forever
% '[no]raise'       Autoraise the figure
% 'movie',M         Save an animation to the movie M
% 'trail',L         Draw a line recording the tip path, with line style L.
%                   L can be a cell array, eg. {'r', 'LineWidth', 2}
%-
% 'scale',S         Annotation scale factor
% 'zoom',Z          Reduce size of auto-computed workspace by Z, makes
%                   robot look bigger
% 'ortho'           Orthographic view
% 'perspective'     Perspective view (default)
% 'view',V          Specify view V='x', 'y', 'top' or [az el] for side elevations,
%                   plan view, or general view by azimuth and elevation
%                   angle.
% 'top'             View from the top.
%-
% '[no]shading'     Enable Gouraud shading (default true)
% 'lightpos',L      Position of the light source (default [0 0 20])
% '[no]name'        Display the robot's name
%-
% '[no]wrist'       Enable display of wrist coordinate frame
% 'xyz'             Wrist axis label is XYZ
% 'noa'             Wrist axis label is NOA
% '[no]arrow'       Display wrist frame with 3D arrows
%-
% '[no]tiles'       Enable tiled floor (default true)
% 'tilesize',S      Side length of square tiles on the floor
% 'tile1color',C   Color of even tiles [r g b] (default [0.5 1 0.5]  light green)
% 'tile2color',C   Color of odd tiles [r g b] (default [1 1 1] white)
%-
% '[no]shadow'      Enable display of shadow (default true)
% 'shadowcolor',C   Colorspec of shadow, [r g b]
% 'shadowwidth',W   Width of shadow line (default 6)
%-
% '[no]jaxes'       Enable display of joint axes (default false)
% '[no]jvec'        Enable display of joint axis vectors (default false)
% '[no]joints'      Enable display of joints
% 'jointcolor',C    Colorspec for joint cylinders (default [0.7 0 0])
% 'pjointcolor',C   Colorspec for prismatic joint boxes (default [0.4 1 .03])
% 'jointdiam',D     Diameter of joint cylinder in scale units (default 5)
%-
% 'linkcolor',C     Colorspec of links (default 'b')
%-
% '[no]base'        Enable display of base 'pedestal'
% 'basecolor',C     Color of base (default 'k')
% 'basewidth',W     Width of base (default 3)
%
% The options come from 3 sources and are processed in order:
% - Cell array of options returned by the function PLOTBOTOPT (if it exists)
% - Cell array of options given by the 'plotopt' option when creating the
%   SerialLink object.
% - List of arguments in the command line.
%
% Many boolean options can be enabled or disabled with the 'no' prefix.  The
% various option sources can toggle an option, the last value encountered is used.
%
% Graphical annotations and options::
%
% The robot is displayed as a basic stick figure robot with annotations
% such as:
% - shadow on the floor
% - XYZ wrist axes and labels
% - joint cylinders and axes
% which are controlled by options.
%
% The size of the annotations is determined using a simple heuristic from
% the workspace dimensions.  This dimension can be changed by setting the
% multiplicative scale factor using the 'mag' option.
%
% Figure behaviour::
%
% - If no figure exists one will be created and the robot drawn in it.
% - If no robot of this name is currently displayed then a robot will
%   be drawn in the current figure.  If hold is enabled (hold on) then the
%   robot will be added to the current figure.
% - If the robot already exists then that graphical model will be found
%   and moved.
%
% Multiple views of the same robot::
%
% If one or more plots of this robot already exist then these will all
% be moved according to the argument Q.  All robots in all windows with
% the same name will be moved.
%
%  Create a robot in figure 1
%         figure(1)
%         p560.plot(qz);
%  Create a robot in figure 2
%         figure(2)
%         p560.plot(qz);
%  Now move both robots
%         p560.plot(qn)
%
% Multiple robots in the same figure::
%
% Multiple robots can be displayed in the same plot, by using "hold on"
% before calls to robot.plot().
%
%  Create a robot in figure 1
%         figure(1)
%         p560.plot(qz);
%  Make a clone of the robot named bob
%         bob = SerialLink(p560, 'name', 'bob');
%  Draw bob in this figure
%         hold on
%         bob.plot(qn)
%
%  To animate both robots so they move together:
%         qtg = jtraj(qr, qz, 100);
%         for q=qtg'
%           p560.plot(q');
%           bob.plot(q');
%         end
%
% Making an animation::
%
% The 'movie' options saves the animation as a movie file or separate frames in a folder
% - 'movie','file.mp4' saves as an MP4 movie called file.mp4
% - 'movie','folder' saves as files NNNN.png into the specified folder
%   - The specified folder will be created
%   - NNNN are consecutive numbers: 0000, 0001, 0002 etc.
%   - To convert frames to a movie use a command like:
%          ffmpeg -r 10 -i %04d.png out.avi
%
% Notes::
% - The options are processed when the figure is first drawn, to make different options come
%   into effect it is neccessary to clear the figure.
% - The link segments do not neccessarily represent the links of the robot, they are a pipe
%   network that joins the origins of successive link coordinate frames.
% - Delay betwen frames can be eliminated by setting option 'delay', 0 or
%   'fps', Inf.
% - By default a quite detailed plot is generated, but turning off labels,
%   axes, shadows etc. will speed things up.
% - Each graphical robot object is tagged by the robot's name and has UserData
%   that holds graphical handles and the handle of the robot object.
% - The graphical state holds the last joint configuration
% - The size of the plot volume is determined by a heuristic for an all-revolute
%   robot.  If a prismatic joint is present the 'workspace' option is
%   required.  The 'zoom' option can reduce the size of this workspace.
%
% See also SerialLink.plot3d, plotbotopt, SerialLink.animate, SerialLink.teach.

% HANDLES:
%
% A robot comprises a bunch of individual graphical elements and these are
% kept in a structure:
%
%   h.link     the graphical elements that comprise each joint/link
%   h.wrist     the coordinate frame marking the wrist frame
%   h.shadow    the robot's shadow
%   h.floorlevel the z-coordinate of the tiled floor
%   h.robot pointer to the robot object
%   h.opt   the final options structure
%
% The h.link graphical element is tagged with the robot's name and has this
% struct as its UserData.
%
%  h.links -> h -> robot
%
%  This enables us to find all robots with a given name, in all figures,
% and update them.

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

% TODO
% deal with base transform and tool
% more consistent option names, scale, mag etc.
function plot(robot, qq, varargin)

    % check the joint angle data matches the robot
    n = robot.n;
    assert(numcols(qq) == n, 'RTB:SerialLink:plot:badarg', 'Insufficient columns in q')

    % process options, these come from:
    %  - passed arguments
    %  - the robot object itself
    %  - the file plotopt.m
    opt = RTBPlot.plot_options(robot, varargin);

    % logic to handle where the plot is drawn, are old figures updated or
    % replaced?
    %  calls create_floor() and create_robot() as required.

    if strcmp(get(gca,'Tag'), 'RTB.plot')
        % this axis is an RTB plot window

        rhandles = findobj('Tag', robot.name);

        if isempty(rhandles)
            % this robot doesnt exist here, create it or add it

            if ishold
                % hold is on, add the robot, don't change the floor
                handle = create_robot(robot, opt);

                % tag one of the graphical handles with the robot name and hang
                % the handle structure off it
%                 set(handle.joint(1), 'Tag', robot.name);
%                 set(handle.joint(1), 'UserData', handle);
            else
                % create the robot and floor
                newplot();

                RTBPlot.create_floor(opt);

                handle = create_robot(robot, opt);
                set(gca, 'Tag', 'RTB.plot');
            end

        end

    else
        % this axis never had a robot drawn in it before, let's use it
        RTBPlot.create_floor(opt);
        handle = create_robot(robot, opt);
        set(gca, 'Tag', 'RTB.plot');
        set(gcf, 'Units', 'Normalized');
        pf = get(gcf, 'Position');
%         if strcmp( get(gcf, 'WindowStyle'), 'docked') == 0
%             set(gcf, 'Position', [0.1 1-pf(4) pf(3) pf(4)]);
%         end
    end

    % deal with a few options that need to be stashed in the SerialLink object
    % movie mode has not already been flagged
    if opt.movie
        robot.movie = Animate(opt.movie);
    end
    robot.delay = opt.delay;
    robot.loop = opt.loop;

    if opt.raise
        % note this is a very time consuming operation
        figure(gcf);
    end

    if strcmp(opt.projection, 'perspective')
        set(gca, 'Projection', 'perspective');
    end

    if isstr(opt.view)
        switch opt.view
            case 'top'
                view(0, 90);
                set(gca, 'OuterPosition', [-0.8 -0.8 2.5 2.5]);
            case 'x'
                view(0, 0);
            case 'y'
                view(90, 0)
            otherwise
                error('rtb:plot:badarg', 'view must be: x, y, top')
        end
    elseif isnumeric(opt.view) && length(opt.view) == 2
        view(opt.view)
    end

    % enable mouse-based 3D rotation
    rotate3d on

    robot.trail = []; % clear the previous trail

    robot.animate(qq);

    if opt.movie
        robot.movie.close();
    end
end

% Create a new graphical robot in the current figure.
% Returns a structure of handles that describe the various graphical entities in the robot model
% Make extensive use of hgtransform, all entities are defined at the origin, then moved to their
% proper pose
%
% The graphical hiearchy is:
%  hggroup: Tag = robot name
%     hgtransform: Tag = 'link#'
%
% The top-level group has user data which is the handle structure.

function h = create_robot(robot, opt)

    %disp('creating new robot');

    links = robot.links;
    s = opt.scale;

    % create an axis
    ish = ishold();
    if ~ishold
        % if hold is off, set the axis dimensions
        axis(opt.workspace);
        hold on
    end

    N = robot.n;

    % create the base
    if opt.base
        bt = robot.base.t;
        bt = [bt'; bt'];
        bt(1,3) = opt.floorlevel;
        line(bt(:,1), bt(:,2), bt(:,3), 'LineWidth', opt.basewidth, 'Color', opt.basecolor);
    end
    view(3)

    % add the robot's name
    if opt.name
        b = robot.base.t;
        bz = 0;
        if opt.base
            bz = 0.5*opt.floorlevel;
        end
        text(b(1), b(2)-s, bz, [' ' robot.name], 'FontAngle', 'italic', 'FontWeight', 'bold')
    end
    group = hggroup('Tag', robot.name);
    h.group = group;

    % create the graphical joint and link elements
    for L=1:N
        if opt.debug
            fprintf('create graphics for joint %d\n', L);
        end

        % create the transform for displaying this element (joint cylinder + link)
        h.link(L) = hgtransform('Tag', sprintf('link%d', L), 'Parent', group);

        % create a joint cylinder
        if opt.joints
            % create the body of the joint
            if links(L).isrevolute
                RTBPlot.cyl('z', opt.jointdiam*s, opt.jointlen*s*[-1 1], opt.jointcolor, [], 'Parent', h.link(L));
            else
                % create an additional hgtransform for positioning and scaling the prismatic
                % element.  The element is created with unit length.
                h.pjoint(L) = hgtransform('Tag', 'prismatic', 'Parent', h.link(L));
                if links(L).mdh
                    % make the box extend in negative z-dir because scaling factor in animate
                    % must be positive
                    RTBPlot.box('z', s, [0 -1], opt.jointcolor, [], 'Parent', h.pjoint(L));
                else
                    RTBPlot.box('z', s, [0 1], opt.jointcolor, [], 'Parent', h.pjoint(L));
                end
            end
        end

        % create the body of the link

        % create elements to represent translation between joint frames, ie. the link itself.
        % This is drawn to resemble orthogonal plumbing.
        if robot.mdh
            % modified DH convention
            if L < N
                A = links(L+1).A(0);
                t = transl(A);
                if t(1) ~= 0
                    RTBPlot.cyl('x', s, [s t(1)], opt.linkcolor, [], 'Parent', h.link(L));
                end
                if t(2) ~= 0
                    RTBPlot.cyl('y', s, [s t(2)], opt.linkcolor, [t(1) 0 0], 'Parent', h.link(L));
                end
                if t(3) ~= 0
                    RTBPlot.cyl('z', s, [s t(3)], opt.linkcolor, [t(1) t(2) 0], 'Parent', h.link(L));
                end
            end
        else
            % standard DH convention
            if L > 1
                Ainv = inv(links(L-1).A(0));
                t = Ainv.t;
                if t(1) ~= 0
                    RTBPlot.cyl('x', s, [0 t(1)], opt.linkcolor, [], 'Parent', h.link(L));
                end
                if t(2) ~= 0
                    RTBPlot.cyl('y', s, [s t(2)], opt.linkcolor, [t(1) 0 0], 'Parent', h.link(L));
                end
                if t(3) ~= 0
                    RTBPlot.cyl('z', s, [s t(3)], opt.linkcolor, [t(1) t(2) 0], 'Parent', h.link(L));
                end
                %line([0 t(1)]', [0 t(2)]', [0 t(3)]', 'Parent', h.link(L));
            end
        end

        assert( ~(opt.jaxes && opt.jvec), 'RTB:plot:badopt', 'Can''t specify ''jaxes'' and ''jvec''')
        % create the joint axis line
        if opt.jaxes
            line('XData', [0 0], ...
                'YData', [0 0], ...
                'ZData', 14*s*[-1 1], ...
                'LineStyle', ':', 'Parent', h.link(L));

            % create the joint axis label
            text(0, 0, 14*s, sprintf('q%d', L), 'Parent', h.link(L))
        end
         % create the joint axis vector
        if opt.jvec
            daspect([1 1 1]);
            ha = arrow3([0 0 -12*s], [0 0 15*s], 'c');
            set(ha, 'Parent', h.link(L));

            % create the joint axis label
            text(0, 0, 20*s, sprintf(' q%d', L), 'Parent', h.link(L))
        end
    end
    if opt.debug
        fprintf('create graphics for tool\n');
    end
    % display the tool transform if it exists
    h.link(N+1) = hgtransform('Tag', sprintf('link%d', N+1), 'Parent', group);
    tool = SE3();
    if ~robot.mdh
        tool = links(L).A(0);
    end
    if ~isempty(robot.tool)
        tool = tool * robot.tool;
    end
    Tinv = inv(tool);
    t = Tinv.t;
    if t(1) ~= 0
        RTBPlot.cyl('x', s, [0 t(1)], 'r', [], 'Parent', h.link(N+1));
    end
    if t(2) ~= 0
        RTBPlot.cyl('y', s, [s t(2)], 'r', [t(1) 0 0], 'Parent', h.link(N+1));
    end
    if t(3) ~= 0
        RTBPlot.cyl('z', s, [s t(3)], 'r', [t(1) t(2) 0], 'Parent', h.link(N+1));
    end

    % display the wrist coordinate frame
    if opt.wrist
        if opt.arrow
            % compute arrow3 scale factor...
            ax = gca; d = [ax.XLim ax.YLim ax.ZLim];
            d = norm( d(4:6)-d(1:3) ) / 72;
            extra = {'arrow', 'width', 1.5*s/d};
        else
            extra = {};
        end
        h.wrist = trplot(eye(4,4), 'labels', upper(opt.wristlabel), ...
                'rgb', 'length', opt.wristlen*s, extra{:});

    else
        h.wrist = [];
    end

    % display a shadow on the floor
    if opt.shadow
        % create the polyline which is the shadow on the floor
        h.shadow = line('LineWidth', opt.shadowwidth, 'Color', opt.shadowcolor);
    end

    if ~isempty(opt.trail)
        if iscell(opt.trail)
            h.trail = plot(0, 0, opt.trail{:});
        else
            h.trail = plot(0, 0, opt.trail);
        end
        robot.trail = [];
    end

    % deal with some display options
    if opt.shading
        lighting gouraud
        light('position', opt.lightpos)
    end
    xlabel('X')
    ylabel('Y')
    zlabel('Z')
    grid on

    % restore hold setting
    if ~ish
        hold off
    end
    h.floorlevel = opt.floorlevel;
    h.robot = robot;
    h.opt = opt;

    % attach the handle structure to the top graphical element
    set(group, 'UserData', h);
end

% draw a cylinder of radius r in the direction specified by ax, with an
% extent from extent(1) to extent(2)

author: petercorke

repo: robotics-toolbox-matlab


text:
%% Wiener Deconvolution for Image Deblurring
%
% Sample shows how DFT can be used to perform
% <https://en.wikipedia.org/wiki/Wiener_deconvolution Weiner deconvolution>
% of an image with user-defined point spread function (PSF).
%
% Use controls to adjust PSF parameters, and swtich between linear/cirular PSF.
%
% See also: |deconvwnr|,
% <https://www.mathworks.com/help/images/image-restoration-deblurring.html>
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/deconvolution.py>
%

function varargout = weiner_deconvolution_demo_gui(im)
    % load grayscale image
    if nargin < 1
        im = fullfile(mexopencv.root(), 'test', 'licenseplate_motion.jpg');
        if exist(im, 'file') ~= 2
            disp('Downloading image...')
            url = 'https://cdn.rawgit.com/opencv/opencv/3.2.0/samples/data/licenseplate_motion.jpg';
            urlwrite(url, im);
        end
        img = cv.imread(im, 'Grayscale',true);
    elseif ischar(im)
        img = cv.imread(im, 'Grayscale',true);
    else
        img = im;
    end
    img = single(img) / 255;  % convert 8-bit to floating-point precision

    % create the UI
    app = initApp(img);
    h = buildGUI(app);
    if nargout > 0, varargout{1} = h; end
end

function app = initApp(img)
    %INIT_APP  Initialize options and state structure

    app = struct();

    % options
    app.snr = 25;       % signal/noise ratio in dB
    app.angle = 135;    % blur angle in degrees [0,180]
    app.d = 22;         % blur diameter [1,ksz]
    app.ksz = 65;       % kernel size

    % process image and compute its DFT
    app.img0 = img;
    app.img = blur_edge(img);
    app.IMG = cv.dft(app.img, 'ComplexOutput',true);
end

function img = blur_edge(img, d)
    %BLUR_EDGE  Blur image edges to reduce ringing effect in deblurred image
    %
    %     img = blur_edge(img)
    %     img = blur_edge(img, d)
    %
    % ## Input
    % * __img__ input image
    % * __d__ gaussian size, default 31
    %
    % ## Output
    % * __img__ output image
    %
    % See also: edgetaper
    %

    if nargin < 2, d = 31; end
    img_blur = cv.copyMakeBorder(img, [d d d d], 'BorderType','Wrap');
    img_blur = cv.GaussianBlur(img_blur, 'KSize',[d d]*2+1, 'SigmaX',-1);
    img_blur = img_blur(d+1:end-d, d+1:end-d);

    [h,w] = size(img);
    [y,x] = ndgrid(1:h, 1:w);
    D = min(cat(3, x, w-x+1, y, h-y+1), [], 3);
    a = min(single(D)/d, 1);

    img = img.*a + img_blur.*(1-a);
end

function kern = motion_kernel(ang, d, sz)
    %MOTION_KERNEL  Create linear motion filter
    %
    %     kern = motion_kernel(ang, d)
    %     kern = motion_kernel(ang, d, sz)
    %
    % ## Input
    % * __ang__ linear motion angle
    % * __d__ linear motion length
    % * __sz__ kernel size, default 65
    %
    % ## Output
    % * __kern__ kernel
    %
    % See also: fspecial (motion)
    %

    if nargin < 3, sz = 65; end
    sz2 = floor(sz / 2);
    A = [cos(ang) -sin(ang); sin(ang) cos(ang)];
    A(:,3) = [sz2; sz2] - A*[(d-1)*0.5; 0];
    kern = ones(1,d,'single');
    kern = cv.warpAffine(kern, A, 'DSize',[sz sz], 'Interpolation','Cubic');
end

function kern = defocus_kernel(d, sz)
    %DEFOCUS_KERNEL  Create circular defocus kernel
    %
    %     kern = defocus_kernel(d)
    %     kern = defocus_kernel(d, sz)
    %
    % ## Input
    % * __d__ circular motion diameter
    % * __sz__ kernel size, default 65
    %
    % ## Output
    % * __kern__ kernel
    %
    % See also: fspecial (gaussian)
    %

    if nargin < 2, sz = 65; end
    kern = zeros(sz,sz,'uint8');
    kern = cv.circle(kern, [sz sz], d, ...
        'Color',255, 'Thickness','Filled', 'LineType','AA', 'Shift',1);
    kern = single(kern) / 255;
end

function onChange(~,~,h,app)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    psf_type = get(h.pop, 'Value');
    psf_viz = get(h.cbox, 'Value') == get(h.cbox, 'Max');
    app.snr = round(get(h.slid(1), 'Value'));
    app.d = round(get(h.slid(2), 'Value'));
    app.angle = round(get(h.slid(3), 'Value'));
    set(h.txt(1), 'String',sprintf('SNR (dB): %d', app.snr));
    set(h.txt(2), 'String',sprintf('Diameter: %d', app.d));
    set(h.txt(3), 'String',sprintf('Angle: %d', app.angle));

    if psf_type == 3
        % show original blurred image
        res = app.img0;
    else
        % linear/cirular PSF
        if psf_type == 1
            psf = motion_kernel(deg2rad(app.angle), app.d, app.ksz);
        else
            psf = defocus_kernel(app.d, app.ksz);
        end

        % deconvolution
        noise = 10^(-0.1 * app.snr);
        [kh,kw] = size(psf);
        PSF = zeros(size(app.img), class(app.img));
        PSF(1:kh, 1:kw) = psf / sum(psf(:));
        PSF = cv.dft(PSF, 'ComplexOutput',true, 'NonzeroRows',kh);
        PSF = bsxfun(@rdivide, PSF, sum(PSF.^2, 3) + noise);
        res = cv.mulSpectrums(app.IMG, PSF);
        res = cv.dft(res, 'Inverse',true, 'Scale',true, 'RealOutput',true);
        res = circshift(res, floor(-[kh kw]/2));

        % visualize PSF kernel (overlaid on top of output image)
        if psf_viz
            res(1:kh,1:kw) = psf;
        end
    end

    % show result
    set(h.img, 'CData',res);
    drawnow;
end

function h = buildGUI(app)
    %BUILDGUI  Creates the UI

    % parameters
    sz = size(app.img0);
    sz(2) = max(sz(2), 300);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.fig = figure('Name','Deconvolution', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+105-1]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 105 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(app.img0, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(app.img0);
    end
    h.txt(1) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',11, 'HorizontalAlignment','left', ...
        'Position',[5 5 100 20], 'String','SNR (dB):');
    h.txt(2) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',11, 'HorizontalAlignment','left', ...
        'Position',[5 30 100 20], 'String','Diameter:');
    h.txt(3) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',11, 'HorizontalAlignment','left', ...
        'Position',[5 55 100 20], 'String','Angle:');
    h.txt(4) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',11, 'HorizontalAlignment','left', ...
        'Position',[5 80 100 20], 'String','PSF:');
    h.slid(1) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',app.snr, 'Min',0, 'Max',50, 'SliderStep',[1 5]./(50-0), ...
        'Position',[105 5 sz(2)-105-5 20]);
    h.slid(2) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',app.d, 'Min',1, 'Max',50, 'SliderStep',[1 5]./(50-1), ...
        'Position',[105 30 sz(2)-105-5 20]);
    h.slid(3) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',app.angle, 'Min',0, 'Max',180, 'SliderStep',[1 10]./(180-0), ...
        'Position',[105 55 sz(2)-105-5 20]);
    h.pop = uicontrol('Parent',h.fig, 'Style','popupmenu', ...
        'Position',[105 80 100 20], 'String',{'Linear','Circular','-None-'});
    h.cbox = uicontrol('Parent',h.fig, 'Style','checkbox', ...
        'Position',[210 80 100 20], 'Value',1, 'String','Show PSF');

    % hook event handlers, and trigger default start
    set([h.slid, h.pop, h.cbox], 'Callback',{@onChange,h,app}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h,app);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [bbox]=boundingBox(targetMask);
% find the bounding box of the binary structure targetMask
% bbox has form [xmin xmax ymin ymax zmin zmax]
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

  n=size(targetMask);

  %zmin
  cntr=1;
  true=sum(sum(targetMask(:,:,cntr)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(:,:,cntr)));
  end
  bbox(5)=cntr;

  %zmax
  cntr=0;
  true=sum(sum(targetMask(:,:,end-cntr)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(:,:,end-cntr)));
  end
  bbox(6)=n(3)-cntr;

  %ymin
  cntr=1;
  true=sum(sum(targetMask(:,cntr,:)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(:,cntr,:)));
  end
  bbox(3)=cntr;

  %ymax
  cntr=0;
  true=sum(sum(targetMask(:,end-cntr,:)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(:,end-cntr,:)));
  end
  bbox(4)=n(2)-cntr;

  %xmin
  cntr=1;
  true=sum(sum(targetMask(cntr,:,:)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(cntr,:,:)));
  end
  bbox(1)=cntr;

  %xmax
  cntr=0;
  true=sum(sum(targetMask(end-cntr,:,:)));

  while(~true)
    cntr=cntr+1;
    true=sum(sum(targetMask(end-cntr,:,:)));
  end
  bbox(2)=n(1)-cntr;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function B = barycentric_coordinates(P,varargin)
  % BARYCENTRIC_COORDINATES Computes barycentric coordinates of point p in
  % simplex (v1,v2,v3)
  %
  % B = barycentric_coordinates(P,V1,V2,V3, ...)
  %
  % Inputs:
  %   P  #P by dim list of query point locations
  %   V1  #P by dim list of simplex corner locations
  %   V2  #P by dim list of simplex corner locations
  %   ...
  % Outputs:
  %   B  #P by dim+1 list of barycentric coordinates
  %

  % SHOULD BE USING VARGIN TO ACCEPT ARBITRARY DIMENSION INPUT
  for varg = varargin
    assert(size(P,1) == size(varg{1},1), 'All inputs should be same length');
    assert(size(P,2) == size(varg{1},2), 'All inputs should be same dimension');
  end

  function v = volume(ad,bd,cd)
    r =[bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), ...
        bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), ...
        bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)];
    v = -sum(ad.*r,2)./6;
  end

  n = size(P,1);
  switch numel(varargin)
  case 4
    V1 = varargin{1};
    V2 = varargin{2};
    V3 = varargin{3};
    V4 = varargin{4};
    V1P = V1-P;
    V2P = V2-P;
    V3P = V3-P;
    V4P = V4-P;
    %T = bsxfun(@plus,(1:n)',(0:3)*n);
    A1 = volume(V2P,V4P,V3P);
    A2 = volume(V1P,V3P,V4P);
    A3 = volume(V1P,V4P,V2P);
    A4 = volume(V1P,V2P,V3P);
    A  = volume(V1-V4,V2-V4,V3-V4);
    if size(P,2)>3 && max(abs(sum([A1 A2 A3 A4],2)-A))>1e-14
      warning('Possibly negative coordinates. Not supported in dim~=3');
    end
    B = bsxfun(@rdivide,[A1 A2 A3 A4],A);
  case 3
    V1 = varargin{1};
    V2 = varargin{2};
    V3 = varargin{3};
    A1 = doublearea([ P;V2;V3],[1:n;n+[1:n;n+(1:n)]]');
    A2 = doublearea([V1; P;V3],[1:n;n+[1:n;n+(1:n)]]');
    A3 = doublearea([V1;V2; P],[1:n;n+[1:n;n+(1:n)]]');
    A  = doublearea([V1;V2;V3],[1:n;n+[1:n;n+(1:n)]]');
    if size(P,2)>2 && max(abs(sum([A1 A2 A3],2)-A))>1e-14
      warning('Possibly negative coordinates. Not supported in dim~=2');
    end
    B = bsxfun(@rdivide,[A1 A2 A3],A);
  case 2
    V1 = varargin{1};
    V2 = varargin{2};
    A1 = edge_lengths([P;V2],[1:n;n+[1:n]]');
    A2 = edge_lengths([V1;P],[1:n;n+[1:n]]');
    A = edge_lengths([V1;V2],[1:n;n+[1:n]]');
    if size(P,2)>1 && max(abs(sum([A1 A2],2)-A))>1e-14
      warning('Possibly negative coordinates. Not supported in dim~=1');
    end
    B = bsxfun(@rdivide,[A1 A2],A);
  otherwise
    error(sprintf('%d-simplices not supported',numel(varargin)));
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [lf] = halfspace_medium_leadfield(rd, elc, vol)

% HALFSPACE_MEDIUM_LEADFIELD calculate the halfspace medium leadfield
% on positions pnt for a dipole at position rd and conductivity cond
% The halfspace solution requires a plane dividing a conductive zone of
% conductivity cond, from a non coductive zone (cond = 0)
%
% [lf] = halfspace_medium_leadfield(rd, elc, cond)

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

siz = size(rd);
if any(siz==1)
  % positions are specified as a single vector
  Ndipoles = prod(siz)/3;
  rd = rd(:)'; % ensure that it is a row vector
elseif siz(2)==3
  % positions are specified as a Nx3 matrix -> reformat to a single vector
  Ndipoles = siz(1);
  rd = rd';
  rd = rd(:)'; % ensure that it is a row vector
else
  ft_error('incorrect specification of dipole locations');
end

Nelc     = size(elc,1);
lf       = zeros(Nelc,3*Ndipoles);

for i=1:Ndipoles
  % this is the position of dipole "i"
  dip1 = rd((1:3) + 3*(i-1));

  % distances electrodes - dipole
  r1 = elc - ones(Nelc,1) * dip1;

  % Method of mirror dipoles:
  % Defines the position of mirror dipoles being symmetric to the plane
  dip2 = get_mirror_pos(dip,vol);

  % distances electrodes - mirror dipole
  r2 = elc - ones(Nelc,1) * dip2;

  % denominator
  R1 =  (4*pi*vol.cond) * (sum(r1' .^2 ) .^ 1.5)';
  % denominator, mirror term
  R2 = -(4*pi*vol.cond) * (sum(r2' .^2 ) .^ 1.5)';

  % condition of dipoles falling in the non conductive halfspace
  condition = get_dip_halfspace(dip1,vol);

  invacuum = acos(dot(ori,(P-pnt)./norm(P-pnt))) < pi/2;

  if invacuum
    ft_warning('dipole lies on the vacuum side of the plane');
    lf(:,(1:3) + 3*(i-1)) = NaN(Nelc,3);
  elseif any(R1)==0
    ft_warning('dipole coincides with one of the electrodes');
    lf(:,(1:3) + 3*(i-1)) = NaN(Nelc,3);
  else
    lf(:,(1:3) + 3*(i-1)) = (r ./ [R1 R1 R1]) + (rp ./ [R2 R2 R2]);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [y,or1,or2,dmse] = predictor(x,dmodel)
%PREDICTOR  Predictor for y(x) using the given DACE model.
%
% Call:   y = predictor(x, dmodel)
%         [y, or] = predictor(x, dmodel)
%         [y, dy, mse] = predictor(x, dmodel)
%         [y, dy, mse, dmse] = predictor(x, dmodel)
%
% Input
% x      : trial design sites with n dimensions.
%          For mx trial sites x:
%          If mx = 1, then both a row and a column vector is accepted,
%          otherwise, x must be an mx*n matrix with the sites stored
%          rowwise.
% dmodel : Struct with DACE model; see DACEFIT
%
% Output
% y    : predicted response at x.
% or   : If mx = 1, then or = gradient vector/Jacobian matrix of predictor
%        otherwise, or is an vector with mx rows containing the estimated
%                   mean squared error of the predictor
% Three or four results are allowed only when mx = 1,
% dy   : Gradient of predictor; column vector with  n elements
% mse  : Estimated mean squared error of the predictor;
% dmse : Gradient vector/Jacobian matrix of mse

% hbn@imm.dtu.dk
% Last update August 26, 2002

    or1 = NaN; or2 = NaN; dmse = NaN;	% Default return values
    if isnan(dmodel.beta)
        error('DMODEL has not been found')
    end
    [m,n] = size(dmodel.S);     % number of design sites and number of dimensions
    sx    = size(x);            % number of trial sites and their dimension
    if min(sx) == 1 && n > 1    % Single trial point
        nx = max(sx);
        if nx == n
            mx = 1;
            x  = x(:).';
        end
    else
        mx = sx(1);
        nx = sx(2);
    end
    if nx ~= n
        error('Dimension of trial sites should be %d',n)
    end
    % Normalize trial sites
    x = (x - repmat(dmodel.Ssc(1,:),mx,1)) ./ repmat(dmodel.Ssc(2,:),mx,1);
    q = size(dmodel.Ysc,2);  % number of response functions
    if mx == 1  % one site only
        dx = repmat(x,m,1) - dmodel.S;  % distances to design sites
        if nargout > 1                  % gradient/Jacobian wanted
            [f,df] = feval(dmodel.regr, x);
            [r,dr] = feval(dmodel.corr, dmodel.theta, dx);
            % Scaled Jacobian
            dy = (df * dmodel.beta).' + dmodel.gamma * dr;
            % Unscaled Jacobian
            or1 = dy .* repmat(dmodel.Ysc(2, :)', 1, nx) ./ repmat(dmodel.Ssc(2,:), q, 1);
            if q == 1
                % Gradient as a column vector
                or1 = or1';
            end
            if nargout > 2  % MSE wanted
                rt = dmodel.C \ r;
                u = dmodel.Ft.' * rt - f.';
                v = dmodel.G \ u;
                or2 = repmat(dmodel.sigma2,mx,1) .* repmat((1 + sum(v.^2) - sum(rt.^2))',1,q);
                if nargout > 3  % gradient/Jacobian of MSE wanted
                    % Scaled gradient as a row vector
                    Gv = dmodel.G' \ v;
                    g = (dmodel.Ft * Gv - rt)' * (dmodel.C \ dr) - (df * Gv)';
                    % Unscaled Jacobian
                    dmse = repmat(2 * dmodel.sigma2',1,nx) .* repmat(g ./ dmodel.Ssc(2,:),q,1);
                    if q == 1
                    % Gradient as a column vector
                    dmse = dmse';
                    end
                end
            end
        else  % predictor only
            f = feval(dmodel.regr, x);
            r = feval(dmodel.corr, dmodel.theta, dx);
        end
        % Scaled predictor
        sy = f * dmodel.beta + (dmodel.gamma*r).';
        % Predictor
        y = (dmodel.Ysc(1,:) + dmodel.Ysc(2,:) .* sy)';
	else  % several trial sites
        % Get distances to design sites
        dx = zeros(mx*m,n);
        kk = 1 : m;
        for k = 1 : mx
            dx(kk,:) = repmat(x(k,:),m,1) - dmodel.S;
            kk = kk + m;
        end
        % Get regression function and correlation
        f = feval(dmodel.regr, x);
        r = feval(dmodel.corr, dmodel.theta, dx);
        r = reshape(r, m, mx);
        % Scaled predictor
        sy = f * dmodel.beta + (dmodel.gamma * r).';
        % Predictor
        y = repmat(dmodel.Ysc(1,:),mx,1) + repmat(dmodel.Ysc(2,:),mx,1) .* sy;
        if nargout > 1	% MSE wanted
            rt  = dmodel.C \ r;
            u   = dmodel.G \ (dmodel.Ft.' * rt - f.');
            or1 = repmat(dmodel.sigma2,mx,1) .* repmat((1 + sum(u.^2,1) - sum(rt.^2,1))',1,q);
            if  nargout > 2
                disp('WARNING from PREDICTOR.  Only  y  and  or1=mse  are computed')
            end
        end
    end
end

function [r,dr] = corrgauss(theta,d)
%CORRGAUSS  Gaussian correlation function,

    [m,n] = size(d);  % number of differences and dimension of data
    if length(theta) == 1
        theta = repmat(theta,1,n);
    elseif length(theta) ~= n
        error('Length of theta must be 1 or %d',n)
    end
    td = d.^2 .* repmat(-theta(:).',m,1);
    r  = exp(sum(td, 2));
	dr = repmat(-2*theta(:).',m,1) .* d .* repmat(r,1,n);
end

function [f,df] = regpoly0(S)
%REGPOLY0  Zero order polynomial regression function

    f  = ones(size(S,1),1);
	df = zeros(size(S,2),1);
end

function [f,df] = regpoly1(S)
%REGPOLY1  First order polynomial regression function

    f  = [ones(size(S,1),1),S];
	df = [zeros(size(S,2),1),eye(size(S,2))];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function view=xformFlatMenu(view)
%
% view = xformFlatMenu(view)
%
% Set up the callbacks for the xformView menu in the Flat view
%
% djh, 1/9/98
% rmk, 1/10/99 added xformParMap
% rmk, 1/15/99 added xformAllROIs

mrGlobals;

xformMenu = uimenu('Label','Xform','separator','on');

% Xform ROI callback:
%   gray=checkSelectedGray;
%   view=vol2flatCurROI(gray,view);
%   view=refreshScreen(view,0);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCurROI(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',0);'];
uimenu(xformMenu,'Label','ROI (selected)','Separator','off',...
    'CallBack',callBackstr);

% Xform AllROIs callback:
%   gray=checkSelectedGray;
%   view=vol2flatAllROIs(gray,view);
%   view=refreshScreen(view,0);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatAllROIs(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',0);'];
uimenu(xformMenu,'Label','ROIs (all)','Separator','off',...
    'CallBack',callBackstr);

% Xform CorAnal all scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnal(gray,view,0);
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnal(gray,',view.name,',0); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (all scans)','Separator','off',...
    'CallBack',callBackstr);

% Xform CorAnal current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnal(gray,view,getCurScan(view));
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnal(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (current scan)','Separator','off',...
    'CallBack',callBackstr);

% Xform CorAnal select scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatCorAnal(gray,view);
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatCorAnal(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','CorAnal (select scans)','Separator','off',...
    'CallBack',callBackstr);

% Xform ParMap all scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMap(gray,view,0);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMap(gray,',view.name,',0); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (all scans)','Separator','off',...
    'CallBack',callBackstr);

% Xform ParMap current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMap(gray,view,getCurScan(view));
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMap(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (current scan)','Separator','off',...
    'CallBack',callBackstr);

% Xform ParMap select scans callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMap(gray,view);
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMap(gray,',view.name,'); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (select scans)','Separator','off',...
    'CallBack',callBackstr);

% Xform ParMap current scan, nearest-neighbor callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMap2(gray,view,getCurScan(view));
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMap2(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (nearest-neighbor)','Separator','off',...
    'CallBack',callBackstr);

% Xform ParMap current scan callback:
%   gray=checkSelectedGray;
%   view=vol2flatParMapLevels(gray,view,getCurScan(view));
%   view=setDisplayMode(view,'map');
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=vol2flatParMapLevels(gray,',view.name,',getCurScan(',view.name,')); ',...
        view.name,'=setDisplayMode(',view.name,',''map''); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Parameter Map (separate gray levels)','Separator','off',...
    'CallBack',callBackstr);

% Xform retModel:
%   gray=checkSelectedGray;
%   view=rmVol2flat(gray,view);
%   view=refreshScreen(view,1);
callBackstr=['gray=checkSelectedGray; ',...
        view.name,'=rmVol2flat(gray,',view.name,'); ',...
        view.name,'=refreshScreen(',view.name,',1);'];
uimenu(xformMenu,'Label','Retinotopic Model','Separator','off',...
    'CallBack',callBackstr);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pd1 = FitNPV(NPV)
%CREATEFIT    Create plot of datasets and fits
%   PD1 = CREATEFIT(NPV)
%   Creates a plot, similar to the plot in the main distribution fitting
%   window, using the data that you provide as input.  You can
%   apply this function to the same data you used with dfittool
%   or with different data.  You may want to edit the function to
%   customize the code and this help message.
%
%   Number of datasets:  1
%   Number of fits:  1
%
%   See also FITDIST.

% This function was automatically generated on 12-Sep-2012 17:35:05

% Output fitted probablility distribution: PD1

% Data from dataset "NPV data":
%    Y = NPV

% Force all inputs to be column vectors
NPV = NPV(:);

% Prepare figure
clf;
hold on;
LegHandles = []; LegText = {};

% --- Plot data originally in dataset "NPV data"
[CdfF,CdfX] = ecdf(NPV,'Function','cdf');  % compute empirical cdf
BinInfo.rule = 1;
[~,BinEdge] = internal.stats.histbins(NPV,[],[],BinInfo,CdfF,CdfX);
[BinHeight,BinCenter] = ecdfhist(CdfF,CdfX,'edges',BinEdge);
hLine = bar(BinCenter,BinHeight,'hist');
set(hLine,'FaceColor','none','EdgeColor',[0.333333 0 0.666667],...
    'LineStyle','-', 'LineWidth',1);
xlabel('Data');
ylabel('Density')
LegHandles(end+1) = hLine;
LegText{end+1} = 'NPV data';

% Create grid where function will be computed
XLim = get(gca,'XLim');
XLim = XLim + [-1 1] * 0.01 * diff(XLim);
XGrid = linspace(XLim(1),XLim(2),100);

% --- Create fit "fit 1"

% Fit this distribution to get parameter values
% To use parameter estimates from the original fit:
%     pd1 = ProbDistUnivParam('normal',[ 3887789.069562, 1315866.846527])
pd1 = fitdist(NPV, 'normal');
YPlot = pdf(pd1,XGrid);
hLine = plot(XGrid,YPlot,'Color',[1 0 0],...
    'LineStyle','-', 'LineWidth',2,...
    'Marker','none', 'MarkerSize',6);
LegHandles(end+1) = hLine;
LegText{end+1} = 'fit 1';

% Adjust figure
box on;
hold off;

% Create legend from accumulated handles and labels
hLegend = legend(LegHandles,LegText,'Orientation', 'vertical', 'Location', 'NorthEast');
set(hLegend,'Interpreter','none');

author: Sable

repo: mcbench-benchmarks


text:
%% Bacterial foraging
% Animiation of bacteria movement to get the global minimum solution every chemotactic
%
% Author: Wael Mansour (wael192@yahoo.com)
%
% MSc Student, Electrical Enginering Dept,
% Faculty of Engineering Cairo University, Egypt

%%
%Initialization
clear all
clc
p=2;                         % dimension of search space
s=26;                        % The number of bacteria
Nc=50;                       % Number of chemotactic steps
Ns=4;                        % Limits the length of a swim
Nre=4;                       % The number of reproduction steps
Ned=2;                       % The number of elimination-dispersal events
Sr=s/2;                      % The number of bacteria reproductions (splits) per generation
Ped=0.25;                    % The probabilty that each bacteria will be eliminated/dispersed
c(:,1)=0.05*ones(s,1);       % the run length
for m=1:s                    % the initital posistions
    P(1,:,1,1,1)= 50*rand(s,1)';
    P(2,:,1,1,1)= .2*rand(s,1)';
   %P(3,:,1,1,1)= .2*rand(s,1)';
end

%%
%Main loop

%Elimination and dispersal loop
for ell=1:Ned

%Reprodution loop

    for K=1:Nre

%  swim/tumble(chemotaxis)loop

        for j=1:Nc

            for i=1:s
                J(i,j,K,ell)=Live_fn(P(:,i,j,K,ell));

% Tumble

                Jlast=J(i,j,K,ell);
                Delta(:,i)=(2*round(rand(p,1))-1).*rand(p,1);
                P(:,i,j+1,K,ell)=P(:,i,j,K,ell)+c(i,K)*Delta(:,i)/sqrt(Delta(:,i)'*Delta(:,i)); % This adds a unit vector in the random direction

% Swim (for bacteria that seem to be headed in the right direction)

                J(i,j+1,K,ell)=Live_fn(P(:,i,j+1,K,ell));
                m=0;         % Initialize counter for swim length
                    while m<Ns
                          m=m+1;
                          if J(i,j+1,K,ell)<Jlast
                             Jlast=J(i,j+1,K,ell);
                             P(:,i,j+1,K,ell)=P(:,i,j+1,K,ell)+c(i,K)*Delta(:,i)/sqrt(Delta(:,i)'*Delta(:,i)) ;
                             J(i,j+1,K,ell)=Live_fn(P(:,i,j+1,K,ell));
                          else
                             m=Ns ;
                          end

                    end
                J(i,j,K,ell)=Jlast;
                sprintf('The value of interation i %3.0f ,j = %3.0f  , K= %3.0f, ell= %3.0f' , i, j, K ,ell );

            end % Go to next bacterium

            x = P(1,:,j,K,ell);
            y = P(2,:,j,K,ell);
            clf
            plot(x, y , 'h')
            axis([-5 5 -5 5]);
            pause(.1)
        end  % Go to the next chemotactic

%Reprodution
        Jhealth=sum(J(:,:,K,ell),2);              % Set the health of each of the S bacteria
        [Jhealth,sortind]=sort(Jhealth);          % Sorts the nutrient concentration in order of ascending
        P(:,:,1,K+1,ell)=P(:,sortind,Nc+1,K,ell);
        c(:,K+1)=c(sortind,K);                    % And keeps the chemotaxis parameters with each bacterium at the next generation

%Split the bacteria (reproduction)
            for i=1:Sr
                P(:,i+Sr,1,K+1,ell)=P(:,i,1,K+1,ell); % The least fit do not reproduce, the most fit ones split into two identical copies
                c(i+Sr,K+1)=c(i,K+1);
            end
        end %  Go to next reproduction

%Eliminatoin and dispersal
        for m=1:s
            if  Ped>rand % % Generate random number
                P(1,:,1,1,1)= 50*rand(s,1)';
                P(2,:,1,1,1)= .2*rand(s,1)';
               %P(3,:,1,1,1)= .2*rand(s,1)';
            else
                P(:,m,1,1,ell+1)=P(:,m,1,Nre+1,ell); % Bacteria that are not dispersed
            end
        end
    end % Go to next elimination and disperstal

%Report
           reproduction = J(:,1:Nc,Nre,Ned);
           [jlastreproduction,O] = min(reproduction,[],2);  % min cost function for each bacterial
           [Y,I] = min(jlastreproduction)
           pbest=P(:,I,O(I,:),K,ell)


text:
function exact = p03_exact ( )

%*****************************************************************************80
%
%% P03_EXACT returns the exact integral for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = 2.0 / 3.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%NMC Nearest Mean Classifier
%
%   W = NMC(A)
%   W = A*NMC
%
% INPUT
%   A    Dataset
%
% OUTPUT
%   W    Nearest Mean Classifier
%
% DESCRIPTION
% Computation of the nearest mean classifier between the classes in the
% dataset A.  The use of soft labels is supported. Prior probabilities are
% not used.
%
% The difference with NMSC is that NMSC is based on an assumption of normal
% distributions and thereby automatically scales the features and is
% sensitive to class priors. NMC is a plain nearest mean classifier for
% which the assigned classes are are sensitive to feature scaling and
% unsensitive to class priors.
%
% The estimated class confidences by B*NMC(A), however, are based on
% assumed spherical gaussian distributions of the same size.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, NMSC, LDC, FISHERC, QDC, UDC

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: nmc.m,v 1.14 2009/12/09 15:49:54 duin Exp $

function W = nmc(varargin)

	mapname = 'NearestMean';
  argin = setdefaults(varargin,[],false);

  if mapping_task(argin,'definition')
    W = define_mapping(argin,'untrained',mapname);

  elseif mapping_task(argin,'training')			% Train a mapping.

    [a,flag] = deal(argin{:});
    % flag = true forces execution by knn_map instead of normal_map
    % flag = false (defaults) finds posteriors by spherical normal
    % distributions

    if nargin < 2, flag = false; end
    % flag=true forces to skip computation of posteriors based on normal
    % densities

    if nargin < 1 | isempty(a)
      W = prmapping(mfilename);
      W = setname(W,'NearestMean');
      return
    end

    islabtype(a,'crisp','soft');
    isvaldfile(a,1,2); % at least 1 object per class, 2 classes

    [m,k,c] = getsize(a);
    if isdatafile(a), a = setfeatsize(a,k); end

    if ~flag & c == 2      % 2-class case: store linear classifier
      u = meancov(a);
      u1 = +u(1,:);
      u2 = +u(2,:);
      R = [u1-u2]';
      offset =(u2*u2' - u1*u1')/2;
      W = affine([R -R],[offset -offset],a,getlablist(a));
      W = cnormc(W,a);
      W = setname(W,'Nearest Mean');
    else
      if ~flag & all(classsizes(a) > 1)
        a = setprior(a,0);  % NMC should be independent of priors: make them equal
        p = getprior(a);
        U = zeros(c,k);
        V = zeros(c,k);
        s = sprintf('NMC: Processing %i classes: ',c);
        prwaitbar(c,s);
        for j=1:c
          prwaitbar(c,j,[s int2str(j)]);
          b = seldat(a,j);
          [v,u] = var(b);
          U(j,:) = +u;
          V(j,:) = +v;
        end
        prwaitbar(0);
        %G = mean(V'*p') * eye(k);
        G = mean(V(:));
        w.mean = +U;
        w.cov = G;
        w.prior =p;
        W = normal_map(w,getlablist(a),k,c);
        W = setname(W,'Nearest Mean');
        W = setcost(W,a);
      else
        u = meancov(a);
        W = knnc(u,1);
        W = setname(W,mapname);
      end

    end

    W = setcost(W,a);

  end

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function RevertDirEdit(topH)

% dirS = RevertDirEdit(topH);
%
% Returns all the editable fields in the dialog to their original values,
% and reverts the directory structure.
%
% DBR 6/99

% Revert the session and directory structures:
uiData = get(topH, 'UserData');
dirS = uiData.origDirS;
uiData.dirS = dirS;
set(topH, 'UserData', uiData)

% Revert the directory dialog fields
nDir = length(uiData.dirData);
for iField=1:nDir
  editFlag = uiData.dirData(iField).edit;
  if editFlag
    evalStr = ['dirS.', uiData.dirData(iField).field];
    set(uiData.dirData(iField).handle, 'string', mat2str(eval(evalStr)));
  end
end


text:
function [Vout,HT_index, HT_values]=make_halfway_vertices(EV_table,ET_table,ETV_index,V,Ne)
% ETV_index_sparse=sparse(size(V,1),size(V,1));
%
% for i=1:size(ETV_index,1)
% 	if(ETV_index(i,1)>0)
% 		ETV_index_sparse(ETV_index(i,1),ETV_index(i,2))=i;
% 	else
% 		ET_table=ET_table(1:(i-1),:);
% 		EV_table=EV_table(1:(i-1));
% 		break
% 	end
% end

% Table to cell array
ETV_index_vall=cell(length(V),1);
for i=1:size(ETV_index,1)
    if(ETV_index(i,1)>0)
        ETV_index_vall{ETV_index(i,1)}=[ETV_index_vall{ETV_index(i,1)} i];
    else
        ET_table=ET_table(1:(i-1),:);
        EV_table=EV_table(1:(i-1));
        break;
    end
end

HT_index=cell(length(V),1);
HT_values=cell(length(V),1);

% Make output V
Vout=zeros(size(V,1)*4,3);
Vout(1:size(V,1),:)=V;
Vindex=size(V,1);

for i=1:length(V)
    Pneig=Ne{i};
    for j=1:length(Pneig);
        % Get the tangent and velocity of the edge P -> Pneig
        index=Ne{i}; vals=ETV_index_vall{i}; select1=vals(index==Pneig(j));
        Va=EV_table( select1); Ea=ET_table( select1,:);
        % Get the tangent and velocity of the edge Pneig -> P
        index=Ne{Pneig(j)}; vals=ETV_index_vall{Pneig(j)}; select2=vals(index==i);
        Vb=EV_table(select2); Eb=ET_table(select2,:);

        % The four points describing the spline
        P0=V(i,:);
        P3=V(Pneig(j),:);
        P1=P0+Ea*Va/3;
        P2=P3+Eb*Vb/3;

        % Spline used to calculated the xyz coordinate of the middle of each edge;
        c = 3*(P1 - P0);
        b = 3*(P2 - P1) - c;
        a = P3 - P0 - c - b;

        halfwayp = a*0.125 + b*0.250 + c*0.500 + P0;

        % Save the edge middle point
        if(sum(HT_index{i}==Pneig(j))==0)
            Vindex=Vindex+1;
            Vout(Vindex,:)=halfwayp;
            HT_index {i}= [HT_index{i} Pneig(j)];
            HT_values{i}=[HT_values{i} Vindex];
            HT_index {Pneig(j)}=[HT_index{ Pneig(j)} i];
            HT_values{Pneig(j)}=[HT_values{Pneig(j)} Vindex];
        end
    end
end
Vout=Vout(1:Vindex,:);


text:
function varargout = process_evt_simple( varargin )
% PROCESS_EVT_SIMPLE: Convert extended events to simple events
%
% USAGE:  OutputFiles = process_evt_simple('Run', sProcess, sInput)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2018

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Convert to simple event';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 63;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/EventMarkers#Other_menus';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Event name
    sProcess.options.eventname.Comment = 'Event names: ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = '';
    % Time offset
    sProcess.options.method.Comment = {'Keep the start of the events', 'Keep the middle of the events', 'Keep the end of the events'};
    sProcess.options.method.Type    = 'radio';
    sProcess.options.method.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFile = Run(sProcess, sInput) %#ok<DEFNU>
    % Return all the input files
    OutputFile = {sInput.FileName};

    % ===== GET OPTIONS =====
    % Time offset
    switch (sProcess.options.method.Value)
        case 1,  Method = 'start';
        case 2,  Method = 'middle';
        case 3,  Method = 'end';
    end
    % Event names
    EvtNames = strtrim(str_split(sProcess.options.eventname.Value, ',;'));
    if isempty(EvtNames)
        bst_report('Error', sProcess, [], 'No events selected.');
        return;
    end

    % ===== LOAD FILE =====
    % Get file descriptor
    isRaw = strcmpi(sInput.FileType, 'raw');
    % Load the raw file descriptor
    if isRaw
        DataMat = in_bst_data(sInput.FileName, 'F');
        sEvents = DataMat.F.events;
        sFreq = DataMat.F.prop.sfreq;
    else
        DataMat = in_bst_data(sInput.FileName, 'Events', 'Time');
        sEvents = DataMat.Events;
        sFreq = 1 ./ (DataMat.Time(2) - DataMat.Time(1));
    end
    % If no markers are present in this file
    if isempty(sEvents)
        bst_report('Error', sProcess, sInput, 'This file does not contain any event.');
        return;
    end
    % Find event names
    iEvtList = [];
    for i = 1:length(EvtNames)
        iEvt = find(strcmpi(EvtNames{i}, {sEvents.label}));
        if isempty(iEvt)
            bst_report('Warning', sProcess, sInput, 'This file does not contain any event.');
        elseif isempty(sEvents(iEvt).times)
            bst_report('Warning', sProcess, sInput, ['Event category "' sEvents(iEvt).label '" is empty.']);
        elseif (size(sEvents(iEvt).times,1) == 1)
            bst_report('Warning', sProcess, sInput, ['Event category "' sEvents(iEvt).label '" already contains simple events.']);
        else
            iEvtList(end+1) = iEvt;
        end
    end
    % No events to process
    if isempty(iEvtList)
        bst_report('Error', sProcess, sInput, 'No events to process.');
        return;
    end

    % ===== PROCESS EVENTS =====
    for i = 1:length(iEvtList)
        switch Method
            case 'start'
                sEvents(iEvtList(i)).times = sEvents(iEvtList(i)).times(1,:);
            case 'middle'
                sEvents(iEvtList(i)).times = mean(sEvents(iEvtList(i)).times,1);
            case 'end'
                sEvents(iEvtList(i)).times = sEvents(iEvtList(i)).times(2,:);
        end
        sEvents(iEvtList(i)).times = round(sEvents(iEvtList(i)).times .* sFreq) / sFreq;
    end

    % ===== SAVE RESULT =====
    % Report results
    if isRaw
        DataMat.F.events = sEvents;
    else
        DataMat.Events = sEvents;
    end
    % Only save changes if something was change
    bst_save(file_fullpath(sInput.FileName), DataMat, 'v6', 1);
end

author: brainstorm-tools


text:
function varargout = centerOfMass(A,varargin)
% CENTEROFMASS finds the center of mass of the N-dimensional input array
%
%   CENTEROFMASS(A) finds the gray-level-weighted center of mass of the
%   N-dimensional numerical array A. A must be real and finite. A warning
%   is issued if A contains any negative values. Any NaN elements of A will
%   automatically be ignored. CENTEROFMASS produces center of mass
%   coordinates in units of pixels. An empty array is returned if the
%   center of mass is undefined.
%
%   The center of mass is reported under the assumption that the first
%   pixel in each array dimension is centered at 1.
%
%   Also note that numerical arrays other than DOUBLE and SINGLE are
%   converted to SINGLE in order to prevent numerical roundoff error.
%
%   Examples:
%       A = rgb2gray(imread('saturn.png'));
%       C = centerOfMass(A);
%
%       figure; imagesc(A); colormap gray; axis image
%       hold on; plot(C(2),C(1),'rx')
%
%   See also:
%
%

%
%   Jered R Wells
%   2013/05/07
%   jered [dot] wells [at] gmail [dot] com
%
%   v1.0
%
%   UPDATES
%       YYYY/MM/DD - jrw - v1.1
%
%

%% INPUT CHECK
narginchk(0,1);
nargoutchk(0,1);
fname = 'centerOfMass';

% Checked required inputs
validateattributes(A,{'numeric'},{'real','finite'},fname,'A',1);

%% INITIALIZE VARIABLES
A(isnan(A)) = 0;
if ~(strcmpi(class(A),'double') || strcmpi(class(A),'single'))
    A = single(A);
end
if any(A(:)<0)
    warning('MATLAB:centerOfMass:neg','Array A contains negative values.');
end

%% PROCESS
sz = size(A);
nd = ndims(A);
M = sum(A(:));
C = zeros(1,nd);
if M==0
    C = [];
else
    for ii = 1:nd
        shp = ones(1,nd);
        shp(ii) = sz(ii);
        rep = sz;
        rep(ii) = 1;
        ind = repmat(reshape(1:sz(ii),shp),rep);
        C(ii) = sum(ind(:).*A(:))./M;
    end
end

% Assemble the VARARGOUT cell array
varargout = {C};

end % MAIN

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function c8vec_uniform_01_test ( )

%*****************************************************************************80
%
%% C8VEC_UNIFORM_01_TEST tests C8VEC_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8VEC_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C8VEC_UNIFORM_01 computes pseudorandom complex values\n' );
  fprintf ( 1, '  in the unit circle.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );

  n = 10;

  [ x, seed ] = c8vec_uniform_01 ( n, seed );

  for i = 1 : n
    fprintf ( 1, '  %6d  ( %f, %f )\n', i, real ( x(i) ), imag ( x(i) ) );
  end

  return
end


text:
function obj = addHolonomicConstraint(obj, constr, load_path)
    % Adds a holonomic constraint to the dynamical system
    %
    %
    % Parameters:
    % constr: the expression of the constraints @type HolonomicConstraint
    % load_path: the path from which the symbolic expressions for the
    % input map can be loaded @type char
    %
    % @note Any holonomic constraint introduces a set of constrained
    % wrenchs to the system (Lagrangian multiplier).

    % validate input argument
    validateattributes(constr, {'HolonomicConstraint'},...
        {},'ContinuousDynamics', 'HolonomicConstraint');

    n_constr = numel(constr);

    for i=1:n_constr
        c_obj = constr(i);
        c_name = c_obj.Name;

        if isfield(obj.HolonomicConstraints, c_name)
            error('The holonomic constraint (%s) has been already defined.\n',c_name);
        else

            % add virtual constraint
            obj.HolonomicConstraints.(c_name) = c_obj;

            % add constant parameters
            obj = addParam(obj, c_obj.ParamName, c_obj.Param);
            Jh = c_obj.ConstrJac;
            if ~exist('load_path','var')
                obj = addInput(obj, 'ConstraintWrench', c_obj.InputName, c_obj.Input, transpose(Jh));
            else
                obj = addInput(obj, 'ConstraintWrench', c_obj.InputName, c_obj.Input, transpose(Jh),'LoadPath',load_path);
            end
        end

    end

end

author: ayonga

repo: frost-dev

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


text:
function [x, funVal, ValueL]=glLeastR(A, y, z, opts)
%
%%
% Function glLeastR
%      Least Squares Loss with the (group) L1/Lq-norm Regularization
%
%% Problem
%
%  min  1/2 || A x - y||^2 + z * sum_j gWeight_j ||x^j||_q
%
%  x is grouped into k groups according to opts.ind
%  The indices of x_j in x is (ind(j)+1):ind(j+1)
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        L1/Lq norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%
%  x-         Solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 19, 2010.
%
%% Related papers
%
% [1]  Jun Liu, Shuiwang Ji, and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, UAI, 2009
%
% [2]  Jun Liu, Lei Yuan, Songcan Chen and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, Technical Report ASU, 2009.
%
%% Related functions:
%
%  sll_opts, initFactor, pathSolutionLeast
%  glLogisticR, eppVectorR
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts.ind should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z<=0)
    error('\n z should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind and q
if (~isfield(opts,'ind'))
    error('\n In glLeastR, the field .ind should be specified');
else
    ind=opts.ind;

    k=length(ind)-1; % the number of groups

    if (ind(k+1)~=n)
        error('\n Check opts.ind');
    end
end

if (~isfield(opts,'q'))
    q=2; opts.q=2;
else
    q=opts.q;
    if (q<1)
        error('\n q should be larger than 1');
    end
end

% gWeight: the weigtht for each group
if (isfield(opts,'gWeight'))
    gWeight=opts.gWeight;
    if (size(gWeight,1)~=k)
        error('\n opts.gWeight should a %d x 1 vector',k);
    end

    if (min(gWeight)<=0)
        error('\n .gWeight should be positive');
    end
else
    gWeight=ones(k,1);
end

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy =A'*y;
elseif (opts.nFlag==1)
    ATy= A'*y - sum(y) * mu';  ATy=ATy./nu;
else
    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
end

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
    if (z<0 || z>1)
        error('\n opts.rFlag=1, and z should be in [0,1]');
    end

    if q==1
        q_bar=Inf;
    elseif q>=1e6
        q_bar=1;
    else
        q_bar=q/(q-1);
    end

    % compute the norm of ATy corresponding to each group
    norm_ATy=zeros(k,1);
    for i=1:k
        norm_ATy(i,1)=norm(  ATy( (1+ind(i)):ind(i+1) ), q_bar );
    end

    % incorporate the gWeight
    norm_ATy=norm_ATy./gWeight;

    % compute lambda_max
    lambda_max=max(norm_ATy);

    % As .rFlag=1, we set lambda as a ratio of lambda_max
    lambda=z*lambda_max;
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
end

if (opts.init==0) % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail

    % the q-norm of x
    x_norm=zeros(k,1);
    for i=1:k
        x_norm(i,1)=norm(  x( (1+ind(i)):ind(i+1) ), q);
    end

    sum_x_norm=x_norm'*gWeight;

    if x_norm>=1e-10
        ratio=initFactor(sum_x_norm, Ax, y, lambda,'glLeastR');
        x=ratio*x;    Ax=ratio*Ax;
    end
end

%% The main program
% The Armijo Goldstein line search schemes + accelearted gradient descent

bFlag=0; % this flag tests whether the gradient step only changes a little

if (opts.mFlag==0 && opts.lFlag==0)
    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;

            % L1/Lq-norm regularized projection
            if (q<1e6)
                x=eppVector(v, ind, k, n, lambda/ L * gWeight, q);
            else % when q>=1e6, we treat q as inf
                x=eppVector(v, ind, k, n, lambda/ L * gWeight, 1e6);
            end

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            Av=Ax -As;
            r_sum=v'*v; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;

        ValueL(iterStep)=L;

        % the q-norm of x
        x_norm=zeros(k,1);
        for i=1:k
            x_norm(i,1)=norm(  x( (1+ind(i)):ind(i+1) ), q);
        end

        % function value = loss + regularizatioin
        funVal(iterStep)=Axy'* Axy/2 + lambda * x_norm'* gWeight;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% Reformulated problem + Nemirovski's line search
%
%  Original problem:
%  min  1/2 || A x - y||^2 + z * sum_j ||x^j||_q
%
%  Reformulated problem
%  min  1/2 || A x - y||^2  + z * gWeight' * t
%  ||x^j||_q <=t_j
%
%  We only deal with q=2

if (opts.mFlag==1 && opts.lFlag==0 && opts.q==2)

    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);

    t=zeros(k,1);
    for i=1:k
        t(i,1)=norm(  x( (1+ind(i)):ind(i+1) ), 2);
    end
    tp=t;

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp; s_t=t + beta * (t-tp);

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        tp=t;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            u=s-g/L; v= s_t - lambda * gWeight / L;

            % projection
            [x, t] = eppVectorR(u, v, ind, n, k);

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            Av=Ax -As;
            r_sum=v'*v + norm(t-s_t)^2; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;

        xxp=x-xp;   Axy=Ax-y;

        % function value = loss + regularizatioin
        funVal(iterStep)=Axy'* Axy/2 + lambda * t'* gWeight;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% adaptive line search

% .mFlag=1, and .lFlag=1
%  refomulate the problem as the constrained convex optimization
%  problem, and then apply adaptive line search scheme

% Problem:
%    min  1/2 || A x - y||^2 + z * t' * gWeight
%    s.t.   |x| <= t

if(opts.mFlag==1 && opts.lFlag==1 && opts.q==2)

    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    gamma=1;
    % we shall set the value of gamma = L,
    % where L is appropriate for the starting point

    xp=x; Axp=Ax;
    % store x and Ax
    xxp=zeros(n,1);
    % the difference of x and xp

    t=zeros(k,1);
    for i=1:k
        t(i,1)=norm(  x( (1+ind(i)):ind(i+1) ), 2);
    end
    tp=t;
    % t is the upper bound of the 2-norm of x

    % compute AT Ax
    if (opts.nFlag==0)
        ATAx=A'*Ax;
    elseif (opts.nFlag==1)
        ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
    else
        invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
    end

    % We begin the adaptive line search in the following
    %
    % Note that, in the line search, L and beta are changing

    for iterStep=1:opts.maxIter

        ATAxp=ATAx;
        % store ATAx to ATAxp

        if (iterStep~=1)
            % compute AT Ax
            if (opts.nFlag==0)
                ATAx=A'*Ax;
            elseif (opts.nFlag==1)
                ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
            else
                invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
            end
        end

        %--------- Line Search for L begins
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;   s_t= t + beta * (t -tp);
                As=Ax + beta* (Ax-Axp);
                ATAs=ATAx + beta * (ATAx- ATAxp);
                % compute the search point s, A * s, and A' * A * s
            else
                alpha= (-1+ sqrt(5)) / 2;
                beta=0; s=x; s_t=t; As=Ax; ATAs=ATAx;
            end

            g=ATAs-ATy;
            % compute the gradient g

            % let s walk in a step in the antigradient of s
            u=s-g/L; v= s_t - lambda * gWeight / L;

            % projection
            [xnew, tnew] = eppVectorR(u, v, ind, n, k);

            v=xnew-s;  % the difference between the new approximate solution x
                            % and the search point s
            v_t=tnew-s_t;

            % compute A xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            Av=Axnew -As;
            r_sum=v'*v + v_t'*v_t; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2
            %                       <= L * (||v||_2^2 + ||v_t|| _2^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end
        %--------- Line Search for L ends

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;    x=xnew; xxp=x-xp;
        Axp=Ax;  Ax=Axnew;
        % update x and Ax with xnew and Axnew
        tp=t; t=tnew;
        % update tp and t

        Axy=Ax-y;
        funVal(iterStep)=Axy'* Axy/2 + lambda * t'* gWeight;
        % compute function value

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp+ norm(t-tp)^2);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp + tp'*tp);    norm_xxp=sqrt(xxp'*xxp+ norm(t-tp)^2);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%%
if(opts.mFlag==0 && opts.lFlag==1)
    error('\n The function does not support opts.mFlag=0 & opts.lFlag=1!');
end

author: thomaskuestner


text:
function [I J K]=sphere_index(r,Ic,Jc,Kc)

% function [I J K]=sphere_index(r)
% ------------------------------------------------------------------------
% This function generates the indices 'I',  'J' and 'K' of voxel centers
% found within a sphere of radius 'r' (in voxels). The indices are not
% 'real' indices but are 'centered around zero'. They can be used to
% generate the indices of voxels found inside a sphere around point 'n' by
% adding In,Jn,Kn to I,J and K.
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 07/08/2008
% ------------------------------------------------------------------------

if nargin == 1
    Ic=0; Jc=0; Kc=0;
end

%Set-up mesh calculate radius
[X,Y,Z] = meshgrid((-round(r+1)):1:(round(r+1)));
X=X+Jc; Y=Y+Ic; Z=Z+Kc;
radius = hypot(hypot(X,Y),Z);
[I,J,K]=ind2sub(size(X),find(radius<=r));
IJK_middle=round(size(X)/2);
I=I-IJK_middle(1);
J=J-IJK_middle(2);
K=K-IJK_middle(3);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function cnt = countcover(sz,blocksize,stepsize)
%COUNTCOVER Covering of signal samples by blocks
%  CNT = COUNTCOVER(SZ,BLOCKSIZE,STEPSIZE) assumes a p-dimensional signal
%  of size SZ=[N1 N2 ... Np] covered by (possibly overlapping) blocks of
%  size BLOCKSIZE=[M1 M2 ... Mp]. The blocks start at position (1,1,..,1)
%  and are shifted between them by steps of size STEPSIZE=[S1 S2 ... Sp].
%  COUNTCOVER returns a matrix the same size as the signal, containing in
%  each entry the number of blocks covering that sample.
%
%  See also IM2COLSTEP, COL2IMSTEP, IM2COL.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2008

cnt = ones(sz);
for k = 1:length(sz)

  % this code is modified from function NDGRID, so it computes one
  % output argument of NDGRID at a time (to conserve memory)
  ids = (1:sz(k))';
  s = sz; s(k) = [];
  ids = reshape(ids(:,ones(1,prod(s))),[length(ids) s]);
  ids = permute(ids,[2:k 1 k+1:length(sz)]);

  cnt = cnt .* max( min(floor((ids-1)/stepsize(k)),floor((sz(k)-blocksize(k))/stepsize(k))) - ...
                    max(ceil((ids-blocksize(k))/stepsize(k)),0) + 1 , 0 );
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function cinh_values_test ( )

%*****************************************************************************80
%
%% CINH_VALUES_TEST demonstrates the use of CINH_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CINH_VALUES_TEST:\n' );
  fprintf ( 1, '  CINH_VALUES stores values of\n' );
  fprintf ( 1, '  the Hyperbolic Cosine Integral function CINH(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            CINH(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

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

  end

  return
end


text:
close all;
clear all;
clc;

data_file_path = 'bin/ols_phase_transition_gaussian_dict_gaussian_data.mat';
options.export = false;
options.export_dir = 'bin';
options.export_name = 'gaussian_dict_gaussian_data';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'OLS', options);

author: indigits

repo: sparse-plex


text:
% Convert edge list to adjacency matrix.
%
% INPUTS: edge list: mx3, m - number of edges
% OUTPUTS: adjacency matrix nxn, n - number of nodes
%
% Note: information about nodes is lost: indices only (i1,...in) remain
% GB: last updated, Sep 25, 2012

function adj=edgeL2adj(el)

nodes=sort(unique([el(:,1) el(:,2)])); % get all nodes, sorted
adj=zeros(numel(nodes));         % initialize adjacency matrix

% across all edges
for i=1:size(el,1); adj(find(nodes==el(i,1)),find(nodes==el(i,2)))=el(i,3); end

author: aeolianine

repo: octave-networks-toolbox


text:
function [planC] = dicomrt_dicomrt2cerr(plan,xmesh,ymesh,zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,voi,optName)
% dicomrt_dicomrt2CERR(dose,xmesh,ymesh,zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,voi)
%
% Convert DICOMRT-TOOLBOX data in CERR format
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(9,10,nargin))

if exist('optName')==0
    optName='CERROptions.m';
end

%Create CERR data Container & templates for certain cells.
planInitC       = initializeCERR;
indexS          = planInitC{end};
headerInitS     = planInitC{indexS.header};
commentInitS    = planInitC{indexS.comment};
scanInitS       = initializeScanInfo;
structureInitS  = planInitC{indexS.structures};
beamsInitS      = planInitC{indexS.beams};
doseInitS       = planInitC{indexS.dose};

tagMapS = initTagMapS;
tags = dicomrt_d2c_rtogtags;

% Load options
optS=opts4Exe(optName);

% Initialize planC
planC = planInitC;
nCells = length(planInitC);
tags.nimages = 0;

% Convert header
disp('(+) Converting DICOM headers');
if isempty(plan)==0
    headerInitS=dicomrt_d2c_header(headerInitS,plan);
    disp('(+) Fetching data from RTPLAN');

    % Convert beam geometry
    disp('(+) Converting DICOM beam geometry');
    [tmp_beam,tags]=dicomrt_d2c_beam(beamsInitS,indexS,plan,tags);
    disp('(=) Conversion DICOM beam geometry completed.');

    % Convert dose
    disp('(+) Converting DICOM dose');
    [tmp_dose,tags]=dicomrt_d2c_dose(doseInitS,indexS,plan,xmesh,ymesh,zmesh,tags);
    disp('(=) Conversion DICOM dose completed.');

    % Store header info
    planC{indexS.header} = headerInitS;

    % Store beams
    planC{indexS.beams}=tmp_beam;

    % Store dose
    planC{indexS.dose}=tmp_dose;
elseif isempty(plan)==1 & isempty(ct)==0
    headerInitS=dicomrt_d2c_header(headerInitS,ct);
    disp('(+) Fetching data from CT scan');
else
    headerInitS=dicomrt_d2c_header(headerInitS,voi);
    disp('(+) Fetching data from RTSTRUCT');
end
disp('(=) Conversion DICOM headers completed.');

% Convert comments
% commentInitS

% Convert scans
if isempty(ct)==0
    disp('(+) Converting DICOM scans');
    [tmp_scan,tags]=dicomrt_d2c_scan(scanInitS,indexS,ct,ct_xmesh,ct_ymesh,ct_zmesh,tags);
    disp('(=) Conversion DICOM scans completed.');

    % Store scan
    planC{indexS.scan}=tmp_scan;
else
    disp('(=) CT scan not available.');
end

% Validating and converting structures
% Validation
if isempty(voi)==0
    voitype=dicomrt_checkvoitype(voi);
    if isequal(voitype,'3D')
        voi=dicomrt_3dto2dVOI(voi);
    end

    if isempty(ct)==0
        disp('(+) Fitting DICOM structures to CT scans to comply with RTOG format');
        voi=dicomrt_fitvoi2ct(voi,ct_zmesh,ct);
        disp('(=) Fitting DICOM structures completed.');

        disp('(+) Validating DICOM structures');
        [voi]=dicomrt_validatevoi(voi,ct_xmesh,ct_ymesh,ct_zmesh);
        [newvoi]=dicomrt_closevoi(voi);
        disp('(=) Validation DICOM structures completed.');

        % Conversion
        disp('(+) Converting DICOM structures');
        [tmp_str,tags]=dicomrt_d2c_voi(structureInitS,indexS,newvoi,ct_xmesh,ct_ymesh,ct_zmesh,tags);
        disp('(=) Conversion DICOM structures completed.');

    else
        % Conversion
        disp('(+) Converting DICOM structures');
        [tmp_str,tags]=dicomrt_d2c_voi(structureInitS,indexS,voi,ct_xmesh,ct_ymesh,ct_zmesh,tags);
        disp('(=) Conversion DICOM structures completed.');
    end

    % Store structures
    planC{indexS.structures}=tmp_str;
else
    disp('(=) RTSTRUCT not available.');
end

% Convert digital films
nfilms=0;

% Convert DVH
% N/A

% Convert seed geometry

% Store options and indexS
planC{indexS.CERROptions} = optS;
planC{indexS.indexS} = indexS;

% Post processing actions

if isempty(ct)==0
    % Determine voxel thicknesses
    planC=dicomrt_d2c_setVoxelThicknesses(planC,indexS);
end

if isempty(ct)==0 & isempty(voi)==0
    % Get structure scan segments
    planC =  getRasterSegs(planC, optS);

    % Get any dose surface points
    planC =  getDSHPoints(planC, indexS, optS);
end

if isempty(ct)==0
    % Create uniformized datasets
    planC=dicomrt_d2c_uniformizescan(planC,optS);
end

% Save plan
save_planC(planC,optS);


text:
function x = emailFeatures(word_indices)
%EMAILFEATURES takes in a word_indices vector and produces a feature vector
%from the word indices
%   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
%   produces a feature vector from the word indices.

% Total number of words in the dictionary
n = 1899;

% You need to return the following variables correctly.
x = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return a feature vector for the
%               given email (word_indices). To help make it easier to
%               process the emails, we have have already pre-processed each
%               email and converted each word in the email into an index in
%               a fixed dictionary (of 1899 words). The variable
%               word_indices contains the list of indices of the words
%               which occur in one email.
%
%               Concretely, if an email has the text:
%
%                  The quick brown fox jumped over the lazy dog.
%
%               Then, the word_indices vector for this text might look
%               like:
%
%                   60  100   33   44   10     53  60  58   5
%
%               where, we have mapped each word onto a number, for example:
%
%                   the   -- 60
%                   quick -- 100
%                   ...
%
%              (note: the above numbers are just an example and are not the
%               actual mappings).
%
%              Your task is take one such word_indices vector and construct
%              a binary feature vector that indicates whether a particular
%              word occurs in the email. That is, x(i) = 1 when word i
%              is present in the email. Concretely, if the word 'the' (say,
%              index 60) appears in the email, then x(60) = 1. The feature
%              vector should look like:
%
%              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
%
%
for i=1:size(word_indices)
    x(word_indices(i)) = 1;
end
% =========================================================================

end

author: UtkarshPathrabe


text:
function [ x, seed ] = r8vec_uniform_ab ( n, a, b, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_AB returns a scaled pseudorandom R8VEC.
%
%  Discussion:
%
%    Each dimension ranges from A to B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A, B, the range of the pseudorandom values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N,1), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_AB - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_AB - Fatal error!' );
  end

  x = zeros ( n, 1 );

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    x(i) = a + ( b - a ) * seed * 4.656612875E-10;

  end

  return
end