text:
function varargout = mri_editMask(varargin)
% MRI_EDITMASK: Edit a 3D binary mask for a MRI.
%
% USAGE:  finalMask = mri_editMask(Mri)
%         finalMask = mri_editMask(Mri, Voxsize)
%         finalMask = mri_editMask(Mri, Voxsize, InitialMask)
%         finalMask = mri_editMask(Mri, Voxsize, InitialMask, InitialPos)
%         finalMask = mri_editMask(Mri, Voxsize, InitialMask, InitialPos, ColormapName)
%
% INPUT:
%    - Mri         : Full 3D matrix
%    - Voxsize     : Voxel size of the Mri volume (x,y,z)
%    - InitialMask : Full 3D matrix representing the initial mask
%                    (optional, if not specified, use an empty volume)
%    - InitialPos  : (dim, sliceIndex), slice to display when opening the window
%    - ColormapName: colormap name to display the anatomy
%
% OUTPUT:
%    - finalMask = final binary 3D mask
%
% A global variable is used to store application MRI and Masks : gEditMaskData
% Structure with following fields :
%    |- mri           : MRI that is currently processed
%    |- voxsize       : size of each mri voxel (x,y,z)
%    |- mask          : currently edited mask (for brain or grey matter)
%    |- isClicked     : structure that handles the mouse clicks and movements
%    |     |- source              : name of the object clicked by the user
%    |     |- clickPositionAxes   : click position in figure coordinates
%    |     |- clickPositionFigure : click position in axis coordinates
%    |- sliderValue   : [1,3] with the index of the current slice in each
%    |                  orientation for preview image (x,y,z)
%    |- orientation   : current orientation of preview slice (x,y,z)
%    |- intensityMax  : maximum intensity in image
%    |- colormap
%    |- colormapMax   : value above which all intensity are displayed as white

% @=============================================================================
% 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, 2006-2010

%% ===== CHECK INPUT =====
global gEditMaskData;
% === EXTERNAL CALL ===
if (nargin >= 1) && isnumeric(varargin{1})
    % === PARSE INPUTS ===
    % Argument #1: MRI
    gEditMaskData.mri = varargin{1};
    % Argument #2: VOXEL SIZE
    if (nargin < 2) || isempty(varargin{2})
        gEditMaskData.voxsize = [1 1 1];
    else
        gEditMaskData.voxsize = varargin{2};
    end
    % Argument #3: MASK
    if (nargin < 3) || isempty(varargin{3})
        gEditMaskData.mask = false(size(gEditMaskData.mri));
    elseif any(size(varargin{3}) ~= size(gEditMaskData.mri))
        error('MRI and mask must have the same dimensions.');
    else
        gEditMaskData.mask = logical(varargin{3});
    end
    % Argument #4: INITIAL POSITION
    if (nargin < 4) || isempty(varargin{4})
        InitialPos = [];
    else
        InitialPos = varargin{4};
    end
    % Argument #5: COLORMAP
    if (nargin < 5) || isempty(varargin{5})
        gEditMaskData.colormap = 'gray';
    else
        gEditMaskData.colormap = varargin{5};
    end

    % === START MASK EDITOR ===
    % Progress bar
    bst_progress('start', 'Edit mask', 'Initialization...');
    % Call again function, to initialize GUI
    hFig = mri_editMask();
    % Define Mouse wheel callback separately (not supported by old versions of Matlab)
    if isprop(hFig, 'WindowScrollWheelFcn')
        set(hFig, 'WindowScrollWheelFcn',  @FigureMouseWheelCallback);
    end
    % Set initial orientation and position
    if ~isempty(InitialPos)
        SetOrientation(InitialPos(1));
        SetSliceLocation(InitialPos(2));
    end
    % Close progress bar
    bst_progress('stop');

    % === WAIT AND RETURN ===
    % Wait for end of edition
    waitfor(hFig);
    % Return edited mask
    if (nargout >= 1)
        varargout{1} = gEditMaskData.mask;
    end
    % Empty global variable
    gEditMaskData = [];

% === INTERNAL CALL : CREATE GUI ===
elseif (nargin == 0)
    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @figureEditMask_OpeningFcn, ...
                       'gui_OutputFcn',  @figureEditMask_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % Update figures layout
    gui_layout('Update');

% === CALLBACKS ===
else
    if (nargout)
        [varargout{1:nargout}] = bst_call(str2func(varargin{1}), varargin{2:end});
    else
        bst_call(str2func(varargin{1}), varargin{2:end});
    end
end

%% ===== OUTPUT FUNCTION =====
function varargout = figureEditMask_OutputFcn(hObject, eventdata, handles)
    % Get default command line output from handles structure
    varargout{1} = hObject;

%% ===== OPENING FUNCTION =====
function figureEditMask_OpeningFcn(hObject, eventdata, handles, varargin)
    global gEditMaskData;
    % Check that MRI and MASK where defined
    if ~isfield(gEditMaskData, 'mri') || isempty(gEditMaskData.mri) || ~isfield(gEditMaskData, 'mask') || isempty(gEditMaskData.mask)
        error('Usage: mri_editMask(volume,mask);');
    end
    drawnow

    % === Replace Matlab slider by Java slider ===
    % Disable the Java-related warnings after 2019b
    if (bst_get('MatlabVersion') >= 907)
        warning('off', 'MATLAB:ui:javacomponent:FunctionToBeRemoved');
    end
%     sliderPos = get(handles.sliderPreview, 'Position');
%     hParent = get(handles.sliderPreview, 'Parent');
    delete(handles.sliderPreview);
    [handles.jSliderPreview, handles.sliderPreview] = javacomponent(javax.swing.JSlider(0,10,0), [0 0 1 1], hObject);
%     set(handles.sliderPreview, 'Units', 'Pixels', 'Position', sliderPos);
    java_setcb(handles.jSliderPreview, 'MouseReleasedCallback', @(h,ev)SetSliceLocation(ev.getSource().getValue()), ...
                                       'KeyPressedCallback',    @(h,ev)SetSliceLocation(ev.getSource().getValue()));

    % === REFRESH DISPLAY ===
    set(handles.mri_editMask, 'WindowStyle', 'normal');

    % ===== ICONS =====
    set(handles.toggleLassoMinus,   'CData', java_geticon('ICON_LASSO_MINUS'));
    set(handles.toggleLassoPlus,    'CData', java_geticon('ICON_LASSO_PLUS'));
    set(handles.toggleCurveMinus,   'CData', java_geticon('ICON_CURVE_MINUS'));
    set(handles.toggleCurvePlus,    'CData', java_geticon('ICON_CURVE_PLUS'));
    set(handles.buttonContrastDown, 'CData', java_geticon('ICON_BRIGHTNESS_MINUS'));
    set(handles.buttonContrastUp,   'CData', java_geticon('ICON_BRIGHTNESS_PLUS'));

    % Store handles in global variable
    handles.output = hObject;
    gEditMaskData.Handles = handles;
    % Keyboard shortcuts
    set([handles.mri_editMask, handles.buttonContrastDown, handles.buttonContrastUp, ...
         handles.buttonOk, handles.buttonCancel, ...
         handles.buttonAxial, handles.buttonCoronal, handles.buttonSagittal, ...
         handles.toggleLassoMinus, handles.toggleLassoPlus, handles.toggleCurveMinus, handles.toggleCurvePlus], ...
        'KeyPressFcn', @figureKeyPress_Callback);
    % Initial preview axes settings
    mriSize = size(gEditMaskData.mri);
    gEditMaskData.sliderValue = [round(mriSize(1)/2), round(mriSize(2)/2), round(mriSize(3)/2)];
    gEditMaskData.orientation = 3; % 1=x, 2=y, 3=z
    gEditMaskData.overlayColor = [1 1 0];
    gEditMaskData.overlayAlpha = .6;
    % Initial action: lasso
    SetAction('lasso+');
    % Mouse clicks and motion management
    gEditMaskData.isClicked.source = '';
    gEditMaskData.isClicked.clickPositionAxes = [];
    gEditMaskData.isClicked.clickPositionFigure = [];
    gEditMaskData.isClicked.roiPolygon = [];
%     % Initial action
%     set(handles.buttongroupAction, 'SelectedObject', handles.toggleLassoPlus);
%     % Update icons
%     gui_update_toggle([handles.toggleLassoPlus, handles.toggleLassoMinus, handles.toggleCurvePlus, handles.toggleCurveMinus]);
    % Initial contrast
    gEditMaskData.colormapMax = 256;
    gEditMaskData.intensityMax = max(gEditMaskData.mri(:));
    % Initialize axes
    SetOrientation(gEditMaskData.orientation);
    % Draw first slice
    UpdateSlice();

%% ===== RESIZE CALLBACK =====
function mri_editMask_ResizeFcn(varargin) %#ok<DEFNU>
    global gEditMaskData;
    if ~isfield(gEditMaskData, 'Handles') || ~ishandle(gEditMaskData.Handles.mri_editMask)
        return
    end
    handles = gEditMaskData.Handles;
    % Get figure and panels position and size
    figPos    = get(handles.mri_editMask,  'Position');
    rightPos  = get(handles.panelRight,    'Position');
    % Set right panel position
    set(handles.panelRight, 'Position', [0, ...
                                         figPos(4) - rightPos(4), ...
                                         rightPos(3), ...
                                         rightPos(4)]);
    % Set "MRI" panel position
    margin = 5;
    set(handles.panelMri, 'Position', [rightPos(3), ...
                                       margin, ...
                                       figPos(3) - rightPos(3) - margin, ...
                                       figPos(4) - 2 * margin]);
    % Set slider position
    set(handles.sliderPreview, 'Position', [11, figPos(4)-290, 140, 24]);

%% =================================================================================================
%  ===== DISPLAY MRI ===============================================================================
%  =================================================================================================
%% ===== SET ORIENTATION =====
function SetOrientation(orientation)
    global gEditMaskData;
    handles = gEditMaskData.Handles;
    % Save new orientation
    gEditMaskData.orientation = orientation;
    % Select right orientation button
    hButtons = [handles.buttonSagittal, handles.buttonCoronal, handles.buttonAxial];
    set(hButtons(orientation), 'Value', 1);

    % === PREPARE DISPLAY ===
    cubeSize = size(gEditMaskData.mri);
    FOV = cubeSize .* gEditMaskData.voxsize;
    % Empty axes
    cla(handles.axesPreview);
    % Set colormap
    colormap(handles.axesPreview, gEditMaskData.colormap);

    % === CREATE IMAGES ===
    % Get image size
    imgSize = cubeSize;
    imgSize(orientation) = [];
    FOV(orientation) = [];
    % MRI image
    handles.img_mri = image('XData',        [1, FOV(1)], ...
                            'YData',        [1, FOV(2)], ...
                            'CData',        zeros(imgSize(2), imgSize(1)), ...
                            'CDataMapping', 'direct', ...
                            'Parent',       gEditMaskData.Handles.axesPreview);
    % Overlay image
    handles.img_overlay = image('XData',     [1, FOV(1)], ...
                                'YData',     [1, FOV(2)], ...
                                'CData',     zeros([imgSize(2), imgSize(2), 3]), ...
                                'AlphaData', 0, ...
                                'Parent',    gEditMaskData.Handles.axesPreview);
    % Axes
    axis(handles.axesPreview, 'image', 'off');
    % Set mouse callbacks
    set([handles.axesPreview, handles.img_mri, handles.img_overlay], 'ButtonDownFcn', @previewMouseClick_Callback);
    % Save handles structure
    gEditMaskData.Handles = handles;

    % === ORIENTATION MARKERS ===
    if (orientation == 2)
        strOrient = 'PA';
    else
        strOrient = 'LR';
    end
    % Coronal
    V = axis(handles.axesPreview);
    fontSize = bst_get('FigFont');
    text(       6, 15, strOrient(1), 'verticalalignment', 'top', 'FontSize', fontSize, 'FontUnits', 'points', 'color','w', 'Parent', gEditMaskData.Handles.axesPreview);
    text(.95*V(2), 15, strOrient(2), 'verticalalignment', 'top', 'FontSize', fontSize, 'FontUnits', 'points', 'color','w', 'Parent', gEditMaskData.Handles.axesPreview);

    % === SLIDER ===
    % Slider bounds
    handles.jSliderPreview.setMinimum(1);
    handles.jSliderPreview.setMaximum(cubeSize(orientation));
    % Set slice location
    SetSliceLocation(gEditMaskData.sliderValue(orientation));

    % === COLORMAP ===
    UpdateColormap();

%% ===== SET SLICE LOCATION =====
function SetSliceLocation(sliceLoc)
    global gEditMaskData;
    handles = gEditMaskData.Handles;
    % Force location inside the MRI volume
    sliceLoc = bst_saturate(sliceLoc, [1, size(gEditMaskData.mri, gEditMaskData.orientation)]);
    % Save new position
    gEditMaskData.sliderValue(gEditMaskData.orientation) = sliceLoc;
    % Update slider position
    handles.jSliderPreview.setValue(sliceLoc);
    % Reset points selection
    gEditMaskData.isClicked.roiPolygon = [];
    delete(findobj(handles.axesPreview, 'Tag', 'ptCurve'));
    % Update text
    set(handles.textValue, 'String', sprintf('%d', sliceLoc));
    strOrient = {'x=','y=','z='};
    set(handles.textTitleValue, 'String', strOrient{gEditMaskData.orientation});
    % Update location axes
    UpdateLocationAxes();
    % Redraw slice
    UpdateSlice();

%% ===== UPDATE SLICE =====
function UpdateSlice()
    global gEditMaskData;
    % Get slice location
    sliceLoc = gEditMaskData.sliderValue(gEditMaskData.orientation);
    % Update MRI slice
    sliceMRI = mri_getslice(gEditMaskData.mri, sliceLoc, gEditMaskData.orientation)';
    set(gEditMaskData.Handles.img_mri, 'CData', double(sliceMRI) ./ double(gEditMaskData.intensityMax) .* 256);
    % Update overlay slice
    if ~isempty(gEditMaskData.mask)
        sliceOverlay = mri_getslice(gEditMaskData.mask, sliceLoc, gEditMaskData.orientation)';
        sliceAlpha   = .6 * double(sliceOverlay);
        sliceOverlay = cat(3, sliceOverlay * gEditMaskData.overlayColor(1), ...
                              sliceOverlay * gEditMaskData.overlayColor(2), ...
                              sliceOverlay * gEditMaskData.overlayColor(3));
        set(gEditMaskData.Handles.img_overlay, 'CData', double(sliceOverlay), ...
                                               'AlphaData', sliceAlpha);
    end
    % Update location line
    DrawLocationLine();
    drawnow();

%% =================================================================================================
%  ===== LOCATION FUNCTIONS ========================================================================
%  =================================================================================================
%% ===== UPDATE LOCATION AXES =====
function UpdateLocationAxes()
    global gEditMaskData;
    % Get a centered slice, orthogonal to the direction of the preview slice
    switch (gEditMaskData.orientation)
        case 1 % Coronal preview => extract a sagittal slice
            locOrient = 3;
            lineOrient = 1;
        case 2 % Axial preview => extract a sagittal slice
            locOrient = 3;
            lineOrient = 2;
        case 3 % Sagittal preview => extract an axial slice
            locOrient = 2;
            lineOrient = 2;
    end
    % Get slice at the middle of the volume
    cubeSize = size(gEditMaskData.mri);
    locPosition = round(cubeSize(locOrient)/2);
    locSlice = mri_getslice(gEditMaskData.mri, locPosition, locOrient);
    % Get pixel size
    FOV = cubeSize .* gEditMaskData.voxsize;
    FOV(locOrient) = [];

    % Empty axes
    cla(gEditMaskData.Handles.axesLocation);
    % Display slice on location axes
    gEditMaskData.img_location = image('XData',         [1, FOV(1)], ...
                                       'YData',         [1, FOV(2)], ...
                                       'CData',         double(locSlice)' ./ double(gEditMaskData.intensityMax) .* 256, ...
                                       'CDataMapping',  'scaled', ...
                                       'Parent',        gEditMaskData.Handles.axesLocation, ...
                                       'ButtonDownFcn', @locationMouseClick_Callback);
    % Configure axes
    axis(gEditMaskData.Handles.axesLocation, 'image', 'off');
    set(gEditMaskData.Handles.axesLocation, 'UserData', lineOrient);

%% ===== DRAW LOCATION LINE =====
function DrawLocationLine()
    global gEditMaskData;
    % Get line orientation (UserData of the location Axes)
    lineOrient = get(gEditMaskData.Handles.axesLocation, 'UserData');
    % Get values from interface controls
    n = gEditMaskData.sliderValue(gEditMaskData.orientation) * gEditMaskData.voxsize(gEditMaskData.orientation);
    XLim = get(gEditMaskData.Handles.axesLocation, 'XLim');
    YLim = get(gEditMaskData.Handles.axesLocation, 'YLim');
    % Display a line that represent the location of the slice
    delete(findobj(gEditMaskData.Handles.axesLocation, 'Tag', 'locationLine'));
    switch(lineOrient)
        case 1 % Vertical line
            line([n,n]+.5, YLim+.5, [1,1], 'Color', [1 .2 .2], ...
                 'Parent', gEditMaskData.Handles.axesLocation, 'Tag', 'locationLine');
        case 2 % Horizontal line
            line(XLim+.5, [n,n], [1,1], 'Color', [1 .2 .2], ...
                 'Parent', gEditMaskData.Handles.axesLocation, 'Tag', 'locationLine');
    end
    drawnow;

%% ===== LOCATION AXES CALLBACK =====
function locationMouseClick_Callback(varargin)
    global gEditMaskData;
    % Get application handles
    handles = gEditMaskData.Handles;
    % Get line orientation (UserData of the location Axes)
    lineOrient = get(handles.axesLocation, 'UserData');
    % Mouse location on the axis
    temp = get(handles.axesLocation, 'CurrentPoint');
    curptAxes = temp(1,1:2);
    % Update slice location
    switch(lineOrient)
        case 1 % Vertical line
            n = round(curptAxes(1) * gEditMaskData.voxsize(1) - .5);
        case 2 % Horizontal line
            n = round(curptAxes(2) * gEditMaskData.voxsize(2) - .5);
    end
    SetSliceLocation(n);

%% =================================================================================================
%  ===== AXES CALLBACKS ============================================================================
%  =================================================================================================
%% ===== KEYBOARD CALLBACK =====
function figureKeyPress_Callback(hObject, eventdata, varargin)
    global gEditMaskData;
    % Get GUI handles
    handles = gEditMaskData.Handles;
    % If the focus is on the slider, set the focus on the axes, to avoid conflits between uicontrol
    % automatic keyboard handling and, and this manual handling
    if (get(handles.mri_editMask, 'CurrentObject') == handles.sliderPreview)
        uicontrol(handles.buttonAxial);
    end
    % Get current slice location
    sliceLoc = handles.jSliderPreview.getValue();
    % Switch between keys
    switch (eventdata.Key)
        % PAGE UP : increase slice indice (10)
        case 'pageup'
            SetSliceLocation(sliceLoc + 10);
        % PAGE DOWN : decrease slice indice (10)
        case 'pagedown'
            SetSliceLocation(sliceLoc - 10);
        % RIGHT ARROW : increase slice indice (1)
        case 'rightarrow'
            SetSliceLocation(sliceLoc + 1);
        % LEFT ARROW : decrease slice indice
        case 'leftarrow'
            SetSliceLocation(sliceLoc - 1);
    end

%% ===== MOUSE CLICK =====
function previewMouseClick_Callback(varargin)
    global gEditMaskData;
    % Get figure handles
    hFig  = gEditMaskData.Handles.mri_editMask;
    hAxes = gEditMaskData.Handles.axesPreview;
    % Mouse location on the figure
    gEditMaskData.isClicked.source = 'preview';
    temp = get(hFig, 'CurrentPoint');
    gEditMaskData.isClicked.clickPositionFigure = temp(1,1:2);
    % Mouse location on the axis
    temp = get(hAxes, 'CurrentPoint');
    curptAxes= temp(1,1:2);
    % Apply changes
    switch(gEditMaskData.action)
        case {'lasso+', 'lasso-'}
            % Reset selection
            gEditMaskData.isClicked.roiPolygon = [];
            % Add first point to the ROI polygon
            gEditMaskData.isClicked.roiPolygon(1,1:2) = curptAxes;
    end

%% ===== MOUSE MOVE =====
function figureMouseMove_Callback(varargin) %#ok<DEFNU>
    global gEditMaskData;
    if ~isfield(gEditMaskData, 'Handles')
        return;
    end
    handles = gEditMaskData.Handles;
    % Mouse location on the figure
    temp = get(handles.mri_editMask, 'CurrentPoint');
    curptFigure = temp(1,1:2);
    % Mouse location on the axis
    temp = get(handles.axesPreview, 'CurrentPoint');
    curptAxes= temp(1,1:2);

    switch(gEditMaskData.isClicked.source)
        case 'preview'
            % Mouse motion
            mouseMotionFigure = gEditMaskData.isClicked.clickPositionFigure - curptFigure;
            gEditMaskData.isClicked.clickPositionFigure = curptFigure;
            % Mouse selection type ?
            switch(lower(get(handles.mri_editMask, 'SelectionType')))
                % Normal left click : Currently selected action
                case 'normal'
                    switch(gEditMaskData.action)
                        case {'lasso-', 'lasso+'}
                            % Add current point to the ROI polygon
                            n_1 = size(gEditMaskData.isClicked.roiPolygon,1);
                            gEditMaskData.isClicked.roiPolygon(n_1+1,1:2) = curptAxes;
                            % Draw line between current polygon vertex and last polygon vertex
                            line('XData',[gEditMaskData.isClicked.roiPolygon(n_1,1), gEditMaskData.isClicked.roiPolygon(n_1+1,1)], ...
                                 'YData',[gEditMaskData.isClicked.roiPolygon(n_1,2), gEditMaskData.isClicked.roiPolygon(n_1+1,2)], ...
                                 'ZData',[.1 .1], 'Color', [1 .2 .2], ...
                                 'Parent', handles.axesPreview, ...
                                 'Tag', 'ptCurve', ...
                                 'ButtonDownFcn', @previewMouseClick_Callback);
                    end

                % (Control + Left click) or right click : Change color scale
                case 'alt'
                    % Increase/decrease image contrast
                    delete(findobj(handles.axesPreview, 'Tag', 'selectionArea'));
                    gEditMaskData.colormapMax = max(8, min(505, gEditMaskData.colormapMax + mouseMotionFigure(2)));
                    UpdateColormap();
            end
    end
    drawnow

%% ===== MOUSE UP =====
function figureMouseUp_Callback(varargin) %#ok<DEFNU>
    global gEditMaskData;
    handles = gEditMaskData.Handles;
    gEditMaskData.isClicked.source = '';
    % Mouse location on the axes
    temp = get(handles.axesPreview, 'CurrentPoint');
    curptAxes= temp(1,1:2);
    % Add or remove new mask
    isAddNewMask = ~isempty(strfind(gEditMaskData.action, '+'));
    % If the mouse pointer is on the preview axes
    switch(gEditMaskData.action)
        case {'lasso+', 'lasso-'}
            % Close the polygon and display the polygon mask
            n = size(gEditMaskData.isClicked.roiPolygon,1);
            if (n > 1)
                gEditMaskData.isClicked.roiPolygon(n+1,1:2) = gEditMaskData.isClicked.roiPolygon(1,1:2);
                % Extract the binary mask from the polygon
                roiMask = getMask(gEditMaskData.isClicked.roiPolygon);
                % Update mask volume
                AddMaskSlice(roiMask, isAddNewMask);
                % Display mask border
                delete(findobj(handles.axesPreview, 'Tag', 'ptCurve'));
            end

        case {'curve+', 'curve-'}
            % Left button : add point
            % Right button : add point and end curve
            % Get new point
            newPt = curptAxes;
            % Add it to list of points
            n = size(gEditMaskData.isClicked.roiPolygon,1);
            gEditMaskData.isClicked.roiPolygon(n+1,1:2) = newPt;
            % Plot point
            hold on
            plot3(handles.axesPreview, newPt(1), newPt(2), 0.2, 'r.', 'Tag', 'ptCurve');
            % Only for right button
            if strcmpi(get(handles.mri_editMask, 'SelectionType'), 'alt') && (n > 1)
                % Add first point to close spline
                gEditMaskData.isClicked.roiPolygon(n+2,1:2) = gEditMaskData.isClicked.roiPolygon(1,1:2);
                % Interpolate with a spline curve and finer spacing.
                t = 1:n+2;
                ts = 1:0.1:n+2;
                splineRoi = spline(t, gEditMaskData.isClicked.roiPolygon', ts);
                % Delete all curve points
                delete(findobj(handles.axesPreview, 'Tag', 'ptCurve'));
                % Plot the interpolated closed curve
                line(splineRoi(1,:), splineRoi(2,:), ones(1,size(splineRoi,2)) * 0.2, ...
                     'Color',     [1 0 0], ...
                     'LineStyle', '-', ...
                     'Tag',       'ptCurve', ...
                     'Parent',    handles.axesPreview);
                % Extract the binary mask from the polygon
                roiMask = getMask(splineRoi');
                % Update mask volume
                AddMaskSlice(roiMask, isAddNewMask);
                % Reset polygon selection
                gEditMaskData.isClicked.roiPolygon = [];
            end
    end

%% ===== MOUSE WHEEL =====
% Move slices
function FigureMouseWheelCallback(hObject, event, varargin)
    global gEditMaskData;
    % Get new slice index
    iSlice = gEditMaskData.sliderValue(gEditMaskData.orientation) - double(event.VerticalScrollCount);
    % Go to new slice
    SetSliceLocation(iSlice);

%% ===== ACTION BUTTONS CALLBACKS =====
function buttongroupAction_SelectionChangeFcn(varargin) %#ok<DEFNU>
    global gEditMaskData;
    handles = gEditMaskData.Handles;
    % Get selected tool button
    selButton = get(handles.buttongroupAction, 'SelectedObject');
    % If no button is selected, for selection of "Lasso+"
    if isempty(selButton)
        selButton = handles.toggleLassoPlus;
        set(handles.buttongroupAction, 'SelectedObject', selButton);
    end
    % Associate an action name with each button
    switch (selButton)
        case handles.toggleLassoMinus
            action = 'lasso-';
        case handles.toggleLassoPlus
            action = 'lasso+';
        case handles.toggleCurveMinus
            action = 'curve-';
        case handles.toggleCurvePlus
            action = 'curve+';
    end
    SetAction(action);

%% ===== SET ACTION =====
function SetAction(action)
    global gEditMaskData;
    handles = gEditMaskData.Handles;
    % Set action
    switch lower(action)
        case 'lasso+'
            set(handles.toggleLassoPlus, 'Value', 1);
        case 'lasso-'
            set(handles.toggleLassoMinus, 'Value', 1);
        case 'curve+'
            set(handles.toggleCurvePlus, 'Value', 1);
        case 'curve-'
            set(handles.toggleCurveMinus, 'Value', 1);
    end
    gEditMaskData.action = action;
    % Update icons
    gui_update_toggle([handles.toggleLassoPlus, handles.toggleLassoMinus, handles.toggleCurvePlus, handles.toggleCurveMinus]);
    % Reset polygon selection
    gEditMaskData.isClicked.roiPolygon = [];
    delete(findobj(handles.axesPreview, 'Tag', 'ptCurve'));

%% ===== GET POLYGON =====
function roiMask = getMask(roiPolygon)
    global gEditMaskData;
    % Get slice size
    sliceSize = size(gEditMaskData.mri);
    sliceSize(gEditMaskData.orientation) = [];
    % Get pixel size
    pixSize = gEditMaskData.voxsize;
    pixSize(gEditMaskData.orientation) = [];
    % Build full coordinates of each pixel
    sliceX = meshgrid(1:sliceSize(1),1:sliceSize(2));
    sliceX = sliceX(:);
    sliceY = meshgrid(1:sliceSize(2),1:sliceSize(1))';
    sliceY = sliceY(:);
    % Test which pixels are inside the polygon
    iInsidePixels = inpolygon(sliceX(:), sliceY, roiPolygon(:,1) / pixSize(1), roiPolygon(:,2) / pixSize(2));
    iInsidePixels = sub2ind(sliceSize, round(sliceX(iInsidePixels)), round(sliceY(iInsidePixels)));
    % Build mask
    roiMask = false(sliceSize);
    roiMask(iInsidePixels) = 1;

%% ===== ADD MASK SLICE =====
function AddMaskSlice(sliceMask, isAddNewMask)
    % Update volume
    global gEditMaskData;
    n = gEditMaskData.sliderValue(gEditMaskData.orientation);
    % Get previous slice from mask
    initSlice = mri_getslice(gEditMaskData.mask, n, gEditMaskData.orientation);
    % Add or remove
    if isAddNewMask
        sliceMask = (initSlice | sliceMask);
    else
        sliceMask = (initSlice & ~sliceMask);
    end
    % Depends on orientation
    switch (gEditMaskData.orientation)
        case 1
            gEditMaskData.mask(n,:,:) = sliceMask;
        case 2
            gEditMaskData.mask(:,n,:) = sliceMask;
        case 3
            gEditMaskData.mask(:,:,n) = sliceMask;
    end
    % Update display
    UpdateSlice();

%% =================================================================================================
%  ===== COLORMAP AND CONTRAST =====================================================================
%  =================================================================================================
%% ===== INCREASE CONTRAST =====
function buttonContrastUp_Callback(varargin) %#ok<DEFNU>
    global gEditMaskData;
    gEditMaskData.colormapMax = max(8, gEditMaskData.colormapMax - 20);
    UpdateColormap();

%% ===== DECREASE CONTRAST =====
function buttonContrastDown_Callback(varargin) %#ok<DEFNU>
    global gEditMaskData;
    gEditMaskData.colormapMax = min(505, gEditMaskData.colormapMax + 20);
    UpdateColormap();

%% ===== UPDATE COLORMAP =====
function UpdateColormap()
    global gEditMaskData;
    mapName = gEditMaskData.colormap;
    % Between 0 and 256 : fill the map with the top color of the colormap (bright)
    if (gEditMaskData.colormapMax <= 256)
        map = eval(sprintf('%s(%d)', lower(mapName), gEditMaskData.colormapMax));
        nmap = repmat(map(gEditMaskData.colormapMax,:),256,1);
        nmap(1:gEditMaskData.colormapMax, :) = map;
    % Between 256 and 512 : fill the map with the bottom color of the colormap (dark)
    else
        map = eval(sprintf('%s(%d)', lower(mapName), 256 - mod(gEditMaskData.colormapMax, 256)));
        nmap = repmat(map(1,:),256,1);
        nmap(257-length(map):256, :) = map;
    end
    set(gEditMaskData.Handles.mri_editMask, 'Colormap', nmap);

%% =================================================================================================
%  ===== VALIDATION BUTTONS ========================================================================
%  =================================================================================================
%% ===== CANCEL BUTTON =====
function buttonCancel_Callback(varargin) %#ok<DEFNU>
    global gEditMaskData;
    % Return an empty mask
    gEditMaskData.mask = [];
    % Delete window
    delete(gEditMaskData.Handles.mri_editMask);
    drawnow;

%% ===== OK BUTTON =====
function buttonOk_Callback(varargin)  %#ok<DEFNU>
    global gEditMaskData;
    % Delete window
    delete(gEditMaskData.Handles.mri_editMask);
    drawnow;

author: brainstorm-tools


text:
function s=maxmod(a,b)
% The maxmod limiter
s=0.5*(sign(a)+sign(b)).*max(abs(a),abs(b));

repo: Aero-matlab


text:
function cellInfo = Mission_Data(varargin)
% MISSION_DATA returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'Mission_Data_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('ms'), ''}; ...
{'valid_items', 1, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'reserved', 1, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'seq', 8, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Start from 0')}; ...
{'command', 8, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'frame', 8, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'current', 8, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'autocontinue', 8, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'mission_type', 8, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'param1', 8, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'param2', 8, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'param3', 8, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'param4', 8, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'x', 8, 'int32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'y', 8, 'int32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'z', 8, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function engine = dv_unrolled_dbn_inf_engine(bnet, T, varargin)
% JTREE_UNROLLED_DBN_INF_ENGINE Unroll the DBN for T time-slices and apply jtree to the resulting static net
% engine = jtree_unrolled_dbn_inf_engine(bnet, T, ...)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% useC      - 1 means use jtree_C_inf_engine instead of jtree_inf_engine [0]
% constrained - 1 means we constrain ourselves to eliminate slice t before t+1 [1]
%
% e.g., engine = jtree_unrolled_inf_engine(bnet, 'useC', 1);

% set default params
N = length(bnet.intra);
useC = 0;
constrained = 1;

if nargin >= 3
  args = varargin;
  nargs = length(args);
  if isstr(args{1})
    for i=1:2:nargs
      switch args{i},
       case 'useC',   useC = args{i+1};
       case 'constrained',  constrained = args{i+1};
       otherwise,
	error(['invalid argument name ' args{i}]);
      end
    end
  else
    error(['invalid argument name ' args{1}]);
  end
end

bnet2 = hodbn_to_bnet(bnet, T);
ss = length(bnet.intra);
engine.ss = ss;

% If constrained_order = 1 we constrain ourselves to eliminate slice t before t+1.
% This prevents cliques containing nodes from far-apart time-slices.
if constrained
  stages = num2cell(unroll_set(1:ss, ss, T), 1);
else
  stages = { 1:length(bnet2.dag) };
end
if useC
  %jengine = jtree_C_inf_engine(bnet2, 'stages', stages);
  %function is not implemented
  assert(0)
else
  jengine = stab_cond_gauss_inf_engine(bnet2);
end

engine.unrolled_engine = jengine;
% we don't inherit from jtree_inf_engine, because that would only store bnet2,
% and we would lose access to the DBN-specific fields like intra/inter

engine.nslices = T;
engine = class(engine, 'stable_ho_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Template script for FieldMap toolbox
% _____________________________________________________________________
%
% This script gives an example of how to run the FieldMap toolbox
% FieldMap.m without using the GUI. It can be expanded using standard
% matlab code to for example create multiple field maps or use a single
% fieldmap to unwarp multiple images.
%
% As it stands the script uses routines from the FieldMap toolbox to
% create a single field map which is matched to an EPI and then used
% to unwarp it. A structural image is loaded and matched to the unwarped
% EPI.
%
% For details about the FieldMap toolbox, see FieldMap.man. For a
% description of the components of the structure IP, see FieldMap.m.
% For an introduction to the theoretcial and practical principles behind
% the toolbox, see principles.man.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson and Chloe Hutton
% $Id: FieldMap_ngui.m 1358 2008-04-10 11:20:26Z guillaume $

%----------------------------------------------------------------------
% Set up default parameters and structures
%----------------------------------------------------------------------

spm('defaults','FMRI');
IP = FieldMap('Initialise'); % Gets default params from pm_defaults

%----------------------------------------------------------------------
% Load measured field map data - phase and magnitude or real and imaginary
%----------------------------------------------------------------------

if IP.uflags.iformat=='PM'
   for index=1:4
      IP.P{index} = FieldMap('LoadFilePM',index);
   end
else
   for index=1:4
      IP.P{index} = FieldMap('LoadFile',index);
   end
end

%----------------------------------------------------------------------
% Or you may want to load a precalculated Hz phase map instead...
%----------------------------------------------------------------------

% [IP.fm, IP.pP] = FieldMap('LoadFieldMap');

%----------------------------------------------------------------------
% Create field map (in Hz) - this routine calls the unwrapping
%----------------------------------------------------------------------

IP.fm = FieldMap('CreateFieldMap',IP);

%----------------------------------------------------------------------
% Write out field map
% Outputs -> fpm_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------

FieldMap('Write',IP.P{1},IP.fm.fpm,'fpm_',64,'Smoothed phase map');

%----------------------------------------------------------------------
% Convert Hz to voxels and write voxel displacement map
% Outputs -> vdm_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------

[IP.vdm, IP.vdmP]=FieldMap('FM2VDM',IP);

%----------------------------------------------------------------------
% Select an EPI to unwarp
%----------------------------------------------------------------------

IP.epiP = FieldMap('LoadEPI');

%----------------------------------------------------------------------
% Match voxel displacement map to image
% Outputs -> mag_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------

IP.vdmP = FieldMap('MatchVDM',IP);

%----------------------------------------------------------------------
% Unwarp EPI
%----------------------------------------------------------------------

IP.uepiP = FieldMap('UnwarpEPI',IP);

%----------------------------------------------------------------------
% Write unwarped EPI
% Outputs -> uNAME-OF-EPI.img
%----------------------------------------------------------------------

IP.uepiP = FieldMap('Write',IP.epiP,IP.uepiP.dat,'u',IP.epiP.dim(4),'Unwarped image');

%----------------------------------------------------------------------
% Load a structural image
%----------------------------------------------------------------------

IP.nwarp=FieldMap('LoadStructural');

%----------------------------------------------------------------------
% Coregister structural with the unwarped image
%----------------------------------------------------------------------

FieldMap('MatchStructural',IP);

%______________________________________________________________________

author: spm


text:
function cab(varargin)

% This function closes all figures currently open EXCEPT for
% those listed as arguments.  'cab' stands for 'close all but'.
%
% Usage:
%       cab figure_handle1 figure_handle2 ...
%       cab(figure_handle1, figure_handle2, ...)
%       cab('last') % or also:  cab last
%
%   - The 'last' option closes all figures except the last one opened.
%   - Calling 'cab' with no arguments is a convenient
%     alternative to 'close all'
%
% Example:
%   figure(5)
%   figure(7)
%   figure(9)
%   figure(11)
%   cab(7, 11)  % or also:  cab 7 11

% all_figs = findall(0, 'type', 'figure');  % Uncomment this to include ALL windows, including those with hidden handles (e.g. GUIs)
all_figs = findobj(0, 'type', 'figure');
figs2keep = [];
for i = 1:nargin
    if ischar(varargin{i})
        if strcmp(varargin{i}, 'last')
            figs2keep = all_figs(1);
            %figs2keep = gcf;
        else
            % In this case, function was called as follows:  cab 1 2 3
            figs2keep = [figs2keep, str2num(varargin{i})];
        end
    else
        % In this case, function was called as follows:  cab(1, 2, 3)
        figs2keep = [figs2keep, varargin{i}];
    end
end

delete(setdiff(all_figs, figs2keep))

author: Sable

repo: mcbench-benchmarks


text:
% 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/>.

cd(getCERRPath)
cd compiled
mcc -B sgl CERR.m
delete *.h
delete *.c

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = ProxLP(y,p,tau)

if length(y(:))>1
    for i=1:length(y(:))
        x(i)=ProxLP(y(i),p,tau);
    end
    x = reshape(x,size(y));
    return
end

x = linspace(0,abs(y)*1.1,200);
[~,i] = min( 1/2*(x-abs(y)).^2 + tau*abs(x).^p );
x = x*sign(x(i));

end

author: optimaltransport


text:
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( x < 0.0 )
    value = -1.0;
  else
    value = 1.0;
  end

  return
end


text:
%APPLYCOLORMAP  Applies a GNU Octave/MATLAB equivalent colormap on a given image
%
%     dst = cv.applyColorMap(src, colormap)
%     dst = cv.applyColorMap(src, userColor)
%     dst = cv.applyColorMap(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ The source image, grayscale or color of type `uint8`.
% * __colormap__ The colormap to apply. A string, one of:
%   * __Autumn__ Shades of red and yellow color map.
%   * __Bone__ Gray-scale with a tinge of blue color map.
%   * __Jet__ Variant of HSV.
%   * __Winter__ Shades of blue and green color map.
%   * __Rainbow__ Red-orange-yellow-green-blue-violet color map.
%   * __Ocean__ Black to white with shades of blue color map.
%   * __Summer__ Shades of green and yellow colormap.
%   * __Spring__ Shades of magenta and yellow color map.
%   * __Cool__ Shades of cyan and magenta color map.
%   * __HSV__ Hue-saturation-value color map.
%   * __Pink__ Pastel shades of pink color map.
%   * __Hot__ Black-red-yellow-white color map.
%   * __Parula__ Blue-green-orange-yellow color map.
% * __userColor__ The colormap to apply of type `uint8` (1 or 3 channels) and
%   length 256.
%
% ## Output
% * __dst__ The result is the colormapped source image. Same row/column size
%   and same type as `src`.
%
% ## Options
% * __FlipChannels__ whether to flip the order of color channels in output
%   `dst`, from OpenCV's BGR to between MATLAB's RGB. default true
%
% The human perception isn't built for observing fine changes in grayscale
% images. Human eyes are more sensitive to observing changes between colors,
% so you often need to recolor your grayscale images to get a clue about them.
% OpenCV now comes with various colormaps to enhance the visualization in your
% computer vision application.
%
% The second variant of the function applies a user-defined colormap on the
% given image.
%
% ## Example
%
% In OpenCV you only need cv.applyColorMap to apply a colormap on a given
% image. The following sample code takes an image and applies a Jet colormap
% on it and shows the result:
%
%     % We need an input image. (can be grayscale or color)
%     img_in = im2uint8(mat2gray(peaks(500)));
%
%     % Apply the colormap
%     %img_color2 = im2uint8(ind2rgb(img_in, jet(256)));
%     img_color = cv.applyColorMap(img_in, 'Jet');
%
%     % Show the result
%     imshow(img_color)
%
% ## Example
%
%     cmaps = {'Autumn', 'Bone', 'Jet', 'Winter', 'Rainbow', 'Ocean', ...
%         'Summer', 'Spring', 'Cool', 'HSV', 'Pink', 'Hot', 'Parula'};
%     img = cell2mat(cellfun(@(cmap) ...
%         cv.applyColorMap(repmat(uint8(0:255), 20, 1), cmap), cmaps(:), ...
%         'UniformOutput',false));
%     image(img)
%     set(gca, 'YTick', 10:20:20*numel(cmaps), 'YTickLabel',cmaps)
%     title('Colormaps')
%
% See also: cv.LUT, ind2rgb, colormap, autumn, bone, jet, winter, summer,
%  spring, cool, hsv, pink, hot, parula
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function testHarness

%% Lane Marking Identification
% One component of future luxury automobile safety
% systems is warning drivers that they are drifting
% between lanes.  To do this, it is first necessary
% to identify where the lanes are.  In this example,
% we will look at identifying the lanes from a image
% of a road.
%
% Copyright 2007-2013 MathWorks, Inc.
%

%% Setup
clear

%% Initialize Objects
% Video Reader

hvfr = vision.VideoFileReader('viplanedeparture.avi', ...
    'VideoOutputDataType', 'uint8');

%% Configure Video Players
hVidSource = vision.VideoPlayer('Name', 'Source Image');

%% Setup Replay Loop
NumLoops = 1;
for loopIdx = 1:NumLoops

    % Loop through video frames
    while ~isDone(hvfr)
        % Read Frame
        frameRGB = step(hvfr);
        % and display
        step(hVidSource, frameRGB);

        % Call Lane Marking Algorithm
        laneMarkingAlgorithmWithVisualization(frameRGB);

    % End While
    end

    reset(hvfr);
% end loop
end

%% Release resources
release(hvfr   );


text:
function ue= f_ptos(signal1)
%  function to be used in the ptos program
% k1=100;%  N=1;
Na=signal1(1);
k1=signal1(2);
e=signal1(3);
k2=sqrt(2*k1/Na);
if  abs(e)<=1/k1
     ue=(k1/k2)*e;
else
    ue=sign(e)*(sqrt(2*Na*abs(e))- 1/k2);
end


text:
function SVP = svp(T)

% T = temperature on Kelvin scale
% SVP = saturation vapor pressure according to Goff & Gratch Equation (1946)
% SVP in Pascal

Ts = 373.15;    % standard temperature at steam point on Kelvin scale (Goff, 1965)
Ps = 101324.6;  % standard atmospheric pressure at steam point (Pascal)

log10SVP = -7.90298*(Ts./T - 1) + 5.02808*log10(Ts./T) -1.3816e-7*(10.^(11.344*(1 - T./Ts)) - 1) + 8.1328e-3*(10.^(3.49149*(1 - Ts./T)) - 1) + log10(Ps);

SVP = 10.^(log10SVP);


text:
function D = computeDistMatrix_AVFC (data,T,options)
%
% It uses KL-divergence to compute a subject-by-subject distance matrix
% of average FC matrices (on fMRI).
%
% INPUT
% data          observations; in this case it has to be a cell, each with
%               the data for one subject
% T             length of series, also a cell.
%
% OUTPUT
% D             (N by N) distance matrix, with the distance between each
%               pair of subjects in "HMM space"
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2020)

if ~iscell(data) || ~iscell(T), error('X and T must both be cells'); end

if nargin<3, options = struct(); end

N = length(data);

for n = 1:N
    if ischar(data{n})
        fsub = data{n};
        loadfile_sub;
    else
        X = data{n};
    end
    X = preprocdata(X,T{n},options);
    if n == 1
        ndim = size(X,2);
        V = zeros(ndim,ndim,N);
    end
    V(:,:,n) = X' * X;
end

D = NaN(N);
try
    for n1 = 1:N-1
        for n2 = n1+1:N
            D(n1,n2) =  ( wishart_kl(V(:,:,n1),V(:,:,n2),sum(T{n1}),sum(T{n2})) + ...
                wishart_kl(V(:,:,n2),V(:,:,n1),sum(T{n2}),sum(T{n1})) ) /2;
            D(n2,n1) = D(n1,n2);
        end
    end
catch
    for n = 1:N
        V(:,:,n) = V(:,:,n) + 0.0001*eye(ndim);
    end
    for n1 = 1:N-1
        for n2 = n1+1:N
            D(n1,n2) =  ( wishart_kl(V(:,:,n1),V(:,:,n2),sum(T{n1}),sum(T{n2})) + ...
                wishart_kl(V(:,:,n2),V(:,:,n1),sum(T{n2}),sum(T{n1})) ) /2;
            D(n2,n1) = D(n1,n2);
        end
    end
end

end

function X = preprocdata(X,T,options)

% Filtering
if isfield(options,'filter') && ~isempty(options.filter)
    X = filterX(X,T,options.Fs,options.filter);
end
% Detrend X
if isfield(options,'detrend') && options.detrend
    X = detrendX(X,T);
end
% Hilbert envelope
if isfield(options,'onpower') && options.onpower
    X = rawsignal2power(X,T);
end
% Embedding
if isfield(options,'embeddedlags') && length(options.embeddedlags) > 1
    X = embeddata(X,T,options.embeddedlags);
end
X = zscore(X);

end

author: OHBA-analysis

repo: HMM-MAR


text:
classdef GemDropDownMenu < GemUserInterfaceObject
    % GemDropDownMenu GEM class for a pop-up menu control
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Callback
        FontSize
        HorizontalAlignment
        MenuItems
        SelectedIndex
        Tag
        ToolTip
    end

    methods
        function obj = GemDropDownMenu(parent, menu_items, tooltip, callback)
            obj = obj@GemUserInterfaceObject(parent);

            obj.MenuItems = menu_items;
            obj.ToolTip = tooltip;
            obj.Callback = callback;
            obj.FontSize = 11;
            obj.HorizontalAlignment = 'left';
            obj.SelectedIndex = 1;
        end

        function delete(obj)
        end

        function CreateGuiComponent(obj, position)
            obj.GraphicalComponentHandle = uicontrol('Style', 'popupmenu', 'Parent', obj.Parent.GetContainerHandle, ...
                'FontAngle', 'normal', 'FontUnits', 'pixels', 'FontSize', obj.FontSize, 'HorizontalAlignment', obj.HorizontalAlignment, ...
                'Units', 'pixels', 'Position', position, 'ToolTipString', obj.ToolTip, 'String', obj.MenuItems, ...
                'Callback', @obj.PopupmenuCallback, 'Value', obj.SelectedIndex);

        end

        function SetMenuItems(obj, menu_items)
            obj.MenuItems = menu_items;
            if obj.ComponentHasBeenCreated
                set(obj.GraphicalComponentHandle, 'String', menu_items);
            end
        end

        function SetSelectedIndex(obj, selected_index)
            obj.SelectedIndex = selected_index;
            if obj.ComponentHasBeenCreated
                current_index = get(obj.GraphicalComponentHandle, 'Value');
                if (current_index ~= obj.SelectedIndex)
                    set(obj.GraphicalComponentHandle, 'Value', obj.SelectedIndex);
                end
            end
        end

    end

    methods (Access = protected)

        function PopupmenuCallback(obj, hObject, ~, ~)
            % Item selected from the pop-up menu

            if ~obj.IsCurrentlyRunning()
                obj.Callback(get(hObject, 'Value'));
            end
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function my_save_file(cache_file, data_obj, force_overide, be_quite)

if nargin<3
    force_overide=false;
end

if nargin<4
    be_quite=false;
end

finished=false;
while ~finished
    try

        if ~force_overide && exist(cache_file, 'file')
            disp('file exist, save to new name.');

            [f1, f2, f3]=fileparts(cache_file);
            while exist(cache_file, 'file')
                f2=[f2 '_new'];
                cache_file=[f1 '/' f2 f3];
            end

        end

        tmp_dir=fileparts(cache_file);
        mkdir_notexist(tmp_dir);

        if ~be_quite
            fprintf('save_file:%s\n',cache_file);
        end

        save(cache_file, 'data_obj', '-v7.3');
        finished=true;
    catch err_info
        disp(err_info);
        pause_sec=60+60*rand(1);
        fprintf('save_file failed, retry in %.f seconds, file:%s\n', pause_sec, cache_file);
        pause(pause_sec);
    end
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function mtexColorMap(arg1,varargin)
% define an MTEX colormap

% get input
if length(arg1)==1 && ishandle(arg1) &&...
    (isa(arg1,'axis') || isa(arg1,'matlab.graphics.axis.Axes') || isa(arg1,'matlab.ui.Figure'))
  ax = arg1;
  name = varargin{1};
else
  ax = gcf;
  name = arg1;
end

% get axes
try
  mtexFig = getappdata(ax,'mtexFig');
  ax = mtexFig.children;
end

% detect colormap
if ischar(name)
  if isempty(which([name '.m']))
    name = [name,'ColorMap'];
    if isempty(which(name))
      error('unknown colormap name');
    end
  end
  try
    map = feval(name);
  catch
    map = feval(lower(name));
  end
else
  map = name;
end

% apply the colormap
for i = 1:length(ax), colormap(ax(i),map); end

end

author: mtex-toolbox

repo: mtex

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


text:
function autogradfunc = autograd(problem, fixedrankflag)
% Apply automatic differentiation to computing the Euclidean gradient
%
% function autogradfunc = autograd(problem)
% function autogradfunc = autograd(problem, fixedrankflag)
%
% Returns an AcceleratedFunction or a function handle which can be used to
% compute Euclidean gradients. See https://ch.mathworks.com/help/
% deeplearning/ref/deep.acceleratedfunction.html for more descriptions
% about AcceleratedFunction.
%
% Note: to evaluate the Euclidean gradient of a certain point x(x should be
% of type dlarray), call dfeval(autogradfunc,x) instead of autogradfunc(x).
%
% See also: manoptAD, egradcompute, costgradcompute

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, Aug. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:
%
% To do: Add AD to fixedTTrankfactory, fixedranktensorembeddedfactory
% and the product manifold which contains fixedrankembeddedfactory
% or anchoredrotationsfactory

    % Check availability
    assert(isfield(problem,'M') && isfield(problem,'cost'),...
    'problem structure must contain the fields M and cost.');
    assert(exist('dlarray', 'file') == 2, ['Deep learning tool box is '...
    'needed for automatic differentiation'])

    % Set fixedrankflag to false if the manifold struct is not
    % fixed(multilinear)-rank matrices or tensors with an embedded geometry
    % or tensors of fixed Tensor Train (TT) rank
    if ~exist('fixedrankflag', 'var')|| isempty(fixedrankflag)
        fixedrankflag = false;
    end

    % Obtain the euclidean gradient function via AD
    costfunction = problem.cost;
    % Set fixedrankflag to true if the manifold is fixed-rank matrices with
    % an embedded geometry. The other two cases are not implemented yet.
    if fixedrankflag
        % AcceleratedFunction can lead to a slow down in this case
        autogradfunc = @(x,A,B) autogradfuncinternelfixedrankembedded(x,A,B);
    else
        func = @autogradfuncinternal;
        % accelerate
        try
            autogradfunc = dlaccelerate(func); % Introduced in Matlab 2021a
            clearCache(autogradfunc);
        catch
            warning('manopt:dlaccelerate', ...
                    ['Function dlaccelerate is not available:\nPlease ' ...
                     'upgrade to Matlab 2021a or later and the latest deep\nlearning ' ...
                     'toolbox version if possible.\nMeanwhile, auto-diff ' ...
                     'may be somewhat slower.\n The hessian is not available as well.\n' ...
                     'To disable this warning: warning(''off'', ''manopt:dlaccelerate'')']);
            autogradfunc = func;
        end
    end

    % define Euclidean gradient function
    function [y, egrad] = autogradfuncinternal(x)

        y = costfunction(x);
        % In case that the user forgot to take the real part of the cost
        % when dealing with complex problems with Matlab R2021a or earlier,
        % take the real part for AD
        if iscstruct(y)
            y = creal(y);
        end

        % Call dlgradient to compute the Euclidean gradient. by default,
        % 'RetainData' and 'EnableHigherDerivatives' are set to false
        egrad = dlgradient(y, x);

        % in case that the user is optimizing over anchoredrotationsfactory
        % egrad of anchors with indices in A should be zero
        problem_name = problem.M.name();
        if (contains(problem_name,'Product rotations manifold') &&...,
            contains(problem_name,'anchors') &&...,
            ~startsWith(problem_name,'Product manifold'))
            A = findA_anchors(problem);
            egrad(:, :, A) = 0;
        end
    end

    % fixedrankembeddedfactory part
    % obtain the product of egrad and V and the product of egrad
    % transpose and U by differentiating g1 and g2 w.r.t A and B
    function [g1, egrad] = autogradfuncinternelfixedrankembedded(x, A, B)
        X1.U = A; X1.S = eye(size(x.S,1)); X1.V = x.V;
        X2.U = x.U; X2.S = eye(size(x.S,1)); X2.V = B;
        g1 = costfunction(X1); g2 = costfunction(X2);
        egrad.A = dlgradient(g1,A);  egrad.B = dlgradient(g2,B);
    end

end

author: NicolasBoumal

repo: manopt


text:
function [cmap_i]=resampleColormap(cmap,n)

% function [cmap_i]=resampleColormap(cmap,n)
% ------------------------------------------------------------------------
%
% This function resamples the input colormap cmap using n steps. Resampling
% is based on linear interpolation.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log:
% 2014/09/25
% 2018/06/13 Added handling of single color colormap
%------------------------------------------------------------------------

%%

if n~=size(cmap,1) %If resampling is required
    if size(cmap,1)==1
        cmap_i=cmap(ones(n,1),:);
    else
        ind=(1:1:size(cmap,1))';
        ind_i=linspace(1,size(cmap,1),n)';
        cmap_i=zeros(n,size(cmap,2));

        %Interpolate color data
        for q=1:1:size(cmap,2)
            cmap_i(:,q)=interp1(ind,cmap(:,q),ind_i,'linear');
        end
    end
else
    cmap_i=cmap;
end

%%
% _*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 vars(vh)

num_masses = vh.userdata.num_masses;

lbx = horzcat(-4*ones(1,num_masses), -inf*ones(1,num_masses)).';
ubx = horzcat(4*ones(1,num_masses), inf*ones(1,num_masses)).';

vh.addState('x', [2*num_masses 1], 'lb', lbx, 'ub', ubx);
vh.addControl('u', [num_masses-1 1], 'lb', -0.5, 'ub', 0.5);

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function geometry_test0235 ( )

%*****************************************************************************80
%
%% TEST0235 tests DMS_TO_RADIANS, RADIANS_TO_DMS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0235\n' );
  fprintf ( 1, '  DMS_TO_RADIANS converts an angle from\n' );
  fprintf ( 1, '    degrees/minutes/seconds to radians;\n' );
  fprintf ( 1, '  RADIANS_TO_DEGREES converts an angle from radians\n' );
  fprintf ( 1, '    to degrees/minutes/seconds;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radians     DMS     Radians\n' );
  fprintf ( 1, '\n' );

  for i = -2 : 15
    angle_rad = pi * i / 7.0;
    [ angle_deg, angle_min, angle_sec ] = radians_to_dms ( angle_rad );
    angle_rad2 = dms_to_radians ( angle_deg, angle_min, angle_sec );
    fprintf ( 1, '  %10f  %4d  %3d  %3d  %10f\n', ...
      angle_rad, angle_deg, angle_min, angle_sec, angle_rad2 );
  end

  return
end


text:
function [x,fval,exitflag,output,optimState,gpstruct] = fixedbads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options,fixidx,nout)
%FIXEDBADS Run BADS with reduced dimensionality due to fixed variables.

nvars = size(x0,2);

% Initialize outputs
fval = [];
exitflag = [];
output = [];
optimState = [];
gpstruct = [];

freeidx = ~fixidx;
fixedvars = LB(fixidx);

% Update periodic variables, if present in OPTIONS
if isfield(options,'PeriodicVars') && ~isempty(eval(options.PeriodicVars))
    pervars = false(1,nvars);
    pervars(eval(options.PeriodicVars)) = true;
    pervars(fixidx) = [];
    options.PeriodicVars = find(pervars);
end

% Fix input FUNVALUES, if present
if isfield(options,'FunValues') && ~isempty(options.FunValues) && ...
        isfield(options.FunValues,'X') && ~isempty(options.FunValues.X)
    options.FunValues.X = options.FunValues.X(:,freeidx);
end

% Downlift
x0 = x0(freeidx);
LB = LB(freeidx);
UB = UB(freeidx);
PLB = PLB(freeidx);
PUB = PUB(freeidx);

switch nout
    case 1
        x_free = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
    case 2
        [x_free,fval] = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
    case 3
        [x_free,fval,exitflag] = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
    case 4
        [x_free,fval,exitflag,output] = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
    case 5
        [x_free,fval,exitflag,output,optimState] = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
    case 6
        [x_free,fval,exitflag,output,optimState,gpstruct] = bads(fun_fix,x0,LB,UB,PLB,PUB,nonbcon_fix,options);
end

% Uplift solution
x = zeros(1,nvars);
x(freeidx) = x_free;
x(fixidx) = fixedvars;

if nout > 4
    % Uplift points in OPTIMSTATE
    nx = size(optimState.X,1);
    Xtemp = zeros(nx,nvars);
    Xtemp(:,freeidx) = optimState.X;
    Xtemp(:,fixidx) = repmat(fixedvars, [nx 1]);
    optimState.X = Xtemp;
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
% Given a set of nonnegative vectors x{i}, i=1,2,...,r
%
% Compute the function g(mu) = sum_i beta_i [ e^T x_i(mu) ]
% where
% beta_i = 1/[ ||w{i}|| - min_j w{i}(j) ], and
% If x{i} - beta_i mu w{i} contains at least one postive entry:
%      x_i(mu) = (x{i} - beta_i mu w{i})_+;
%      x_i(mu) = x_i(mu)/||x_i(mu)||_2;
% Otherwise
%     x_i(mu) is 1-sparse with nonzero entry equal to one at position
%     corresponding to the largest entry of x{i} - beta_i mu w{i}

function [vgmu,xp,gradg] = wgmu(x,w,mu);

vgmu = 0;
gradg = 0;
for i = 1 : length(x)
    ni = length(x{i});
    betai = 1/(norm(w{i})-min(w{i}));
    xp{i} = x{i} - mu*betai*w{i};
    indtp = find(xp{i} > 0);
    % Gradient g(mu)
    if ~isempty(indtp)
        xp{i} = max(0,xp{i});
        f2 = norm(xp{i});
        if nargout >= 3
            nip = w{i}(indtp)'*w{i}(indtp);

            gradg = gradg + betai^2 * ( - nip * f2^(-1) ...
                + (w{i}(indtp)'*xp{i}(indtp))^2 * f2^(-3)  );
        end
        xp{i} = xp{i}/norm(xp{i},2);
        vgmu = vgmu + betai*sum(xp{i}.*w{i});
    else
        [~,im] = max(xp{i});
        xp{i} = zeros( ni , 1 );
        xp{i}(im) = 1;
        vgmu = vgmu + betai*w{i}(im);
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
clear, clc;

% This is an example for running the function LogisticR
%
%  Problem:
%
%  min  f(x,c) = - weight_i * log (p_i) + 1/2 * rsL2 * ||x||_2^2
%                 + rho * \|x\|_1
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
% [3]  Jun Liu, Jianhui Chen, and Jieping Ye,
%      Large-Scale Sparse Logistic Regression, KDD, 2009.
%
%% ------------   History --------------------
%
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

m=1000;  n=1000;     % The data matrix is of size m x n

randNum=2;           % a random number

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,1);

randn('state',(randNum-1)*3+3);
y=[ones(n/2,1);...
    -ones(n/2, 1)];  % the response

rho=0.1;             % the regularization parameter
                     % it is a ratio between (0,1), if .rFlag=1

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;    % maximum number of iterations

% Normalization
opts.nFlag=0;       % without normalization

% Regularization
opts.rFlag=1;       % the input parameter 'rho' is a ratio in (0, 1)
%opts.rsL2=0.01;     % the squared two norm term

% Group Property
opts.sWeight=[1,1]; % set the weight for positive and negative samples

%----------------------- Run the code LeastR -----------------------
fprintf('\n mFlag=0, lFlag=0 \n');
opts.mFlag=0;       % treating it as compositive function
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, c1, funVal1, ValueL1]= LogisticR(A, y, rho, opts);
toc;

opts.maxIter=1000;

fprintf('\n mFlag=1, lFlag=0 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=0;       % Nemirovski's line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, c2, funVal2, ValueL2]= LogisticR(A, y, rho, opts);
toc;

fprintf('\n mFlag=1, lFlag=1 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x3, c3, funVal3, ValueL3]= LogisticR(A, y, rho, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
hold on;
plot(funVal3,':g');
legend('mFlag=0, lFlag=0', 'mFlag=1, lFlag=0', 'mFlag=1, lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% --------------------- compute the pathwise solutions ----------------
opts.fName='LogisticR';    % set the function name to 'LogisticR'
Z=[0.5, 0.2, 0.1, 0.01];   % set the parameters

% run the function pathSolutionLogistic
fprintf('\n Compute the pathwise solutions, please wait...');
[X,C]=pathSolutionLogistic(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
function gui_close(hObject, eventdata, handles)

end

author: andrewssobral

repo: lrslibrary


text:
function varargout = kdtreebuild(varargin)
% VL_KDTREEBUILD  Build randomized kd-tree
%   FOREST = VL_KDTREEBUILD(X) returns a structure FOREST containing
%   the kd-tree indexing the data X. X is a M x N dimensional matrix
%   of class DOUBLE or SINGLE with one data point per column. Use
%   VL_KDTREEQUERY() to query the index.
%
%   VL_KDETREEBUILD() accepts the following options:
%
%   NumTrees:: [1]
%     Specifies the number of semi-randomized KD-trees to add to the
%     forest. Multiple trees can be used to improve approximate
%     nearest-neighbors searches with VL_KDTREEQUERY().
%
%   ThresholdMethod:: [MEDIAN]
%     Specifies the thresholding method used to split the data.  Two
%     methods are supported: split around the median (MEDIAN) or split
%     around the mean (MEAN).
%
%   Verbose::
%     Increases the verbosity level (can be repeated).
%
%   The FOREST structure has the following fields:
%
%   FOREST.TREES::
%     A structure array with one element per tree.
%
%   FOREST.DIMENSION::
%     Dimensionality of the indexed data.
%
%   FOREST.NUMDATA::
%     Number of indexed data points.
%
%   The TREE structure has the following fields:
%
%   TREE.NODES::
%      A structure array representing the nodes of the tree.
%
%   TREE.DATAINDEX::
%      A 1 x NUMDATA vector of class UINT32 representing a permutation
%      of the data.
%
%   Nodes are numbered from 1 to NUMNODES.  The NODES structure array
%   has the following fields:
%
%   NODES.LOWERCHILD and NODES.UPPERCHILD::
%      1 x NUMNODES vectors of class INT32. A positive value is the
%      index of the lower/upper child node. A negative value denotes a
%      leaf and is (after negation) is the first or last element plus
%      one of a range of entries in the permutation TREE.DATAINDEX.
%      Such entries are in turn indexes of the data points that belong
%      to that leaf. Typically there is one point per leaf.
%
%   NODES.SPLITDIMENSION and NODES.SPLITTHRESHOLD::
%      1 x NUMNODES vector of class UINT32 and DOUBLE, respectively,
%      with the index of the splitting dimension and the threshold for
%      each node.
%
%   See also: VL_KDTREEQUERY(), VL_HELP().
[varargout{1:nargout}] = vl_kdtreebuild(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
function trace = r8mat_trace ( n, a )

%*****************************************************************************80
%
%% R8MAT_TRACE computes the trace of an R8MAT.
%
%  Discussion:
%
%    The trace of a square matrix is the sum of the diagonal elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix A.
%
%    Input, real A(N,N), the matrix whose trace is desired.
%
%    Output, real TRACE, the trace of the matrix.
%
  trace = 0.0;
  for i = 1 : n
    trace = trace + a(i,i);
  end

  return
end


text:
function V = solveV_eig(U, Us, Vs, Ss, reg_smooth, reg_l2)
    % solve S analytically using eigen-decomposition
    %
    % TODO: check with optimization form.

    n = length(Ss);

    UtU = U' * U;

    [Q1, Lu] = eig(UtU);

    V = cell(n, 1);

    for ii = 1:n
        V{ii} = solveVi(U, Us, Q1, diag(Lu) + reg_l2, Vs{ii}, Ss{ii}, reg_smooth);
    end

end

function Vi = solveVi (U, Us, Q1, Lu, Vsi, Ssi,  reg_smooth)

ti = size(Ssi, 2);

iy  = [1: (ti - 1), 1: (ti - 1)];
ix  = [1: (ti - 1), 2: ti];
val = [ones(1, ti-1), -1 * ones(1, ti-1)];
Ri  = sparse(ix, iy, val, ti, ti-1);

%[Q2, Lr] = eig(reg_smooth * (Ri * Ri'), ti);
[Q2, Lr] = eig(reg_smooth * full(Ri * Ri'));
Lr = diag(Lr);

UtSi = Q1' * (U' * Us) * Vsi * Q2 + Q1' * (U' * Ssi) * Q2;

V_hat = zeros(size(UtSi));

for i = 1:size(V_hat, 1)
    for j = 1: size(V_hat, 2)
        V_hat(i, j) = UtSi(i, j) / (Lu(i) + Lr(j));
    end
end

Vi = Q1 * V_hat * Q2';

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function imdb=my_gen_imdb(train_opts, ds_info)

task_info_train=[];
task_info_train.data_type='train';
task_info_train.task_idxes=ds_info.train_idxes;

task_info_val=[];
task_info_val.data_type='val';
task_info_val.task_idxes=ds_info.test_idxes;

imdb=[];
imdb.ds_info=ds_info;

imdb.task_info_train=task_info_train;
imdb.task_info_val=task_info_val;

task_idxes=[];
task_idxes=cat(1, task_idxes, task_info_train.task_idxes);
task_idxes=cat(1, task_idxes, task_info_val.task_idxes);

imdb.max_task_idx=max(task_idxes);
imdb.task_idxes=task_idxes;

make_ref_obj(imdb);

end


text:
function X = readbmpdata(info)
%READBMPDATA Read bitmap data
%   X = readbmpdata(INFO) reads image data from a BMP file.  INFO is a
%   structure returned by IMBMPINFO. X is a uint8 array that is 2-D for
%   1-bit, 4-bit, and 8-bit image data.  X is M-by-N-by-3 for 24-bit and
%   32-bit image data.

%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.1.6.5 $  $Date: 2011/05/13 17:34:37 $

offset = info.ImageDataOffset;
width = info.Width;
height = info.Height;
filename = info.Filename;

switch info.CompressionType

 case 'none'

  switch info.BitDepth
   case 1
    X = logical(bmpReadData1(filename, offset, width, height));

   case 4
    X = bmpReadData4(filename, offset, width, height);

   case 8
    X = bmpReadData8(filename, offset, width, height);

   case 16
    X = bmpReadData16(filename, offset, width, height);

   case 24
    X = bmpReadData24(filename, offset, width, height);

   case 32
    X = bmpReadData32(filename, offset, width, height);

  end

 case '8-bit RLE'
  X = bmpReadData8RLE(filename, offset, width, height);

 case '4-bit RLE'
  X = bmpReadData4RLE(filename, offset, width, height);

 case 'bitfields'
  error(message('MATLAB:audiovideo:readbmpdata:bitfieldUnsupported'));

 case 'Huffman 1D'
  error(message('MATLAB:audiovideo:readbmpdata:huffmanUnsupported'));

 case '24-bit RLE'
  error(message('MATLAB:audiovideo:readbmpdata:rleUnsupported'));

end

%%%
%%% bmpReadData8 --- read 8-bit bitmap data
%%%
function X = bmpReadData8(fname, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
paddedWidth = 4*ceil(width/4);

fid = fopen(fname,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
X = fread(fid,paddedWidth*abs(height),'*uint8');
status = fclose(fid);

count = length(X);
if (count ~= paddedWidth*abs(height))
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(paddedWidth*abs(height)) = 0;
end

if height>=0
  X = rot90(reshape(X, paddedWidth, height));
else
  X = [reshape(X, paddedWidth, abs(height))]';
end

if (paddedWidth ~= width)
    X = X(:,1:width);
end

%%%
%%% bmpReadData8RLE --- read 8-bit RLE-compressed bitmap data
%%%
function X = bmpReadData8RLE(fname, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
paddedWidth = 4*ceil(width/4);

fid = fopen(fname,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
inBuffer = fread(fid,'*uint8');
status = fclose(fid);

X = bmpdrle(inBuffer, paddedWidth, abs(height), 'rle8');

if height>=0
  X = rot90(X);
else
  X = X';
end
if (paddedWidth ~= width)
    X = X(:,1:width);
end

%%%
%%% bmpReadData4 --- read 4-bit bitmap data
%%%
function X = bmpReadData4(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
paddedWidth = 8*ceil(width/8);
numBytes = paddedWidth * abs(height) / 2; % evenly divides because of padding

fid = fopen(filename,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
XX = fread(fid,numBytes,'*uint8');
status = fclose(fid);

count = length(XX);
if (count ~= numBytes)
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(numBytes) = 0;
end
XX = reshape(XX, paddedWidth / 2, abs(height));

X = repmat(uint8(0), paddedWidth, abs(height));
X(1:2:end,:) = bitslice(XX,5,8);
X(2:2:end,:) = bitslice(XX,1,4);

if height>=0
  X = rot90(X);
else
  X = X';
end
if (paddedWidth ~= width)
    X = X(:,1:width);
end

%%%
%%% bmpReadData4RLE --- read 4-bit RLE-compressed bitmap data
%%%
function X = bmpReadData4RLE(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
paddedWidth = 8*ceil(width/8);
numBytes = paddedWidth * abs(height) / 2; % evenly divides because of padding

fid = fopen(filename,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
inBuffer = fread(fid,'*uint8');
status = fclose(fid);

if height>=0
  X = rot90(bmpdrle(inBuffer, paddedWidth, abs(height), 'rle4'));
else
  X = [bmpdrle(inBuffer, paddedWidth, abs(height), 'rle4')]';
end
if (paddedWidth ~= width)
    X = X(:,1:width);
end

%%%
%%% bmpReadData1 --- read 1-bit bitmap data
%%%
function X = bmpReadData1(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
paddedWidth = 32*ceil(width/32);
numPixels = paddedWidth * abs(height);  % evenly divides because of padding

% 1-bit BMP data has big-endian byte ordering
fid = fopen(filename,'r','ieee-be');
status = fseek(fid,offset,'bof');
if status==-1
  status = fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
[X, count] = fread(fid,paddedWidth*abs(height),'*ubit1');
status = fclose(fid);

if (count ~= numPixels)
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(numPixels) = 0;
end
X = reshape(X, paddedWidth, abs(height));

if height>=0
  X = rot90(X);
else
  X = X';
end

if (paddedWidth ~= width)
    X = X(:,1:width);
end

%%%
%%% bmpReadData16 --- read 16-bit bitmap data
%%%
function RGB = bmpReadData16(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
byteWidth = width;
paddedByteWidth = 4*ceil(byteWidth/4);
numBytes = paddedByteWidth * abs(height);

fid = fopen(filename,'r','ieee-le');
status = fseek(fid,offset,0);
%Check if seek is ok
X = fread(fid,numBytes,'*uint16');
status = fclose(fid);

count = length(X);
if (count ~= numBytes)
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(numBytes) = 0;
end

if height>=0
  X = rot90(reshape(X, paddedByteWidth, abs(height)));
else
  X = [reshape(X, paddedByteWidth, abs(height))]';
end

if (paddedByteWidth ~= byteWidth)
    X = X(:,1:byteWidth);
end

RGB(1:abs(height), 1:width, 1) = uint8(bitslice(X,11,15));
RGB(:,:,2) = uint8(bitslice(X,6,10));
RGB(:,:,3) = uint8(bitslice(X,1,5));

%Scale data for display
RGB = bitor(bitshift(RGB,3),bitshift(RGB,-2));

%%%
%%% bmpReadData24 --- read 24-bit bitmap data
%%%
function RGB = bmpReadData24(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
byteWidth = 3*width;
paddedByteWidth = 4*ceil(byteWidth/4);
numBytes = paddedByteWidth * abs(height);

fid = fopen(filename,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
X = fread(fid,numBytes,'*uint8');
status = fclose(fid);

count = length(X);
if (count ~= numBytes)
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(numBytes) = 0;
end

if height>=0
  X = rot90(reshape(X, paddedByteWidth, abs(height)));
else
  X = [reshape(X, paddedByteWidth, abs(height))]';
end

if (paddedByteWidth ~= byteWidth)
    X = X(:,1:byteWidth);
end

RGB(1:abs(height), 1:width, 3) = X(:,1:3:end);
RGB(:, :, 2) = X(:,2:3:end);
RGB(:, :, 1) = X(:,3:3:end);

%%%
%%% bmpReadData32 --- read 32-bit bitmap data
%%%
function RGB = bmpReadData32(filename, offset, width, height)

% NOTE: BMP files are stored so that scanlines use a multiple of 4 bytes.
byteWidth = 4*width;
paddedByteWidth = 4*ceil(byteWidth/4);
numBytes = paddedByteWidth * abs(height);

fid = fopen(filename,'r','ieee-le');
status = fseek(fid,offset,'bof');
if status==-1
  fclose(fid);
  error(message('MATLAB:audiovideo:readbmpdata:invalidDataOffset'));
end
X = fread(fid,numBytes,'*uint8');
status = fclose(fid);

count = length(X);
if (count ~= numBytes)
    warning(message('MATLAB:audiovideo:readbmpdata:truncatedImageData'));
    % Fill in the missing values with zeros.
    X(numBytes) = 0;
end

if height>=0
  X = rot90(reshape(X, paddedByteWidth, abs(height)));
else
  X = [reshape(X, paddedByteWidth, abs(height))]';
end

if (paddedByteWidth ~= byteWidth)
    X = X(:,1:byteWidth);
end

RGB(1:abs(height), 1:width, 3) = X(:,1:4:end);
RGB(:, :, 2) = X(:,2:4:end);
RGB(:, :, 1) = X(:,3:4:end);

author: kristinbranson

repo: JAABA


text:
function [F,stdF,results] = nips2011_experiment_mohsst5(model, testset, seed)

% Seed for the random number generator
if nargin < 3
  seed = 10;
end
fprintf('Seed: %g\n', seed);
rand('state', seed);
randn('state', seed);

% Load the data
data = mohsst5_loaddata();
sea = sum(~isnan(data.observations),2) > 0;
Y = data.observations(sea,:);
[M,N] = size(Y);

switch testset
 case 1 % randomly (uniformly) selected
  disp('Test set: Uniform');

  I = rand(size(Y)) < 0.2;
  Itest = I & ~isnan(Y);
  Itrain = ~I & ~isnan(Y);
  Ytest = Y;
  Ytrain = Y;
  Ytest(~Itest) = nan;
  Ytrain(~Itrain) = nan;

  test_percent = sum(Itest(:)) / sum(~isnan(Y(:)));
  fprintf('Size of the test set: %.1f percent\n', 100*test_percent);

% $$$   figure
% $$$   spy(Itest)
% $$$   figure
% $$$   spy(Itrain)
% $$$   figure
% $$$   spy(~isnan(Y));
% $$$   figure
% $$$   spy(Itrain&Itest)
% $$$   figure
% $$$   spy((Itrain|Itest)==~isnan(Y));
% $$$   return

  testset_string = 'uniform';

 case 2 % pattern from earliest years used for latest years
  disp('Test set: Pattern');

  Imv = isnan(Y);
  Imv(:,(20*12+1):end) = false;
  Itest = Imv(:,end:-1:1) & ~isnan(Y);
  Itrain = ~Imv(:,end:-1:1) & ~isnan(Y);

  Ytest = Y;
  Ytrain = Y;
  Ytest(~Itest) = nan;
  Ytrain(~Itrain) = nan;

  test_percent = sum(Itest(:)) / sum(~isnan(Y(:)));
  fprintf('Size of the test set: %.1f percent\n', 100*test_percent);

% $$$   figure
% $$$   spy(Itest)
% $$$   figure
% $$$   spy(Itrain)
% $$$   figure
% $$$   spy(~isnan(Y));
% $$$   figure
% $$$   spy(Itrain&Itest)
% $$$   figure
% $$$   spy((Itrain|Itest)==~isnan(Y));
% $$$   return

  testset_string = 'pattern';

 otherwise
  error('Unknown test set requested');
end

switch model

 case 1

  %
  % Local GP
  %

  %disp('Model: Local GP');
  disp('Model: Local GP with two covfuncs per domain');

  % Temporal covariance function
  d = abs(1-(1:length(data.time)));
% $$$   covfunc1 = gp_cov_toeplitz(gp_cov_pp(d,1));
  covfunc1 = gp_cov_toeplitz(gp_cov_sum(gp_cov_pp(d,1), ...
                                        gp_cov_scale(gp_cov_pp(d,1))));
% $$$   theta_temporal = [7];   % length scale
  theta_temporal = [7;   % length scale 1
                    1.0; % magnitude 2
                    3];  % length scale 2

  % Spatial covariance function
  [LON,LAT] = meshgrid(data.longitude,data.latitude);
  X = geographic_to_euclidean([LON(:)';LAT(:)']);

  % Use block-Toeplitz structure for the covariance function
  [lat,lon0] = meshgrid(data.latitude,data.longitude(1));
  X0 = geographic_to_euclidean([lon0(:)';lat(:)']);
  D = sqrt(sq_dist(X0,X));
% $$$   covfunc2 = gp_cov_toeplitz_block(gp_cov_pp(D,3));
  covfunc2 = gp_cov_toeplitz_block(gp_cov_sum(gp_cov_pp(D,3), ...
                                              gp_cov_scale(gp_cov_pp(D,3))));

  % Select sea areas
  covfunc2 = gp_cov_select(covfunc2, sea);

% $$$   theta_spatial = [3000];  % length scale
  theta_spatial = [3000;  % length scale 1
                   1.0;   % magnitude 2
                   2000]; % length scale 2

  % Initial guess for covariance parameters
  theta_init = [0.5; ...               % total magnitude
                theta_temporal(:); ... % temporal parameters
                theta_spatial(:); ...  % spatial parameters
                0.5]';                 % noise magnitude

  % Prior for the hyperparameters
  a = 1e-3;
  b = 1e-3;
  logprior_theta = @(theta) sum(gamma_logpdf(theta, a, b));
  dlogprior_theta = @(theta) gamma_dlogpdf(theta, a, b);

  % Noise scaling
  w = 1./sqrt(cosd(LAT));
  noise_scale = repmat(w(sea), [1,size(Y,2)]);

  % Inference
  N_samples = 1000;
  burnin = floor(N_samples/2);
  res = gp_kron(Ytrain, covfunc1, covfunc2, N_samples, theta_init, ...
                logprior_theta, dlogprior_theta, burnin, 'noise_scale', ...
                noise_scale);

  F = res.F;
  stdF = sqrt(res.FF - res.F.^2);
  results = res;

  %model_string = 'gpkron';
  model_string = 'gpkron-cov2';

 case 2

  %
  % VB PCA
  %

  disp('Model: VB PCA');

  D = 150;

  % PCA module for X (one constant component for modeling bias)
  prior.mu = [1; zeros(D-1,1)];
  prior.CovX = diag([1e-6; ones(D-1,1)]);
  X_module = factor_module_iid(D, N, 'prior', prior);

  % ARD module for W
  W_module = factor_module_ard(D, M);

  % Isotropic noise (precisions weighted proportionally to grid size)
  [LON,LAT] = meshgrid(data.longitude, data.latitude);
  weights = cosd(LAT(:));
  weights = weights(sea); %metoffice_remove_bins(weights,maskfile);
  weights = repmat(weights, [1, N]);
  noise_module = noise_module_isotropic(M, N, 1e-3, 1e-3, ...
                                        'init', 10, ...
                                        'weights', weights);

  % Run VB PCA
  Q = vbfa(D, Ytrain, W_module, X_module, noise_module, ...
           'maxiter', 500, ...
           'rotate', true);
% $$$            'autosavefile', filename, ...
% $$$            'autosave', [1 20:20:2000]);

  % Results
  F = Q.W'*Q.X;
  varF = reshape(Q.CovW, [D*D, M])' * reshape(Q.CovX, [D*D,N]);
  for m=1:M
    for n=1:N
      varF(m,n) = varF(m,n) + ...
          Q.W(:,m)'*Q.CovX(:,:,n)*Q.W(:,m) + ...
          Q.X(:,n)'*Q.CovW(:,:,m)*Q.X(:,n);
    end
  end
  stdF = sqrt(varF);
  results = Q;

  model_string = 'vbpca';

 otherwise
  error('Unknown model requested');
end

% Error measures
[LON,LAT] = meshgrid(data.longitude, data.latitude);
weights = cosd(LAT(:));
weights = weights(sea);
weights = repmat(weights, [1, N]);
rmse_train = rmsew(Y(Itrain)-F(Itrain),weights(Itrain));
rmse_test = rmsew(Y(Itest)-F(Itest),weights(Itest));
rmse_zero = rmsew(Y(Itest),weights(Itest));
fprintf('Training WRMSE=%.4f and testing WRMSE=%.4f\n', rmse_train, rmse_test);
fprintf('Testing WRMSE=%.4f with zero predictions\n', rmse_zero);

filename = sprintf(['/home/jluttine/matlab/publications/nips2011/' ...
                    'results_nips2011_mohsst5_%s_%s_%s'], model_string, ...
                   testset_string, datestr(now,'yyyymmdd'));

save(filename, 'F', 'stdF', 'results', 'Ytest', 'Ytrain');
fprintf('Saved results to %s\n', filename);

if nargout < 1
  clear F
  clear stdF
  clear results
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function test_suite = test_createSoccerBall
%TESTCREATESoccerBall  One-line description here, please.
%
%   output = testCreateSoccerBall(input)
%
%   Example
%   testCreateSoccerBall
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-12-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testCreation(testCase) %#ok<*DEFNU>

createSoccerBall();

function testVEFCreation(testCase)

[v, e, f] = createSoccerBall();
testCase.assertTrue(~isempty(v));
testCase.assertTrue(~isempty(e));
testCase.assertTrue(~isempty(f));

[nv, ne, nf] = getMeshElementsNumber;
testCase.assertEqual([nv 3], size(v));
testCase.assertEqual([ne 2], size(e));
testCase.assertEqual(nf, length(f));

function testVFCreation(testCase)

[v, f] = createSoccerBall();
testCase.assertTrue(~isempty(v));
testCase.assertTrue(~isempty(f));

[nv, ne, nf] = getMeshElementsNumber; %#ok<ASGLU>
testCase.assertEqual([nv 3], size(v));
testCase.assertEqual(nf, length(f));

function testMeshCreation(testCase)

mesh = createSoccerBall();
testCase.assertTrue(isstruct(mesh));
testCase.assertTrue(isfield(mesh, 'vertices'));
testCase.assertTrue(isfield(mesh, 'edges'));
testCase.assertTrue(isfield(mesh, 'faces'));

[nv, ne, nf] = getMeshElementsNumber;
testCase.assertEqual([nv 3], size(mesh.vertices));
testCase.assertEqual([ne 2], size(mesh.edges));
testCase.assertEqual(nf, length(mesh.faces));

function [nv, ne, nf] = getMeshElementsNumber

nv = 60;
ne = 90;
nf = 32;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [x, objV] = proxF_l1(x,t)
%% Usage proxF_l1(x,t)
% ProxF_l1(x,t) applies L1 shrinkage to n-dimensional matrix x for
% a shrinkage factor of t.  Returns the shrunk object.
%
%
% INPUTS:
%
% x: n-dimensional matrix
% t: shrinkage threshold
%
% OUTPUTS:
%
% x: n-dimensional matrix after shrinkage.
%
% Authors: G. Ely, S. Aeron, Z. Zhang, ECE, Tufts Univ. 03/16/2015

tq = t * 1;
s  = 1 - min( tq./abs(x), 1 );
x  = x .* s;
objV = sum(x(:));
end


text:
% Helper method for quickly creating bounding boxes.
%
% Args:
%   basis2d - 2x2 matrix for the basis in the XY plane
%   centroid - 1x3 vector for the 3D centroid of the bounding box.
%   coeffs - 1x3 vector for the radii in each dimension (x, y, and z)
%
% Returns:
%   bb - a bounding box struct.
%
% Author: Nathan Silberman (silberman@cs.nyu.edu)
function bb = create_bounding_box_3d(basis2d, centroid, coeffs)
  assert(all(size(basis2d) == [2, 2]));
  assert(numel(centroid) == 3);
  assert(numel(coeffs) == 3);

  centroid = centroid(:)';
  coeffs = coeffs(:)';

  bb = struct();
  bb.basis = zeros(3,3);
  bb.basis(3,:) = [0 0 1];
  bb.basis(1:2,1:2) = basis2d;

  bb.centroid = centroid;
  bb.coeffs = coeffs;
%   bb.volume = prod(2 * bb.coeffs);
end

author: thusiyuan


text:
function [ti,xi] = keepdata(x,p)
%
% Keeps arbitrary fraction p of equidistant data x at times ti.
% Output contains data xi at time ti

if p > 1, error('p must be between 0 and 1'), end

N=length(x);
Nover=floor(p.*N);

% over is the array of remaining points with value 1
over=zeros(1,N);

ti=zeros(1,Nover);
xi=zeros(1,Nover);

% Select at random Nover points
i = 0;
while (i < Nover)
   r = rand(1);
   k = round((N - 1) .* r) + 1;
   if (over(k) ~= 1)
      over(k) = 1;
      i = i + 1;
   end
end

% Generate output data
j = 1;
for i=1:N
   if (over(i) == 1)
      xi(j) = x(i);
      ti(j) = i;
      j = j + 1;
   end
end
fprintf('Over: %d\n', j-1);

author: Sable

repo: mcbench-benchmarks