text:
function varargout = panel_opticalflow( varargin )
% PANEL_OPTICALFLOW:

% @=============================================================================
% 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: Syed Ashrafulla, 2010-2013

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(TimeVector) %#ok<DEFNU>
    % CREATE_PANEL Setup Java panel for calculating optical flow
    % INPUTS:
    %   TimeVector    - Vector of times in reconstructed map, used to
    %                   automatically correct window chosen by user
    % OUTPUTS:
    %   bstPanelNew   - handle to panel figure
    %   panelName     - name of panel for registration into BrainStorm
    panelName = 'OpticalFlow';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    % Constants
    DEFAULT_HEIGHT = 20;
    TEXT_WIDTH     = 50;

    % Create tool panel
    jPanelNew = gui_river([4,4], [14,5,0,10]);

    % ==== PANEL: SETUP (Horn-Schunck + HHD) ====
    jPanelSetup = gui_river([1,2], [5,15,15,10], 'Setup');

    % Horn-Schunck parameter (label)
    jLabelHornSchunck = JLabel('Horn-Schunck: ');
    jPanelSetup.add('br', jLabelHornSchunck);
    % Horn-Schunck parameter (Text)
    jTextHornSchunck = JTextField(num2str(0.01));
    jTextHornSchunck.setPreferredSize(Dimension(TEXT_WIDTH, DEFAULT_HEIGHT));
    jTextHornSchunck.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
    jPanelSetup.add('tab', jTextHornSchunck);

    % Calculate HHD as well as optical flow
    jCheckHHD = JCheckBox('Calculate HHD', 0);
    jCheckHHD.setEnabled(false); % ======>>>>> RIGHT NOW, HHD IS NOT ENABLED
    java_setcb(jCheckHHD, 'ActionPerformedCallback', @(h,ev)UpdatePanel());
    jPanelSetup.add('br', jCheckHHD);

    % Recursive Depth of HHD (label)
    jLabelRecursiveDepth = JLabel('Recursive Depth: ');
    jLabelRecursiveDepth.setEnabled(false); % Initially, no HHD
    jPanelSetup.add('br', jLabelRecursiveDepth);
    % Recursive Depth of HHD (Text)
    jTextRecursiveDepth = JTextField(num2str(2));
    jTextRecursiveDepth.setPreferredSize(Dimension(TEXT_WIDTH, DEFAULT_HEIGHT));
    jTextRecursiveDepth.setHorizontalAlignment(javax.swing.JLabel.RIGHT);
    jTextRecursiveDepth.setEnabled(false); % Initially, no HHD
    jPanelSetup.add('tab', jTextRecursiveDepth);

    jPanelNew.add('br hfill', jPanelSetup); % Add to main panel (jPanelNew) //

    % ==== PANEL: OPTIONS ====
    jPanelOptions = gui_river([1,2], [5,15,15,10], 'Setup');

    % Show optical flow after calculating
    jCheckVisualize = JCheckBox('Show results', 1);
    jPanelOptions.add('br hfill', jCheckVisualize);

    % Rotate normals of vertex to radius, so flows may fall in line better
    jCheckRotate = JCheckBox('Include inflated-sphere results', 0);
    jPanelOptions.add('br hfill', jCheckRotate);

    jPanelNew.add('br hfill', jPanelOptions); % Add to main panel (jPanelNew) //

    % ==== PANEL: OPTIONS ====
    jPanelStates = gui_river([1,2], [5,15,15,10], 'Flow States');

    % Include microstate results
    jCheckStates = JCheckBox('Calculate stable/transition states', 0);
    java_setcb(jCheckStates, 'ActionPerformedCallback', @(h,ev)UpdatePanel());
    jPanelStates.add('br hfill', jCheckStates);

    % Publish microstate results
    jCheckPublish = JCheckBox('Print stable/transition states', 0);
    jCheckPublish.setEnabled(false);
    jPanelStates.add('br hfill', jCheckPublish);

    % Only calculate microstates
    jCheckStatesOnly = JCheckBox('Only calculate microstates', 0);
    java_setcb(jCheckStatesOnly, 'ActionPerformedCallback', @(h,ev)UpdatePanel());
    jPanelStates.add('br hfill', jCheckStatesOnly);

    jPanelNew.add('br hfill', jPanelStates); % Add to main panel (jPanelNew) //

    % ==== PANEL: TIME INTERVAL ====
    jPanelTime = gui_river([1,1], [0,6,6,6], 'Interval');
    jPanelTime.add('p', JLabel('Start: '));
    jTextTimeStart = JTextField(''); % Flow: Time START
    jTextTimeStart.setPreferredSize(Dimension(TEXT_WIDTH, DEFAULT_HEIGHT));
    jTextTimeStart.setHorizontalAlignment(JTextField.RIGHT);
    jPanelTime.add('tab', jTextTimeStart);

    jPanelTime.add(JLabel('    End: '));
    jTextTimeStop = JTextField(''); % Flow: Time STOP
    jTextTimeStop.setPreferredSize(Dimension(TEXT_WIDTH, DEFAULT_HEIGHT));
    jTextTimeStop.setHorizontalAlignment(JTextField.RIGHT);
    jPanelTime.add(jTextTimeStop);

    % Set time controls callbacks
    TimeUnit = gui_validate_text(jTextTimeStart, [], jTextTimeStop, TimeVector, 'time', [], TimeVector(1),   []);
    TimeUnit = gui_validate_text(jTextTimeStop, jTextTimeStart, [], TimeVector, 'time', [], TimeVector(end), []);
    jLabelTimeUnit = JLabel(TimeUnit); % Display time unit
    jPanelTime.add(jLabelTimeUnit);

    jPanelNew.add('br hfill', jPanelTime); % Add to main panel (jPanelNew) //

    % Separator
    jPanelNew.add('br', JLabel(' '));

    % ===== VALIDATION BUTTONS =====
    jButtonCancel = JButton('Cancel'); % Cancel button
    java_setcb(jButtonCancel, 'ActionPerformedCallback', @ButtonCancel_Callback);
    jPanelNew.add('right', jButtonCancel);

    jButtonRun = JButton('Run'); % Run button
    java_setcb(jButtonRun, 'ActionPerformedCallback', @ButtonRun_Callback);
    jPanelNew.add(jButtonRun);

    jPanelNew.add('br',JLabel(' '));  % Add to main panel (jPanelNew) //

    % Insert panel in a scroll panel
    jScrollPanelNew = JScrollPane(jPanelNew);
    jScrollPanelNew.setBorder([]);

    % ===== PANEL CREATION =====
    % Controls list
    ctrl = struct('jScrollPanelTop',     jScrollPanelNew, ...
        'jPanelTop',           jPanelNew, ...
        ... ==== OPTICAL FLOW SETUP ====
        'jPanelSetup',         jPanelSetup, ...
        'jTextHornSchunck',    jTextHornSchunck, ...
        'jCheckHHD',           jCheckHHD, ...
        'jLabelRecursiveDepth',jLabelRecursiveDepth, ...
        'jTextRecursiveDepth', jTextRecursiveDepth, ...
        ... ==== OPTIONS ====
        'jPanelOptions',       jPanelOptions, ...
        'jCheckVisualize',     jCheckVisualize, ...
        'jCheckRotate',        jCheckRotate, ...
        ... ==== STABLE/TRANSIENT STATES ====
        'jCheckStates',        jCheckStates, ...
        'jCheckPublish',       jCheckPublish, ...
        'jCheckStatesOnly',    jCheckStatesOnly, ...
        ... ==== TIME INTERVAL PANEL =====
        'jPanelTime',          jPanelTime, ...
        'jTextTimeStart',      jTextTimeStart, ...
        'jTextTimeStop',       jTextTimeStop, ...
        'jLabelTimeUnit',      jLabelTimeUnit,  ...
        ... ==== Validation ====
        'jButtonRun',          jButtonRun, ...
        'jButtonCancel',       jButtonCancel);
    bst_mutex('create', panelName); % Return a mutex to wait for panel close
    bstPanelNew = BstPanel(panelName, jPanelNew, ctrl); % Form panel
    UpdatePanel(); % Update panel

%% ======================================================================
% === INTERNAL CALLBACKS ================================================
% =======================================================================

%% ===== CANCEL BUTTON =====
    function ButtonCancel_Callback(varargin)
        gui_hide(panelName); % Close panel
        bst_mutex('release', panelName); % Release the MUTEX
    end

%% ===== RUN BUTTON =====
    function ButtonRun_Callback(varargin)
        bst_mutex('release', panelName); % Release the MUTEX
    end

%% ===== UPDATE HHD =====
    function UpdatePanel()
        isDepth = jCheckHHD.isSelected(); % Does the user want HHD?
        jLabelRecursiveDepth.setEnabled(isDepth); % If so, allow them to ...
        jTextRecursiveDepth.setEnabled(isDepth); % choose recursive depth.

        isStates = jCheckStates.isSelected(); % Does the user want microstates?
        jCheckPublish.setEnabled(isStates); % If so, allow them to publish results
        jCheckStatesOnly.setEnabled(isStates); % If so, allow them to do states only

        isStatesOnly = jCheckStatesOnly.isSelected(); % Does user already have flow?
        jCheckRotate.setEnabled(~(isStates && isStatesOnly)); % If so, no rotation necessary
        jTextTimeStart.setEnabled(~(isStates && isStatesOnly)); % If so, no time start necessary
        jTextTimeStop.setEnabled(~(isStates && isStatesOnly)); % If so, no time end necessary
        jLabelHornSchunck.setEnabled(~(isStates && isStatesOnly)); % If so, no Horn-Schunk parameter label necessary
        jTextHornSchunck.setEnabled(~(isStates && isStatesOnly)); % If so, no Horn-Schunk parameter necessary

    end

end

%% ===== PROCESS PANEL INPUTS =====
function inputs = GetPanelContents() %#ok<DEFNU>
    % GETPANELCONTENTS Get inputs from panel to send to computation
    %
    % INPUTS:
    %   inputs              - structure containing all the inputs from the GUI
    %     .hornSchunck      - regularization parameter
    %     .HHDAvailable     - whether user wants HHD
    %     .depthHHD         - if user wants HHD, recursive depth
    %     .showResults      - immediately show optical flow overlaid after
    %                         computation
    %     .rotate           - calculate flow after normals to vertices are
    %                         rotated to point towards the center of the brain
    %     .segment          - calculate states of stable and fast flow
    %     .publishStates    - if states calculated, whether user wants to
    %                         publish states to PDF, etc.
    %     .tStart           - first time point to calculate optical flow
    %     .tEnd             - last time point to calculate optical flow
    ctrl = bst_get('PanelControls', 'OpticalFlow');

    % Horn-Schunck regularization parameter
    inputs.hornSchunck = str2double(char(ctrl.jTextHornSchunck.getText()));

    % HHD choice & amount of recursive depth
    if ctrl.jCheckHHD.isSelected
        inputs.HHDAvailable = true;
        inputs.depthHHD = str2double(char(ctrl.jTextRecursiveDepth.getText()));
    else
        inputs.HHDAvailable = false;
        inputs.depthHHD = -1;
    end

    % Options (show results, rotate results, overwrite results)
    inputs.showResults = ctrl.jCheckVisualize.isSelected();
    inputs.rotate = ctrl.jCheckRotate.isSelected();
    inputs.segment = ctrl.jCheckStates.isSelected();
    inputs.publishStates = ctrl.jCheckPublish.isSelected();
    inputs.statesOnly = ctrl.jCheckStatesOnly.isSelected();
    TimeUnit = char(ctrl.jLabelTimeUnit.getText);

    % First time point to calculate optical flow
    inputs.tStart = str2double(char(ctrl.jTextTimeStart.getText()));
    if strcmpi(TimeUnit, 'ms')
        inputs.tStart = inputs.tStart/1000;
    end
    % Last time point to calculate optical flow
    inputs.tEnd = str2double(char(ctrl.jTextTimeStop.getText()));
    if strcmpi(TimeUnit, 'ms')
        inputs.tEnd = inputs.tEnd/1000;
    end
end

%% ===== COMPUTE OPTICAL FLOW =====
function Compute(ResultsFile, inputs) %#ok<DEFNU>
    % COMPUTE       Compute optical flow and save to results
    % INPUTS:
    %   ResultsFile   - filename of results that have activity already
    %                   calculated, so we can add optical flow
    %   inputs        - structure of inputs from GUI or otherwise.
    %                   If otherwise, the structure must have
    %     .hornSchunck      - regularization parameter
    %     .HHDAvailable     - whether user wants HHD
    %     .depthHHD         - if user wants HHD, recursive depth
    %     .showResults      - immediately show optical flow overlaid after
    %                         computation
    %     .rotate           - calculate flow after normals to vertices are
    %                         rotated to point towards the center of the brain
    %     .segment          - calculate states of stable and fast flow
    %     .publishStates    - if states calculated, whether user wants to
    %                         publish states to PDF, etc.
    %     .statesOnly       - if states calculated and flow pre-calculated,
    %                         only calculate and display states.
    %     .overwrite        - overwrite old optical flow results
    %                         otherwise, add a second "result" (for testing
    %                         regularization, etc.)
    %     .tStart           - first time point to calculate optical flow
    %     .tEnd             - last time point to calculate optical flow
    % Read input file
    ResultsMat = in_bst_results(ResultsFile, 1, 'Time', 'SurfaceFile', 'ImageGridAmp');
    Time        = ResultsMat.Time;
    SurfaceFile = ResultsMat.SurfaceFile;
    SamplingInterval = Time(2)-Time(1);
    F = abs(double(ResultsMat.ImageGridAmp));
    clear ResultsMat;

    % Get inputs for optical flow
    if (nargin < 2) % Run GUI to get user inputs
        inputs = gui_show_dialog('Compute optical flow', @panel_opticalflow, 0, [], Time);
        if ~isfield(inputs, 'hornSchunck')% No inputs --> Cancel button was hit
            return
        end
    elseif ~isfield(inputs, 'hornSchunck') || ~isfield(inputs, 'HHDAvailable') ...
            || ~isfield(inputs, 'depthHHD') || ~isfield(inputs, 'showResults') ...
            || ~isfield(inputs, 'rotate') || ~isfield(inputs, 'overwrite') ...
            || ~isfield(inputs, 'tStart') || ~isfield(inputs, 'tEnd')

        bst_error(['Inputs not well-defined, structure must include:\n' ...
            'Horn-Schunck regularization parameter (hornSchunck)\n' ...
            'Desire to calculate HHD (HHDAvailable)\n' ...
            'HHD recursive depth parameter (depthHHD)\n' ...
            'Boolean to show results (showResults)\n' ...
            'Boolean to rotate results (rotate)\n' ...
            'Boolean to overwrite previous results (overwrite)\n' ...
            'First time to calculate optical flow (tStart)\n' ...
            'Last time to calculate optical flow (tEnd)']);
    end

    % Read tesselation
    FV = in_tess_bst(SurfaceFile);

    % Calculate flow
    ResultsMat = in_bst_results(ResultsFile, 1, 'OpticalFlow');
    if isempty(ResultsMat.OpticalFlow) || ~inputs.statesOnly % Calculate flow
        % Evaluate optical flow ...
        [flowField, int_dF, errorData, errorReg, poincare, timeInterval] = ...
            bst_call(@bst_opticalflow, F, FV, Time, ...
            inputs.tStart, inputs.tEnd, inputs.hornSchunck);
        if isfield(inputs, 'depthHHD') %  ... and, optionally, HHD
            % [U A H Vcurl Vdiv index] = ...
            %   HHD(dataFile, FV, Time, inputs.tStart, inputs.tEnd, inputs.depthHHD);
        end

        if inputs.rotate % Rotate flow so that tangent bundle is for circumsphere
            flowFieldRotated = rotate_optical_flow(flowField, FV.Vertices', FV.VertNormals');
        else
            flowFieldRotated = [];
        end

        % Save results into Results file as latest optical flow calculation
        save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ...
            Time, int_dF, errorData, errorReg, poincare);
    else
        flowField = ResultsMat.OpticalFlow.flowField;
        timeInterval = ResultsMat.OpticalFlow.timeInterval;
    end

    % Calculate states
    if inputs.segment
        bst_progress('start', 'Optical Flow', 'Segmenting into stable and transition states ...');

        interval = timeInterval(1) : SamplingInterval : timeInterval(2)+2*eps;
        [stableStates, transientStates, stablePoints, transientPoints, dEnergy] = ...
            bst_opticalflow_states(flowField, FV.Faces, FV.Vertices, 3, interval, SamplingInterval, true);

        % Save states
        Results = in_bst_results(ResultsFile);
        protocol = bst_get('ProtocolInfo');
        Results.OpticalFlow.dEnergy = dEnergy;
        Results.OpticalFlow.microstates = [];
        Results.OpticalFlow.microstates.stableStates = stableStates;
        Results.OpticalFlow.microstates.transientStates = transientStates;
        Results.OpticalFlow.microstates.stablePoints = stablePoints;
        Results.OpticalFlow.microstates.transientPoints = transientPoints;
        bst_save(bst_fullfile(protocol.STUDIES,ResultsFile), Results, 'v6');

        bst_progress('stop');
        if inputs.publishStates
            publish_states(ResultsFile, Results, Time, FV);
        end

    end

    % Show optical flow results on reconstructed map
    if inputs.showResults
        bst_memory('UnloadAll'); % Unload all previous datasets
        bst_call(@view_surface_data, SurfaceFile, ResultsFile);
        panel_time('SetCurrentTime', timeInterval(1)+eps);
    end
end

%% ===== ROTATE OPTICAL FLOW TO INFLATED SPHERE =====
function flowFieldRotated = rotate_optical_flow(flowField, Vertices, VertNormals)
    % ROTATE_OPTICAL_FLOW     Rotate optical flow to preserve angle against
    %                         normal-to-vertex if the normal is rotated
    %                         to point towards the center of the brain.
    %                         Equivalently, rotate optical flow so result
    %                         is as if brain was sphere (for addition of
    %                         flow between neighboring vertices, perhaps.)
    % INPUTS:
    %   flowField           - optical flow at every vertex
    %   Vertices            - locations of vertices
    %   VertNormals         - normals to each vertex
    % OUTPUTS:
    %   flowFieldRotated    - Rotate normals-to-vertex to point to sphere,
    %                         and then rotate optical flow using same
    %                         transformation to get this result
    nVertices = size(Vertices,1);
    centeredVertices = Vertices - repmat(mean(Vertices), nVertices, 1);
    flowFieldRotated = zeros(size(flowField));

    bst_progress('start', 'Optical Flow', ...
        'Rotating flows for visualization ... ', 0, nVertices);
    for m = 1:nVertices
        c = VertNormals(m,:); d = -centeredVertices(m,:);
        current = c/norm(c); desired = d/norm(d);
        perpCurrent = cross(desired,current)/norm(cross(desired,current));
        perpDesired = cross(perpCurrent,desired);
        frameChange = [desired' perpDesired' perpCurrent'];
        rotation = [dot(current,desired) -dot(current,perpDesired) 0; ...
            dot(current,perpDesired) dot(current,desired) 0; ...
            0 0 1];
        transform = (frameChange*rotation)/frameChange;
        flow = squeeze(flowField(m,:,:))';
        flowFieldRotated(m,:,:) = transpose(flow*transform);
        if mod(m,20) == 0
            bst_progress('inc', 20);
        end
    end
    bst_progress('stop');
end

%% ===== SAVE OPTICAL FLOW INTO BRAINSTORM =====
function save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ...
        Time, int_dF, errorData, errorReg, poincare)
    % SAVE_FLOW     Save flow results (and possibly publish microstates)
    % INPUTS:
    %   ResultsFile       - File containing results (and surface file)
    %   inputs            - inputs to error-check whether results are written
    %   flowField         - Optical flow field
    %                       dimension (# of vertices) X length(tStart:tEnd)
    %   flowFieldRotated  - flow field rotated for faux spherical brain
    %   Time              - time when activity was reconstructed (including
    %                       times for which optical flow is not calculated)
    %   int_dF            - Constant term in variational formulation
    %   errorData         - Error in fit to data
    %   errorReg          - Energy in regularization
    %   poincare          - Poincar� index

    opticalFlow.flowField = flowField; % Optical flow results
    opticalFlow.flowFieldRotatedAvailable = inputs.rotate;
    if inputs.rotate % Results with surface normal rotated towards center of boundary's volume
        opticalFlow.flowFieldRotated = flowFieldRotated;
    end
    opticalFlow.timeInterval = [inputs.tStart inputs.tEnd]; % Time interval
    opticalFlow.samplingInterval = Time(2)-Time(1); % Time interval
    opticalFlow.hornSchunck = inputs.hornSchunck; % Regularization
    opticalFlow.int_dF = int_dF;
    opticalFlow.errorData = errorData; % Error in fit to data
    opticalFlow.errorReg = errorReg; % Error from smooth regularization
    opticalFlow.poincare = poincare;
    opticalFlow.HHDAvailable = inputs.HHDAvailable;
    if inputs.HHDAvailable
        opticalFlow.depthHHD = inputs.depthHHD;
    end

    % Save optical flow results in original results file
    Results = in_bst_results(ResultsFile);
    if opticalFlow(end).HHDAvailable
        Results = bst_history('add', Results, 'compute', ...
            ['Optical flow estimated: [' int2str(inputs.tStart*1000) ...
            ', ' int2str(inputs.tEnd*1000) ']ms']);
    else
        Results = bst_history('add', Results, 'compute', ...
            ['Optical flow & HHD estimated: [' ...
            int2str(inputs.tStart*1000) ...
            ', ' int2str(inputs.tEnd*1000) ']ms']);
    end
    Results.OpticalFlow = opticalFlow;

    Protocol = bst_get('ProtocolInfo');
    bst_save(bst_fullfile(Protocol.STUDIES, ResultsFile), Results, 'v6');

end

%% ===== PUBLISH STATES OF FLOW =====
function publish_states(ResultsFile, Results, Time, FV)
    % PUBLISH_STATES    Publish states to an output file (PDF, etc.)
    %
    % INPUTS:
    %   ResultsFile         - File containing results (and surface file)
    %   Results             - results structure (so we don't have to re-load it)
    %   Time                - time points of reconstructed data
    %   FV                  - tesselation for plotting

    opticalFlow = Results.OpticalFlow;

    % Get extrema
    tStartIndex = find(Time < opticalFlow.timeInterval(1)-eps, 1, 'last')+1; % Index of first time point for flow calculation
    tEndIndex = find(Time < opticalFlow.timeInterval(2)-eps, 1, 'last')+1; % Index of last time point for flow calculation
    activity = Results.ImageGridAmp(:,tStartIndex:tEndIndex); % Activity in flow-calculated interval
    extrema = sortrows([ ...
        opticalFlow.microstates.stableStates zeros(length(opticalFlow.microstates.stableStates), 1); ...
        opticalFlow.microstates.transientStates ones(length(opticalFlow.microstates.transientStates),1) ...
        ]); % Sort extrema for timeline visual

    % Setup data for visualization
    [opticalFlowDir, iDS, FigureId] = ...
        publish_setup(ResultsFile, Results.Options.DataTypes);
    if isempty(opticalFlowDir) || ~ischar(opticalFlowDir) % Cancel button!
        return
    end
    if ispc && opticalFlowDir(end) ~= '\'
        opticalFlowDir(end+1) = '\';
    elseif opticalFlowDir(end) ~= '/'
        opticalFlowDir(end+1) = '/';
    end

    % Pre-assigned inputs for demonstration
    viewChoice = [-1.5 0.5 1]; zoomIn = 2;
    bst_progress('start', 'Optical Flow', ['Saving images of states to ...' opticalFlowDir], 0, length(extrema));
    for m = 1:length(extrema)

        if extrema(m,3) > 1-eps % Transition state
            opticalFlowTag = ['transient_' ...
                sprintf('%0.2d', sum(extrema(1:m,3) > 1-eps)) '_' ...
                sprintf('%0.3d', round(Time(extrema(m,1) + tStartIndex - 1)*1000)) 'ms_to_' ...
                sprintf('%0.3d', round(Time(extrema(m,2) + tStartIndex - 1)*1000)) 'ms'];
        else % Stable state
            opticalFlowTag = ['stable_' ...
                sprintf('%0.2d', sum(extrema(1:m,3) > 1-eps)) '_' ...
                sprintf('%0.3d', round(Time(extrema(m,1) + tStartIndex - 1)*1000)) 'ms_to_' ...
                sprintf('%0.3d', round(Time(extrema(m,2) + tStartIndex - 1)*1000)) 'ms'];
        end
        bst_progress('text', ['Writing: ' opticalFlowTag]); % inform user

        % Initialize figure
        [hState, iSurfaceState] = ...
            publish_initialize(iDS, FigureId, ResultsFile, ...
            Results.SurfaceFile, max(abs(activity(:))), viewChoice);
        set(hState, 'PaperPositionMode', 'auto');

        % Plot figure
        if extrema(m,3) > 1-eps % Transition state
            publish_transient(hState, FV.Vertices, opticalFlow, extrema(m,1), zoomIn);
        else % Stable state
            publish_stable(hState, iSurfaceState, activity(:, extrema(m,1):extrema(m,2)));
        end

        % Save and close figure
        print(hState, '-dpng', '-r300', [opticalFlowDir opticalFlowTag]);
        if m < length(extrema)
            bst_figures('DeleteFigure', hState, 'NoUnload') % Delete figure
        else
            bst_figures('DeleteFigure', hState) % Delete figure and unload dataset
        end

        bst_progress('inc', 1); % Update waitbar
    end

    bst_progress('stop'); % Finish waitbar
end

function [opticalFlowDir, iDS, FigureId] = publish_setup(ResultsFile, DataTypes)
    % PUBLISH_SETUP   Setup for publishing stable and transition states
    % INPUTS:
    %   ResultsFile   - filename of reconstructed source signals
    %   DataTypes     - list of data types in results
    % OUTPUTS:
    %   iDS           - index of results in BrainStorm
    %   FigureId      - ID for figure

    bst_memory('UnloadAll'); % Unload all previous datasets

    % Select folder to drop images off
    protocol = bst_get('ProtocolInfo');
    if isempty(strfind(ResultsFile, '/'))
        defaultDir = [protocol.STUDIES '\' ResultsFile(1:(find(ResultsFile == '\', 1, 'last')))];
    else
        defaultDir = [protocol.STUDIES '/' ResultsFile(1:(find(ResultsFile == '/', 1, 'last')))];
    end
    opticalFlowDir = uigetdir(defaultDir, ... % Open 'Select directory' dialog
        'Please select directory to save state maps (default = location of results)');
    if isempty(opticalFlowDir) || ~ischar(opticalFlowDir)
        iDS = -1;
        FigureId = -1;
        return
    end

    iDS = bst_memory('LoadResultsFile', ResultsFile); % Dataset index

    % Modality for classifying figure
    AllModalities = DataTypes;
    if all(ismember({'MEG GRAD', 'MEG MAG'}, AllModalities))
        AllModalities{end+1} = 'MEG';
        AllModalities = setdiff(AllModalities, {'MEG GRAD', 'MEG MAG'});
    end
    Modality = AllModalities{1};

    % Create figure ID setup
    FigureId = db_template('FigureId');
    FigureId.Type     = '3DViz';
    FigureId.SubType  = '';
    FigureId.Modality = Modality;
end

function [hState, iSurfaceState] = publish_initialize(iDS, FigureId, ResultsFile, SurfaceFile, DataMinMax, viewChoice)
    % PUBLISH_INITIALIZE    Figure w/ surface for plotting flows or activities
    % INPUTS:
    %   DataTypes     - list of data types in results
    %   iDS           - index of results in BrainStorm
    %   FigureId      - ID for figure
    %   ResultsFile   - filename of reconstructed source signals
    %   SurfaceFile   - filename of surface to plot results on
    %   DataMinMax    - maximum of magnitude of activity for colorbar
    %   viewChoice    - angle to view surface
    % OUTPUTS:
    %   hState          - Figure containing surface to plot flows on
    %   iSurfaceState   - index of surface for sulci, etc.

    % Create and initialize figure with surface
    hState = bst_figures('CreateFigure', iDS, FigureId, 'AlwaysCreate');
    iSurfaceState = panel_surface('AddSurface', hState, SurfaceFile); % Results.SurfaceFile
    setappdata(hState, 'ResultsFile', ResultsFile);

    % Add colormap for the colorbar
    TessInfo = getappdata(hState, 'Surface');
    TessInfo(iSurfaceState).DataMinMax = DataMinMax; % max(abs(activity(:)));
    TessInfo(iSurfaceState).ColormapType = 'source';
    TessInfo(iSurfaceState).SurfSmoothValue = 1; % Completely smoothed surface
    TessInfo(iSurfaceState).SurfShowSulci = true; % Show sulci
    bst_colormaps('AddColormapToFigure', hState, TessInfo(iSurfaceState).ColormapType);

    % Set view just like in BrainStorm (ish)
    hAxes = findobj(hState, '-depth', 1, 'Tag', 'Axes3D'); % Get Axes handle
    view(hAxes, viewChoice); % Update view
    if abs(viewChoice(3)) > abs(viewChoice(2)) + abs(viewChoice(1))
        camup(hAxes, [1 0 0]); % Update camera position for bottom/top
    else
        camup(hAxes, [0 0 1]); % Update camera position for others
    end
    camlight(findobj(hAxes, '-depth', 1, 'Tag', 'FrontLight'), 'headlight'); % Update head light position
end

function publish_transient(hState, Vertices, opticalFlow, start, zoomIn)
    % PUBLISH_TRANSIENT   Show flows in transient state
    % INPUTS:
    %   hState          - Figure containing surface to plot flows on
    %   iSurfaceState   - index of surface for sulci, etc.
    %   Vertices        - locations of vertices in tesselation
    %   opticalFlow     - flow results
    %   start           - flow at start of state

    % Hold axes to plot on top of surface
    hAxes = findobj(hState, '-depth', 1, 'Tag', 'Axes3D'); % Get Axes handle
    hold(hAxes,'on');
    flowField = opticalFlow.flowField(:,:,start);
    useful = sum(flowField.^2, 2) > max(sum(flowField.^2, 2))*0.1;
    flowField(~useful,:) = 0;
    quiver3(hAxes, ...
        Vertices(:,1), Vertices(:,2), Vertices(:,3), ...
        flowField(:,1), flowField(:,2), flowField(:,3), ...
        6, 'c', 'LineWidth', 2);
    hold(hAxes,'off');

    zoom(zoomIn); % Zoom in plz!

    % % Title of state
    % tagLocation = [mean(FV.Vertices(:,1))-0.025 max(FV.Vertices(:,2))+0.01 mean(FV.Vertices(:,3))-0.03];
    % text(tagLocation(1), tagLocation(2), tagLocation(3), 'TRANSIENT STATE', 'fontsize', 24, 'color', [1 0 0])

end

function publish_stable(hState, iSurfaceState, activity)
    % PUBLISH_STABLE      Show mean activity in stable state
    % INPUTS:
    %   hState          - Figure containing surface to plot flows on
    %   iSurfaceState   - index of surface for sulci, etc.
    %   activity        - reconstructed activity in state
    %   DataMinMax      - maximum of colorbar (as used in initialization)

    % Add overlay of mean activity
    TessInfo = getappdata(hState, 'Surface');
    TessInfo(iSurfaceState).Data = mean(abs(activity), 2); % Plot mean activity
    TessInfo(iSurfaceState).DataMinMax = max(TessInfo(iSurfaceState).Data);
    TessInfo(iSurfaceState).DataLimitValue = [0 TessInfo(iSurfaceState).DataMinMax];
    setappdata(hState, 'Surface', TessInfo);

    % Set maximum value of activity
    hAxes = [findobj(hState, '-depth', 1, 'Tag', 'Axes3D'), ...
        findobj(hState, '-depth', 1, 'Tag', 'axc'), ...
        findobj(hState, '-depth', 1, 'Tag', 'axa'), ...
        findobj(hState, '-depth', 1, 'Tag', 'axs')]; % Figure axes
    set(hAxes, 'CLim', TessInfo(iSurfaceState).DataLimitValue);

    % Add colorbar
    ColormapInfo = getappdata(hState, 'Colormap'); % Default colormap
    sColormap = bst_colormaps('GetColormap', ColormapInfo.Type); % Get figure colormap
    set(hState, 'Colormap', sColormap.CMap); % Set figure colormap
    bst_colormaps('SetColorbarVisible', hState, sColormap.DisplayColorbar); % Create/Delete colorbar
    bst_colormaps('ConfigureColorbar', hState, ColormapInfo.Type, 'opticalflow', ColormapInfo.DisplayUnits); % Display only one colorbar (preferentially the results colorbar)

    % Reveal plot!
    figure_3d('UpdateSurfaceColor', hState, iSurfaceState);

    % % Title of state
    % tagLocation = [mean(FV.Vertices(:,1)) mean(FV.Vertices(:,2)) max(FV.Vertices(:,3)) + 0.01];
    % text(tagLocation(1)-0.04, tagLocation(2), tagLocation(3), 'STABLE STATE', 'fontsize', 24, 'color', [0 1 0])
end

%% ===== PLOT OPTICAL FLOW ON BRAINSTORM CURRENT DENSITY FIGURE =====
function PlotOpticalFlow(hFig, opticalFlow, currentTime, sSurf)
    % PLOTOPTICALFLOW     From BrainStorm figure handler, plot optical flow
    %                     results on top of surface with data
    % INPUTS:
    %   hFig          - handle to current figure
    %   opticalFlow   - optical flow results
    %   currentTime   - time point shown in BrainStorm
    %   sSurf         - surface file (to plot arrows in the right spots)

    % Process figure (removing old flows if necessary + getting surface axes)
    nVertices = size(sSurf.Vertices, 1);
    [ax, currentName] = process_surface(hFig);

    % First check if we need to do anything
    flagPlotFlow = 0;
    for n = 1:length(opticalFlow)
        if currentTime > (opticalFlow(n).timeInterval(1)-(opticalFlow(n).samplingInterval/2)) ...
                && currentTime < (opticalFlow(n).timeInterval(2)+(opticalFlow(n).samplingInterval/2))
            flagPlotFlow = 1;
        end
    end
    if ~flagPlotFlow
        return
    end

    % Add ability to see sphere-rotated results if desired
    if opticalFlow.flowFieldRotatedAvailable
        hRotated = checkbox_rotated(hFig, flagPlotFlow, opticalFlow, currentTime, sSurf);
        plotRotatedResults = get(hRotated, 'Value');
    else
        plotRotatedResults = false;
        hRotated = -1;
    end

    button_movie(hFig, hRotated, ... % Push button to play movie!
        opticalFlow.timeInterval(1) : opticalFlow.samplingInterval : opticalFlow.timeInterval(2)+eps);

    % Plot optical flow in a shell outside the brain for rotated results
    if plotRotatedResults
        meanVertex = repmat(mean(sSurf.Vertices), nVertices, 1);
        centeredVertices = sSurf.Vertices - meanVertex;
        radii = sqrt(sum(centeredVertices.^2, 2));
        maxRadius = max(radii) * 1.1;
        Vertices = centeredVertices ./ repmat(radii, 1, 3) * maxRadius + meanVertex;
    else
        Vertices = sSurf.Vertices + sSurf.VertNormals * 0.005;
    end

    % Get time index
    timeIdx = round((currentTime-opticalFlow.timeInterval(1))/opticalFlow.samplingInterval)+1;
    if timeIdx > size(opticalFlow.flowField, 3)
        return
    end

    % Remove old quivers (arrows)
    oldQuivers = findobj(ax, 'Type', 'quiver');
    for iQuiver = 1:length(oldQuivers)
        delete(oldQuivers(iQuiver));
    end

    % Hold axes to plot on top of surface
    hold(ax,'on');
    if plotRotatedResults
        flowField = opticalFlow.flowFieldRotated(:,:,timeIdx);
    else
        flowField = opticalFlow.flowField(:,:,timeIdx);
    end
    useful = sum(flowField.^2, 2) > max(sum(flowField.^2, 2))*0.1;
    flowField(~useful,:) = 0;
    quiver3(ax, ...
        Vertices(:,1), Vertices(:,2), Vertices(:,3), ...
        flowField(:,1), flowField(:,2), flowField(:,3), ...
        6, 'c', 'LineWidth', 2); % Color is cyan, works well with hot colormap
    hold(ax,'off');

    % Modify figure name if we are in stable/transient state
    if isfield(opticalFlow, 'microstates')
        stablePoints = opticalFlow.microstates.stablePoints;
        transientPoints = opticalFlow.microstates.transientPoints;
        if find(abs(stablePoints-timeIdx) <= eps)
            currentName = ['STABLE: ' currentName];
        elseif find(abs(transientPoints-timeIdx) <= eps)
            currentName = ['TRANSIENT: ' currentName];
        end
    end
    set(hFig, 'Name', currentName); % Set figure name to microstate label

end

%% ===== PROCESS SURFACE AND CLEAN PREVIOUS FLOWS =====
function [ax, currentName] = process_surface(hFig)
    % PROCESS_SURFACE     Clean surface:
    %                     * get rid of old results
    %                     * get rid of old state results in name of figure
    %                     * get axis handle for surface so we can overlay
    %                       optical flow results
    % INPUTS:
    %   hFig          - handle to figure containing activity
    % OUTPUTS:
    %   ax            - axis containing surface (for plotting flow)
    %   currentName   - name on title of figure

    % Get axis handle for surface
    axes = get(hFig, 'children');
    ax = [];
    for n = 1:length(axes)
        if isprop(axes(n), 'CLim')
            cLim = get(axes(n), 'CLim');
            if cLim(1) ~= 1 || cLim(2) ~= 256
                ax = axes(n);
                break;
            end
        end
    end

    % Clean off previous vector fields
    if ~isempty(ax)
        hOld = get(ax, 'Children');
        for n = 1:length(hOld)
            if strcmp(get(hOld(n), 'Type'), 'hggroup') % hggroup is type name for quiver plot
                delete(hOld(n));
            end
        end
    end

    % Clean off stable/transient labeling
    currentName = get(hFig, 'Name');
    if strfind(currentName, 'STABLE: ')
        currentName(strfind(currentName, 'STABLE: ') + (0:7)) = [];
    elseif strfind(currentName, 'TRANSIENT: ')
        currentName(strfind(currentName, 'TRANSIENT: ') + (0:10)) = [];
    end

end

%% ===== MOVIE FOR ITERATING THROUGH FLOWS =====
function hMovieButton = button_movie(hFig, hRotated, flowInterval)
    % BUTTON_MOVIE    Play a loop of flow within some part of interval
    % INPUTS:
    %   hFig          - Figure containing reconstructions and flows
    %   hRotated      - Handle to rotated-flows checkbox for repositioning
    %   flowInterval  - interval over which flow has been calculated
    % OUTPUTS:
    %   hMovieButton  - Handle to button for Movie? or Stop Movie

    % Find if button already exists on figure
    hMovieButton = -1;
    figureParts = get(hFig, 'children');
    for m = 1:length(figureParts)
        if strcmp(get(figureParts(m),'Type'), 'uicontrol') && ...
                strcmp(get(figureParts(m),'Tag'), 'MovieButton')
            hMovieButton = figureParts(m);
        end
    end

    % Make new button if necessary
    if hMovieButton < 0
        hMovieButton = uicontrol(hFig, 'Style', 'pushbutton', 'String', 'Movie?', ...
            'Position', [20 10 70 20], 'Tag', 'MovieButton', ...
            'Callback', {@ButtonMovie_Callback, hFig, flowInterval, hRotated});
    end

    function ButtonMovie_Callback(hMovieButton, ev, hFig, flowInterval, hRotated)
        % BUTTONMOVIE_CALLBACK  Callback function to the movie button
        %
        % INPUTS:
        %   hMovieButton  - Handle to button
        %   event         - ~
        %   hFig          - Handle to figure showing reconstructions and flows
        %   flowInterval  - interval over which flow has been calculated
        %   hRotated      - Handle to rotated-flows checkbox for repositioning

        if strcmp(get(hMovieButton, 'String'), 'Movie?') % Start movie controls

            set(hMovieButton,  'String', 'Stop Movie'); % Change button to stop movie
            set(hMovieButton, 'Value', 1); % HACK ALERT: current time index stored here
            halt_others(hFig); % Disallow all other figures from trying to play a movie

            if hRotated > -1 % If rotated-flows checkbox available, shift it to the right
                set(hRotated, 'Position', [370 10 130 20]);
            end

            % Start and end of movie loop
            hMovieStart = uicontrol(hFig, 'Style', 'edit', ...
                'String', 'Start', ...
                'Position', [100 10 60 20], ...
                'Tag', 'MovieStart', ...
                'BackgroundColor', [1 1 1], ...
                'HorizontalAlignment', 'right', ...
                'Callback', {@validate_time, flowInterval});
            uicontrol(hFig, 'Style', 'text', 'String', 'ms', 'Tag', 'TextMovieStart', 'Position', [165 10 20 20]);
            hMovieEnd = uicontrol(hFig, 'Style', 'edit', ...
                'String', 'End', ...
                'Position', [190 10 60 20], ...
                'Tag', 'MovieEnd', ...
                'BackgroundColor', [1 1 1], ...
                'HorizontalAlignment', 'right', ...
                'Callback', {@validate_time, flowInterval});
            uicontrol(hFig, 'Style', 'text', 'String', 'ms', 'Tag', 'TextMovieEnd', 'Position', [255 10 20 20]);

            % Movie controls: play and pause
            hMoviePlay = uicontrol(hFig, 'Style', 'pushbutton', ...
                'String', 'Play', ... % BOO NO PLAY BUTTON YOU SUCK MATLAB
                'Tag', 'MoviePlay', ...
                'Position', [280 10 40 20]);
            hMoviePause = uicontrol(hFig, 'Style', 'pushbutton', ...
                'String', 'Pause', ... % BOO NO PAUSE BUTTON YOU SUCK MATLAB
                'Tag', 'MoviePause', ...
                'Position', [325 10 40 20], ...
                'Enable', 'off');
            set(hMoviePlay, 'Callback', ...
                {@ButtonMoviePlay_Callback, hMovieButton, hMoviePause, ...
                hMovieStart, hMovieEnd, flowInterval(2)-flowInterval(1)});
            set(hMoviePause, 'Callback', {@ButtonMoviePause_Callback, hMoviePlay});

        else % Stop movie
            set(hMovieButton, 'FontName', 'arial', 'String', 'Movie?'); % Change label to play movie

            % Clean everything
            figureParts = get(hFig, 'children');
            for n = 1:length(figureParts)
                if strcmp(get(figureParts(n),'Type'), 'uicontrol') && ...
                        (strcmp(get(figureParts(n),'Tag'), 'MovieStart') || ...
                        strcmp(get(figureParts(n),'Tag'), 'TextMovieStart') || ...
                        strcmp(get(figureParts(n),'Tag'), 'MovieEnd') || ...
                        strcmp(get(figureParts(n),'Tag'), 'TextMovieEnd') || ...
                        strcmp(get(figureParts(n),'Tag'), 'MoviePlay') || ...
                        strcmp(get(figureParts(n),'Tag'), 'MoviePause'))
                    delete(figureParts(n));
                end
            end

            allow_others(hFig); % Allow other figures to start movie

            if hRotated > -1 % Move rotated-flows checkbox back
                set(hRotated, 'Position', [100 10 130 20]);
            end

            % Set current time back to first time of flow
            panel_time('SetCurrentTime', flowInterval(1)+eps);
        end

    end

    function validate_time(hText, ev, interval)
        % VALIDATE_TIME   Check if time asked for by user is legitimately in
        %                 interval for which flow is calculated, and correct
        %                 if necessary
        %
        % INPUTS:
        %   hText     - handle containing text element for verificatino
        %   event     - ~
        %   interval  - hText must be an element in interval

        chosen = str2double(get(hText, 'String'))/1000; % Get time
        if isnan(chosen) % NaN -> bad input -> go back to default
            if strcmp(get(hText,'Tag'), 'MovieStart')
                chosen = interval(1);
            elseif strcmp(get(hText,'Tag'), 'MovieEnd')
                chosen = interval(end);
            end
        else
            [tmp,idx] = min(abs(interval-chosen)); % Find closest time point ...
            chosen = interval(idx); % ... and use that as updated value
        end
        set(hText, 'String', sprintf('%0.2f', chosen*1000)); % Show user update
    end

    function ButtonMoviePlay_Callback(hMoviePlay, ev, ...
            hMovieButton, hMoviePause, hMovieStart, hMovieEnd, movieStep)
        set(hMoviePlay, 'Enable', 'off'); % Can't play twice
        set(hMoviePause, 'Enable', 'on'); % Allow pausing
        flowInterval = ... % Interval of times that will be looping
            str2double(get(hMovieStart,'String'))/1000 : ... % Start here
            movieStep : ... % Go this fast
            (str2double(get(hMovieEnd,'String'))/1000)+2*eps; % End here

        while ishandle(hMoviePause) && strcmp(get(hMoviePause, 'Enable'), 'on') % Pause only when pause button is enabled

            currentTimeIdx = get(hMovieButton, 'Value'); % HACK ALERT: current time index stored here
            if currentTimeIdx == 0
                currentTimeIdx = 1; % Correction for first time (sometimes)
            end
            panel_time('SetCurrentTime', flowInterval(currentTimeIdx)); % Update
            drawnow; pause(0.5); % 120fps
            currentTimeIdx = currentTimeIdx+1; % Advance
            if currentTimeIdx > length(flowInterval) % Wrap around for loop
                currentTimeIdx = 1;
            end

            if ishandle(hMovieButton)
                set(hMovieButton, 'Value', currentTimeIdx); % HACK ALERT: current time index stored here after pausing
            end
        end

    end

    function ButtonMoviePause_Callback(hMoviePause,ev,hMoviePlay)
        % BUTTONMOUSEPAUSE_CALLBACK   Pauses looping of flows
        %
        % INPUTS
        %   hMoviePlay    - Button for playing it. Enable it so user can
        %                   continue as desired
        %   hMoviePause   - Button for pausing. Disable it to pause loop
        set(hMoviePlay, 'Enable', 'on');
        set(hMoviePause, 'Enable', 'off');
    end

    function halt_others(hFig)
        % HALT_OTHERS   Prevent other figures from trying to play a movie
        %
        % INPUTS:
        %   hFig  - Figure whose movie controls should NOT be touched
        h = findall(0, 'Type', 'figure'); % All figures
        h = setdiff(h, hFig); % Do not touch hFig!
        for o = 1:length(h)
            if ~strcmp(get(h(o), 'Tag'), '3DViz') % Not a reconstruction
                continue
            end

            figureParts = get(h(o), 'children');
            for p = 1:length(figureParts) % Delete all movie controls
                if strcmp(get(figureParts(p),'Type'), 'uicontrol') && ...
                        (strcmp(get(figureParts(p),'Tag'), 'MovieStart') || ...
                        strcmp(get(figureParts(p),'Tag'), 'TextMovieStart') || ...
                        strcmp(get(figureParts(p),'Tag'), 'MovieEnd') || ...
                        strcmp(get(figureParts(p),'Tag'), 'TextMovieEnd') || ...
                        strcmp(get(figureParts(p),'Tag'), 'MoviePlay') || ...
                        strcmp(get(figureParts(p),'Tag'), 'MoviePause'))
                    delete(figureParts(n));
                elseif strcmp(get(figureParts(p),'Type'), 'uicontrol') && ...
                        (strcmp(get(figureParts(p),'Tag'), 'MovieButton'))
                    set(figureParts(p), 'Enable', 'off'); % Disable movie button
                end
            end

        end

    end

    function allow_others(hFig)
        % ALLOW_OTHERS   Turn on all other movie buttons (when movie is stopped)
        %
        % INPUTS:
        %   hFig  - Figure whose movie controls should NOT be touched

        h = findall(0, 'Type', 'figure'); % All figures ...
        h = setdiff(h, hFig); % ... except the one that just stopped
        for o = 1:length(h)
            if ~strcmp(get(h(o), 'Tag'), '3DViz') % Not a reconstruction
                continue
            end

            figureParts = get(h(o), 'children');
            for p = 1:length(figureParts) % Enable all movie buttons
                if strcmp(get(figureParts(p),'Type'), 'uicontrol') && ...
                        (strcmp(get(figureParts(p),'Tag'), 'MovieButton'))
                    set(figureParts(p), 'Enable', 'on');
                end
            end

        end
    end

end

%% ===== CHECKBOX FOR PLOTTING FLOWS ON SPHERE OR SURFACE =====
function hRotated = checkbox_rotated(hFig, flagPlotFlow, opticalFlow, currentTime, sSurf)
    % CHECKBOX_ROTATED  Create checkbox on surface figure with data, so that
    %                   we can choose whether to see the flows on the surface
    %                   or rotated for addition
    % Inquire whether checkbox exists
    figureParts = get(hFig, 'children');
    for n = 1:length(figureParts)
        if strcmp(get(figureParts(n),'Type'), 'uicontrol') && ...
                strcmp(get(figureParts(n),'Tag'), 'CheckboxRotated')
            hRotated = figureParts(n);
        end
    end

    if ~flagPlotFlow
        if exist('hRotated', 'var') % No flow plotted so ...
            delete(hRotated); % ... remove checkbox
        else
            hRotated = -1; % Send some value back
        end
    elseif flagPlotFlow && ~exist('hRotated', 'var') % Flow plotted so ...
        hRotated = uicontrol(hFig, 'Style', 'checkbox', ... % ... add checkbox
            'String', 'Display rotated flows', ...
            'Value', 0, ...
            'Tag', 'CheckboxRotated', ...
            'Position', [100 10 130 20], ...
            'Callback', {@CheckboxRotated_Callback, hFig, opticalFlow, currentTime, sSurf});
    end

        function CheckboxRotated_Callback(h, ev, hFig, opticalFlow, currentTime, sSurf)
            PlotOpticalFlow(hFig, opticalFlow, currentTime, sSurf);
        end
end

author: brainstorm-tools


text:
function CPD = reset_ess(CPD)
% RESET_ESS Reset the Expected Sufficient Statistics for a CPD (mlp)
% CPD = reset_ess(CPD)

CPD.W1 = [];
CPD.W2 = [];
CPD.b1 = [];
CPD.b2 = [];
CPD.parent_vals = [];
CPD.eso_weights=[];
CPD.self_vals = [];
CPD.nsamples = 0;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ nin, pi ] = triangle_line_imp_int_2d ( a, b, c, t, nin, pi )

%*****************************************************************************80
%
%% TRIANGLE_LINE_IMP_INT_2D finds where an implicit line intersects a triangle in 2D.
%
%  Discussion:
%
%    An implicit line is the set of points ( X, Y ) satisfying
%
%      A * X + B * Y + C = 0
%
%    where at least one of A and B is not zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, determine the equation of the line:
%    A*X + B*Y + C = 0.
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, integer NIN, the number of points of intersection
%    of the line with the triangle.  NIN may be 0, 1, 2 or 3.
%
%    Output, real PI(2,3), contains the intersection points.
%
  dim_num = 2;

  nin = 0;

  for i = 1 : 3

    j = i4_wrap ( i+1, 1, 3 );
%
%  Get the implicit form of the line through vertices I and I+1.
%
    [ a1, b1, c1 ] = line_exp2imp_2d ( t(1:2,i), t(1:2,j) );
%
%  Seek an intersection with the original line.
%
    [ ival, p ] = lines_imp_int_2d ( a, b, c, a1, b1, c1 );
%
%  If there is an intersection, determine if it lies between the two vertices.
%
    if ( ival == 1 )

      test1 = ( p(1:dim_num)'  - t(1:dim_num,i) )' ...
        * ( t(1:dim_num,j) - t(1:dim_num,i) );
      test2 = ( t(1:dim_num,j) - t(1:dim_num,i) )' ...
        * ( t(1:dim_num,j) - t(1:dim_num,i) );

      if ( 0 <= test1 & test1 <= test2 )
        nin = nin + 1;
        pi(1:dim_num,nin) = p(1:dim_num)';
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dist = imDistanceMap(img, varargin)
%IMDISTANCEMAP Compute chamfer distance using scanning algorithm
%
%   Usage:
%   RES = imDistanceMap(BIN);
%   where BIN is a binary images, computes for each foreground pixel the
%   minimum distance to the background.
%   The function propagates distances to orthogonal and diagonal pixels,
%   using weights equal to 1 for orthogonal pixels, and sqrt(2) for
%   diagonal markers. The result RES is given in a double array the same
%   size as BIN.
%
%   RES = imDistanceMap(LBL);
%   The function can be applied to a label image LBL, where the values of
%   LBL correspond to region indices. In that case, the result RES
%   corresponds to the superposition of the results on computed on each
%   label.
%
%
%   Methods:
%   The function uses scanning algorithm, by applying one forward and one
%   backward scan.
%
%   RES = imDistanceMap(..., WEIGHTS);
%   Specifies different weights for computing distance between 2 pixels.
%   WEIGHTS is a 2 elements array, with WEIGHTS(1) corresponding to the
%   distance between two orthonal pixels, and WEIGHTS(2) corresponding to
%   the distance between two diagonal pixels.
%   Possible choices
%   WEIGHTS = [1 sqrt(2)]   -> quasi-euclidean distance
%   WEIGHTS = [1 Inf]       -> "Manhattan" or "cityblock" distance
%   WEIGHTS = [1 1]         -> "Chessboard" distance
%   WEIGHTS = [3 4]         -> Borgerfors' weights
%   WEIGHTS = [5 7]         -> close approximation of sqrt(2)
%   WEIGHTS = [5 7 11]      -> Uses an additional weight for chess-knight
%                              shifts around each pixel( default)
%
%   Note: when specifying weights, the result has the same class/data type
%   than the array of weights. It is possible to balance between speed and
%   memory usage:
%   - if weights are double (the default), the memory usage is larger, but
%       the result can be given in pixel units
%   - if weights are integer (for Borgefors weights, for example), the
%       memory usage is reduced, but representation limit of datatype can
%       be easily reached. One needs to divide by the first weight to get
%       result comparabale with natural distances.
%       For uint8, using [3 4] weigths, the maximal computable distance is
%       around 255/3 = 85 pixels. Using 'int16'  seems to be a good
%       tradeoff, the maximal distance with [3 4] weights is around 11000
%       pixels.
%
%   RES = imDistanceMap(..., 'normalize', TF);
%   If TF is true, normalizes the resulting map by the first weight.
%   Default is true.
%
%   RES = imDistanceMap(..., 'verbose', true);
%   Displays info on iterations.
%
%   Examples:
%   % Computes distance map on closed circles, with Borgefors Metric
%     img = imread('circles.png');
%     se = strel('disk', 6);
%     img2 = imclose(img, se);
%     dist = imDistanceMap(img2, [3 4]);
%     imshow(dist, []); colormap jet;
%
%   % uses the examples from bwdist with different distances
%     img = ones(255, 255);
%     img(126, 126) = 0;
%     res1 = imDistanceMap(img);
%     res2 = imDistanceMap(img, [1 inf]);
%     res3 = imDistanceMap(img, [1 1]);
%     res4 = imDistanceMap(img, [1 1.5]);
%     figure;
%     subplot(221); subimage(mat2gray(res1));
%     hold on; imcontour(res1); title('quasi-euclidean');
%     subplot(222); subimage(mat2gray(res2));
%     hold on; imcontour(res2); title('city-block');
%     subplot(223); subimage(mat2gray(res3));
%     hold on; imcontour(res3); title('chessboard');
%     subplot(224); subimage(mat2gray(res4));
%     hold on; imcontour(res4); title('approx euclidean');
%
%
%   See also:
%     bwdist, imGeodesicDistanceMap, imSeparateParticles
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2012-08-20,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2012 INRA - Cepia Software Platform.

%   HISTORY
%   2010-08-25 fix memory allocation for large images, add vebosity option
%   2012-08-20 adapt imChamferDistance to create imDistanceMap
%   2018-02-22 add management of contiguous labels, and of three weights

%% Process input arguments

% default weights for orthogonal or diagonal
weights = [5 7 11];

normalize = true;

% extract user-specified weights
if ~isempty(varargin)
    weights = varargin{1};
    varargin(1) = [];
end

% extract verbosity option
verbose = false;
if length(varargin) > 1
    varName = varargin{1};
    if ~ischar(varName)
        error('Require options as name-value pairs');
    end

    if strcmpi(varName, 'normalize')
        normalize = varargin{2};
    elseif strcmpi(varName, 'verbose')
        verbose = varargin{2};
    else
        error(['unknown option: ' varName]);
    end
end

%% Initialisations

% determines type of output from type of weights
outputType = class(weights);

% small check up to avoid degenerate cases
w1 = weights(1);
w2 = weights(2);
if w2 < w1
    w2 = 2 * w1;
end

% shifts in directions i and j for (1) forward and (2) backward iterations
if length(weights) == 2
    nShifts = 4;
    di1 = [-1 -1 -1  0];
    dj1 = [-1  0  1 -1];
    di2 = [+1 +1 +1  0];
    dj2 = [-1  0  1 +1];
    ws =  [w2 w1 w2 w1];

elseif length(weights) == 3
    nShifts = 8;
    w3 = weights(3);
    di1 = [-2 -2 -1 -1 -1 -1 -1  0];
    dj1 = [-1 +1 -2 -1  0  1 +2 -1];
    di2 = [+2 +2 +1 +1 +1 +1 +1  0];
    dj2 = [-1 +1 +2 +1  0 -1 -2 +1];
    ws =  [w3 w3 w3 w2 w1 w2 w3 w1];
end

% allocate memory for result
dist = ones(size(img), outputType);

% init result: either max value, or 0 for marker pixels
if isinteger(w1)
    dist(:) = intmax(outputType);
else
    dist(:) = inf;
end
dist(~img) = 0;

% size of image
[D1, D2] = size(img);

%% Forward iteration

if verbose
    disp('Forward iteration %d');
end

for i = 1:D1
    for j = 1:D2
        % computes only for pixels within a region
        if img(i, j) == 0
            continue;
        end

        % compute minimal propagated distance
        newVal = dist(i, j);
        for k = 1:nShifts
            % coordinate of neighbor
            i2 = i + di1(k);
            j2 = j + dj1(k);

            % check bounds
            if i2 < 1 || i2 > D1 || j2 < 1 || j2 > D2
                continue;
            end

            % compute new value
            if img(i2, j2) == img(i, j)
                % neighbor in same region
                % -> add offset weight to neighbor distance
                newVal = min(newVal, dist(i2, j2) + ws(k));
            else
                % neighbor in another region
                % -> initialize with the offset weight
                newVal = min(newVal, ws(k));
            end

        end

        % if distance was changed, update result
        dist(i,j) = newVal;
    end

end % iteration on lines

%% Backward iteration

if verbose
    disp('Backward iteration');
end

for i = D1:-1:1
    for j = D2:-1:1
        % computes only for foreground pixels
        if img(i, j) == 0
            continue;
        end

        % compute minimal propagated distance
        newVal = dist(i, j);
        for k = 1:nShifts
            % coordinate of neighbor
            i2 = i + di2(k);
            j2 = j + dj2(k);

            % check bounds
            if i2 < 1 || i2 > D1 || j2 < 1 || j2 > D2
                continue;
            end

            % compute new value
            if img(i2, j2) == img(i, j)
                % neighbor in same region
                % -> add offset weight to neighbor distance
                newVal = min(newVal, dist(i2, j2) + ws(k));
            else
                % neighbor in another region
                % -> initialize with the offset weight
                newVal = min(newVal, ws(k));
            end

        end

        % if distance was changed, update result
        dist(i,j) = newVal;
    end

end % line iteration

if normalize
    dist(dist>0) = dist(dist>0) / w1;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function range = p15_range ( )

%*****************************************************************************80
%
%% P15_RANGE returns an interval bounding the root for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = - 10.0;
  range(2) =   10.0;

  return
end


text:
function reducedisp_wfmeas(varargin);

% REDUCEDISP_WFMEAS(STARTTIME,ENDTIME,SCNL,DIST,DS,DBOUT,FILT,TIMESTEP,ALGORITHM)
% This function reads calibrated waveforms and write out a wfmeas
% database table containing reduced displacement values. The reduced displacment
% functionality requres the Antelope toolbox, the waveform object toolbox and the
% filterobject toolbox.
%
% Example:
%    starttime = '12/01/2006 19:40:00';	% formatted start time (string)
%    endtime = '12/31/2006 24:00:00';	% formatted end data (string)
%    scnl = scnlobject('BEZB','HHZ')
%    dist = 5.224;			% distance to assumed source (km)
%    ds = datasource('antelope','/home/admin/databases/PIRE/wf/pire_2006');
%    dbout = 'databases/BEZB_p5-10Hz';   	% output db (string)
%    filt = filterobject('B',[0.5 10],4);	% filter to apply to traces
%    Tstep = 600;				% time step and window wdith
%    algorithm = 'BODY';    % algorithm (BODY or SURF)
%    reducedisp_wfmeas(starttime,endtime,scnl,dist,ds,dbout,filt,Tstep,algorithm);

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if numel(varargin) ~= 9
    error('Incorrect number of input arguments');
else
    T1 = varargin{1};
    T2 = varargin{2};
    scnl = varargin{3};
    distkm = varargin{4};
    ds = varargin{5};
    dbout = varargin{6};
    flt = varargin{7};
    Tstep = varargin{8};
    algorithm = varargin{9};
end

% PREP TIME PARAMETERS
Tstep = Tstep/86400;	% time step between measurements in seconds
Twin = Tstep;           % could be used later to have different time steps and windows.
TWstep = 6*3600/86400;	% load 6 hr. waveforms at a time
gapThresh = 0.5;        % [<1.0] maximum allowable gap size (ratio to data length)

% GET WAVEFORM TIMES (LARGE CHUCKS OF TIME)
% round *up* the waveform start times the nearest Tstep
tmod = mod(datenum(T1),Tstep);
if tmod == 0
	tmod = Tstep;
end;
t1 = datenum(T1) + ( Tstep -  tmod);
Nwf = floor( (datenum(T2)-t1(1)) / TWstep );
t1 = t1(1) + TWstep*[0:Nwf-1];
t2 = t1 + TWstep;

for i = 1:Nwf

	% LOAD WAVEFORM AND TEST FOR GOOD DATA
    	PROCESS = 1;
	try
		W = waveform(ds,scnl,t1(i),t2(i));
		if isempty(W)
        		disp(['Skipping ' datestr(t1(i),31) ' to ' datestr(t2(i),31) ' . Empty waveform.']);
			PROCESS = 0;
    		end
	catch
		disp(['Skipping ' datestr(t1(i),31) ' to ' datestr(t2(i),31) ' . Problem loading waveform.']);
    		PROCESS = 0;
	end
    	if PROCESS
		if (get(W,'DURATION_EPOCH') < gapThresh*Tstep)
        		disp(['Skipping ' datestr(t1(i),31) ' to ' datestr(t2(i),31) '. Data is more than ' num2str(100*gapThresh) '% gaps.']);
    			PROCESS = 0;
		end
	end

if PROCESS

        % FILTER DATA
        % TODO: should be changed to be gap aware
        W = fillgaps(W,'meanAll');
        W = demean(detrend(W));
		W = filtfilt(flt,W);

        % CREATE SUBSETS OF DATA
        Nwin = get(W,'DURATION') / Twin;
        To = get(W,'START_MATLAB');
        for ii = 1:Nwin
            Tbegin(ii) = To+(ii-1)*Tstep;
            Tend(ii) = To+(ii-1)*Tstep+Twin;
        end;
        w = extract(W,'TIME',Tbegin,Tend);

        % GET REDUCED DISPLACEMENT
        disp(['Calculating reduced displacements on ' datestr(t1(i),31) ' to ' datestr(t2(i),31) ' ...']);
        w = reducedisp_calc(w,distkm,algorithm);

        % WRITE DB TABLE
        disp(['Writing reduced displacements to ' dbout ' ...']);
        reducedisp_write_wfmeas(dbout,w,flt,algorithm);

    end;

end;

author: geoscience-community-codes

repo: GISMO


text:
function [data] = ft_redefinetrial(cfg, data)

% FT_REDEFINETRIAL allows you to adjust the time axis of your data, i.e. to
% change from stimulus-locked to response-locked. Furthermore, it allows
% you to select a time window of interest, or to resegment your long trials
% into shorter fragments.
%
% Use as
%   [data] = ft_redefinetrial(cfg, data)
% where the input data should correspond to the output of FT_PREPROCESSING and the
% configuration should be specified as explained below. Note that some options are
% mutually exclusive. If you want to use both,  you neew two calls to this function
% to avoid confusion about the order in which they are applied.
%
% For selecting a subset of trials you can specify
%   cfg.trials    = 'all' or a selection given as a 1xN vector (default = 'all')
%
% For selecting trials with a minimum length you can specify
%   cfg.minlength = length in seconds, can be 'maxperlen' (default = [])
%
% For realiging the time axes of all trials to a new reference time
% point (i.e. change the definition for t=0) you can use the following
% configuration option
%   cfg.offset    = single number or Nx1 vector, expressed in samples relative to current t=0
%
% For selecting a specific subsection of (i.e. cut out a time window
% of interest) you can select a time window in seconds that is common
% in all trials
%   cfg.toilim    = [tmin tmax] to specify a latency window in seconds, can be Nx2 vector
%
% Alternatively you can specify the begin and end sample in each trial
%   cfg.begsample = single number or Nx1 vector, expressed in samples relative to the start of the input trial
%   cfg.endsample = single number or Nx1 vector, expressed in samples relative to the start of the input trial
%
% Alternatively you can specify a new trial definition, expressed in
% samples relative to the original recording
%   cfg.trl       = Nx3 matrix with the trial definition, see FT_DEFINETRIAL
%
% Alternatively you can specify the data to be cut into (non-)overlapping
% segments, starting from the beginning of each trial. This may lead to loss
% of data at the end of the trials
%   cfg.length    = number (in seconds) that specifies the length of the required snippets
%   cfg.overlap   = number between 0 and 1 (exclusive) specifying the fraction of overlap between snippets (0 = no overlap)
%
% Alternatively you can merge or stitch pseudo-continuous segmented data back into a
% continuous representation. This requires that the data has a valid sampleinfo field
% and that there are no jumps in the signal in subsequent trials (e.g. due to
% filtering or demeaning). If there are missing segments (e.g. due to artifact
% rejection), the output data will have one trial for each section where the data is
% continuous.
%   cfg.continuous = 'yes'
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_DEFINETRIAL, FT_RECODEEVENT, FT_PREPROCESSING

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% store original datatype
dtype = ft_datatype(data);

% deal with the special case of timelock rpt_chan_time with 1 trial
oneRptTimelock = (strcmp(dtype, 'timelock') && ...
  strcmp(data.dimord, 'rpt_chan_time') && ...
  size(data.trial, 1) == 1);

% check if the input data is valid for this function, this will convert it to raw if needed
data = ft_checkdata(data, 'datatype', {'raw+comp', 'raw'}, 'feedback', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.offset       = ft_getopt(cfg, 'offset',     []);
cfg.toilim       = ft_getopt(cfg, 'toilim',     []);
cfg.begsample    = ft_getopt(cfg, 'begsample',  []);
cfg.endsample    = ft_getopt(cfg, 'endsample',  []);
cfg.minlength    = ft_getopt(cfg, 'minlength',  []);
cfg.trials       = ft_getopt(cfg, 'trials',     'all', 1);
cfg.feedback     = ft_getopt(cfg, 'feedback',   'yes');
cfg.trl          = ft_getopt(cfg, 'trl',        []);
cfg.length       = ft_getopt(cfg, 'length',     []);
cfg.overlap      = ft_getopt(cfg, 'overlap',    0);
cfg.continuous   = ft_getopt(cfg, 'continuous', 'no');

% select trials of interest
if ~strcmp(cfg.trials, 'all')
  if islogical(cfg.trials)
    ft_info('selecting %d trials\n', sum(cfg.trials));
  else
    ft_info('selecting %d trials\n', length(cfg.trials));
  end

  % select trials of interest
  tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  data   = ft_selectdata(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  if length(cfg.offset)>1 && length(cfg.offset)~=length(cfg.trials)
    cfg.offset = cfg.offset(cfg.trials);
  end
  if length(cfg.begsample)>1 && length(cfg.begsample)~=length(cfg.trials)
    cfg.begsample = cfg.begsample(cfg.trials);
  end
  if length(cfg.endsample)>1 && length(cfg.endsample)~=length(cfg.trials)
    cfg.endsample = cfg.endsample(cfg.trials);
  end
end

Ntrial = numel(data.trial);

% check the input arguments, only one method for processing is allowed
numoptions = ~isempty(cfg.toilim) + ~isempty(cfg.offset) + (~isempty(cfg.begsample) || ~isempty(cfg.endsample)) + ~isempty(cfg.trl) + ~isempty(cfg.length) + istrue(cfg.continuous);
if numoptions>1
  ft_error('you should specify only one of the options for redefining the data segments');
end
if numoptions==0 && isempty(cfg.minlength) && strcmp(cfg.trials, 'all')
  ft_error('you should specify at least one configuration option');
end

% start processing
if ~isempty(cfg.toilim)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % select a latency window from each trial
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if numel(cfg.toilim) == 2
    % specified as single [tstart tend] vector
    % expand into Ntrial X 2
    cfg.toilim = repmat(cfg.toilim(:)', Ntrial, 1);
  end

  begsample = zeros(Ntrial,1);
  endsample = zeros(Ntrial,1);
  skiptrial = false(Ntrial,1);
  for i=1:Ntrial
    if cfg.toilim(i,1)>data.time{i}(end) || cfg.toilim(i,2)<data.time{i}(1)
      begsample(i) = nan;
      endsample(i) = nan;
      skiptrial(i) = true;
    else
      begsample(i) = nearest(data.time{i}, cfg.toilim(i,1));
      endsample(i) = nearest(data.time{i}, cfg.toilim(i,2));
      data.trial{i} = data.trial{i}(:, begsample(i):endsample(i));
      data.time{i}  = data.time{i} (   begsample(i):endsample(i));
    end
  end

  % also correct the sample information
  if isfield(data, 'sampleinfo')
    data.sampleinfo(:, 1) = data.sampleinfo(:, 1) + begsample - 1;
    data.sampleinfo(:, 2) = data.sampleinfo(:, 1) + endsample - begsample;
  end

  data.time     = data.time(~skiptrial);
  data.trial    = data.trial(~skiptrial);
  if isfield(data, 'sampleinfo'),  data.sampleinfo  = data.sampleinfo(~skiptrial, :); end
  if isfield(data, 'trialinfo'),   data.trialinfo   = data.trialinfo(~skiptrial, :);  end
  ft_info('removing %d trials in which no data was selected\n', sum(skiptrial));

elseif ~isempty(cfg.offset)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % shift the time axis from each trial
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  offset = cfg.offset(:);
  offset = round(offset); % this is in samples and hence it must be expressed as integers
  if length(cfg.offset)==1
    offset = repmat(offset, Ntrial, 1);
  end
  for i=1:Ntrial
    data.time{i} = data.time{i} + offset(i)/data.fsample;
  end

elseif ~isempty(cfg.begsample) || ~isempty(cfg.endsample)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % select a latency window from each trial based on begin and/or end sample
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  begsample = cfg.begsample(:);
  endsample = cfg.endsample(:);
  if length(begsample)==1
    begsample = repmat(begsample, Ntrial, 1);
  end
  if length(endsample)==1
    endsample = repmat(endsample, Ntrial, 1);
  end
  for i=1:Ntrial
    data.trial{i} = data.trial{i}(:, begsample(i):endsample(i));
    data.time{i}  = data.time{i} (   begsample(i):endsample(i));
  end

  % also correct the sampleinfo
  if isfield(data, 'sampleinfo')
    sampleinfo = data.sampleinfo(:, 1);
    data.sampleinfo(:, 1) = sampleinfo(:, 1) + begsample - 1;
    data.sampleinfo(:, 2) = sampleinfo(:, 1) + endsample - 1;
  end

elseif ~isempty(cfg.trl)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % select new trials from the existing data
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ischar(cfg.trl)
    % load the trial information from file
    newtrl = loadvar(cfg.trl, 'trl');
  else
    newtrl = cfg.trl;
  end

  % ensure that sampleinfo is present, otherwise ft_fetch_data will crash
  data = ft_checkdata(data, 'hassampleinfo', 'yes');

  % make a copy of the old data
  dataold = data;

  % make the header
  hdr = ft_fetch_header(dataold);

  % start with a completely new data structure
  data          = keepfields(dataold, {'cfg' 'fsample' 'label' 'topo' 'topolabel' 'unmixing' 'mixing' 'grad' 'elec' 'opto'}); % account for all potential fields to be copied over
  data.hdr      = hdr;
  data.trial    = cell(1,size(newtrl,1));
  data.time     = cell(1,size(newtrl,1));

  if isfield(dataold, 'trialinfo')
    ft_warning('Original data has trialinfo, using user-specified trialinfo instead');
  end

  if ~istable(newtrl)
    begsample = newtrl(:,1);
    endsample = newtrl(:,2);
    offset    = newtrl(:,3);
  else
    begsample = newtrl.begsample;
    endsample = newtrl.endsample;
    offset    = newtrl.offset;
  end
  trllength = endsample - begsample + 1;

  for iTrl=1:size(newtrl, 1)

    data.trial{iTrl} = ft_fetch_data(dataold, 'header', hdr, 'begsample', begsample(iTrl), 'endsample', endsample(iTrl), 'chanindx', 1:hdr.nChans, 'skipcheckdata', 1);
    data.time{iTrl}  = offset2time(offset(iTrl), dataold.fsample, trllength(iTrl));

    % The following ensures correct handling of trialinfo.

    % Determine which old trials are present in new trials
    iTrlorig = find(begsample(iTrl) <= dataold.sampleinfo(:,2) & endsample(iTrl) >= dataold.sampleinfo(:,1));

    if size(newtrl,2)>3 % In case user specified additional trialinfo
      data.trialinfo(iTrl,:) = newtrl(iTrl,4:end);
    elseif isfield(dataold,'trialinfo') % If old data has trialinfo
      if (numel(iTrlorig) == 1 ...      % only 1 old trial to copy trialinfo from, or
          || size(unique(dataold.trialinfo(iTrlorig,:),'rows'),1)) ... % all old trialinfo rows are identical
          && ~any(diff(dataold.sampleinfo(:,1))<=0) % and the trials are consecutive segments
        data.trialinfo(iTrl,:) = dataold.trialinfo(iTrlorig(1),:);
      else
        ft_error('Old trialinfo cannot be combined into new trialinfo, please specify trialinfo in cfg.trl(:,4)');
      end
    end
  end % for iTrl

  % adjust the sampleinfo in the output
  if isfield(dataold, 'sampleinfo')
    % adjust the sample information
    data.sampleinfo  = [begsample endsample];
  end

elseif ~isempty(cfg.length)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % cut the existing trials into segments of the specified length
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  data = ft_checkdata(data, 'hassampleinfo', 'yes');

  % create dummy trl-matrix and recursively call ft_redefinetrial
  nsmp    = round(cfg.length*data.fsample);
  nshift  = round((1-cfg.overlap)*nsmp);

  if nshift<=0
    ft_error('the overlap is too large');
  end

  newtrl = zeros(0,4);
  for k = 1:numel(data.trial)
    begsample = data.sampleinfo(k,1);
    endsample = data.sampleinfo(k,2);
    offset    = time2offset(data.time{k}, data.fsample);
    thistrl   = (begsample:nshift:(endsample+1-nsmp))';
    if ~isempty(thistrl) % the trial might be too short
      thistrl(:,2) = thistrl(:,1) + nsmp - 1;
      thistrl(:,3) = thistrl(:,1) + offset - thistrl(1,1);
      thistrl(:,4) = k; % keep the trial number in the 4th column, this is needed further down
      newtrl = cat(1, newtrl, thistrl);
    end
  end
  clear begsample endsample offset

  tmpcfg = keepfields(cfg, {'feedback', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.trl = newtrl;

  if isfield(data, 'trialinfo') && ~istable(data.trialinfo)
    % replace the trial number with the original trial information
    tmpcfg.trl = [newtrl(:,1:3) data.trialinfo(newtrl(:,4),:)];
  elseif isfield(data, 'trialinfo') && istable(data.trialinfo)
    % construct the trl matrix as a table
    begsample = newtrl(:,1);
    endsample = newtrl(:,2);
    offset    = newtrl(:,3);
    tmpcfg.trl = [table(begsample, endsample, offset) data.trialinfo(newtrl(:,4),:)];
  elseif ~isfield(data, 'trialinfo')
    % discard the trial number
    tmpcfg.trl = newtrl(:,1:3);
  end

  data   = removefields(data, {'trialinfo'}); % these are in the additional columns of tmpcfg.trl
  data   = ft_redefinetrial(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

elseif istrue(cfg.continuous)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % identify consecutive segments that can be glued back together
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  data = ft_checkdata(data, 'hassampleinfo', 'yes');

  boolvec = artifact2boolvec(data.sampleinfo);
  newtrl = boolvec2trl(boolvec);

  % In general: An offset of 0 means that the first sample of the trial corresponds
  % to the trigger. A positive offset indicates that the first sample is later than
  % the trigger.

  % here we want to use the start of the recording as t=0
  newtrl(:,3) = newtrl(:,1) - 1;

  tmpcfg = keepfields(cfg, {'feedback', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.trl = newtrl;

  data   = removefields(data, {'trialinfo'}); % the trialinfo does not apply any more
  data   = ft_redefinetrial(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

end % processing the realignment or data selection

if ~isempty(cfg.minlength)
  Ntrial    = length(data.trial);
  trllength = zeros(Ntrial, 1);
  % determine the length of each trial
  for i=1:Ntrial
    trllength(i) = size(data.trial{i},2) * 1/data.fsample; % this the the DURATION of the selected samples
  end
  if ischar(cfg.minlength) && strcmp(cfg.minlength, 'maxperlen')
    minlength = max(trllength);
  else
    minlength = cfg.minlength;
  end
  % remove trials that are too short
  skiptrial = (trllength<minlength);
  %if ~isempty(trl), trl = trl(~skiptrial,:); end
  data.time  = data.time(~skiptrial);
  data.trial = data.trial(~skiptrial);
  if isfield(data, 'sampleinfo'), data.sampleinfo  = data.sampleinfo(~skiptrial, :); end
  if isfield(data, 'trialinfo'),  data.trialinfo   = data.trialinfo (~skiptrial, :); end
  ft_info('removing %d trials that are too short\n', sum(skiptrial));
end

% convert back to input type if necessary
switch dtype
  case 'timelock'
    data = ft_checkdata(data, 'datatype', 'timelock');
    if oneRptTimelock
      % deal with the special case of rpt_chan_time timelock data with one
      % repetition
      data.trial = reshape(data.avg, [1 size(data.avg)]);
      data.dimord = 'rpt_chan_time';
    end
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [keypointsAll, imgidxs_missing] = det2keypointsImgMulticut(expidx,parts,annolist_gt,bAllPartsOn)

fprintf('det2keypointsImgMulticut()\n');

if (nargin < 4)
    bAllPartsOn = false;
end

p = rcnn_exp_params(expidx);

predDir = p.multicutDir;
%{
if (isfield(p,'predDir'))
    predDir = [p.predDir '/multicut/'];
else
    conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', [p.expDir '/' p.shortName]);
    predDir = [conf.cache_dir '/multicut/'];
end
%}
keypointsAll = repmat(struct('imgname','','det',nan(16,3)), length(annolist_gt), 1);

if (isfield(p,'oracleClusEval'))
    oracleClusEval = p.oracleClusEval;
else
    oracleClusEval = false;
end

if (isfield(p,'oracleEval'))
    oracleEval = p.oracleEval;
else
    oracleEval = false;
end

if (isfield(p,'oracleSymmEval'))
    oracleSymmEval = p.oracleSymmEval;
else
    oracleSymmEval = false;
end

sim_class_id = [
    6
    5
    4
    3
    2
    1
    12
    11
    10
    9
    8
    7
    13
    14
];

cidxs = p.cidxs;

nmissing = 0;
imgidxs_missing = [];
for imgidx = 1:length(annolist_gt)

    fprintf('.');

    keypointsAll(imgidx).imgname = annolist_gt(imgidx).image.name;

    fname = [predDir '/imgidx_' padZeros(num2str(imgidx),4) '_cidx_' num2str(cidxs(1)) '_' num2str(cidxs(end))];
    try
        load(fname,'unLab','unPos','unProb');
    catch
        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            jidx = parts(pidx+1).pos(1);
            keypointsAll(imgidx).det(jidx+1,:) = [inf inf -1];
        end
        nmissing = nmissing + 1;
        imgidxs_missing = [imgidxs_missing; imgidx];
        continue;
    end

    img = imread(annolist_gt(imgidx).image.name);
    [Y,X,~] = size(img);
    Y = 0.5*Y; X = 0.5*X;

%     X = annolist_gt(imgidx).annorect.objpos.x;
%     Y = annolist_gt(imgidx).annorect.objpos.y;

    rect = annolist_gt(imgidx).annorect;
    refDist = util_get_torso_size(rect);
    for i = 1:length(p.pidxs)
        if (~ismember(i,cidxs))
            continue;
        end
        pidx = p.pidxs(i);
        assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
        jidx = parts(pidx+1).pos(1);

%         detidxs = find(unLab(:,1) == i-1);
        detidxs = find(unLab(:,1) == find(cidxs == i)-1);
        if (isempty(detidxs))
            if (bAllPartsOn)
%                 [val,idx] = max(unProb(:,i));
                [val,idx] = max(unProb(:,find(cidxs == i)));
                det = unPos(idx,:);
            else
                det = [inf inf];
            end
        else
            clusidxs = unLab(detidxs,2);
            clusidxsUniq = unique(clusidxs);
            nExClus = zeros(length(clusidxsUniq),1);
            clusMean = zeros(length(clusidxsUniq),2);
            for j = 1:length(clusidxsUniq)
                nExClus(j) = length(find(unLab(:,2) == clusidxsUniq(j)));
                clusMean(j,:) = mean(unPos(unLab(:,2) == clusidxsUniq(j),:),1);
            end
            d = sqrt(sum((clusMean - repmat([X Y],size(clusMean,1),1)).^2,2));
            [val,id] = sort(d,'ascend');
%             detidxs = find(unLab(:,1) == i-1 & unLab(:,2) == clusidxsUniq(id(1)));
            detidxs = find(unLab(:,1) == find(cidxs == i)-1 & unLab(:,2) == clusidxsUniq(id(1)));

%             mean between the two closest
%             d = inf(length(detidxs),length(detidxs));
%             for j = 1:length(detidxs)
%                 d(j,:) = sqrt(sum((repmat(unPos(detidxs(j),:),length(detidxs),1) - unPos(detidxs,:)).^2,2));
%             end
%             for j = 1:length(detidxs)
%                 d(j,j) = inf;
%             end
%             idx = find(d==min(min(d)));
%             [rx,cx] = ind2sub(size(d),idx);
%             det = mean(unPos(detidxs([rx(1) cx(1)]),:),1);

%             weighted sum excluding outliers
%             if (length(detidxs) > 2)
%                 d = inf(length(detidxs),length(detidxs));
%                 for j = 1:length(detidxs)
%                     d(j,:) = sqrt(sum((repmat(unPos(detidxs(j),:),length(detidxs),1) - unPos(detidxs,:)).^2,2));
%                 end
%                 dd = sum(d,2);
%                 [val,id2] = sort(dd,'descend');
%                 detidxs = detidxs(2:end);
%             end
%             w = unProb(detidxs,i);
%             w = w./sum(w);
%             det = sum(unPos(detidxs,:).*[w w],1);

%             medoid
%             m = mean(unPos(detidxs,:),1);
%             d = sqrt(sum((repmat(m,length(detidxs),1) - unPos(detidxs,:)).^2,2));
%             [val,id2] = min(d);
%             det = unPos(detidxs(id2),:);

%             maximum
%             [val,idx] = max(unProb(detidxs));
%             det = unPos(detidxs(idx),:);

%             weighted sum
%             w = unProb(detidxs,i);
            w = unProb(detidxs,find(cidxs == i));
            w = w./sum(w);
            det = sum(unPos(detidxs,:).*[w w],1);

            if (oracleEval || oracleClusEval || oracleSymmEval)
                if (oracleEval)
                    detidxs = 1:size(unPos,1);
                elseif (oracleSymmEval)
%                     detidxs = [detidxs; find(unLab(:,1) == sim_class_id(i)-1)];
                    detidxs = [detidxs; find(unLab(:,1) == sim_class_id(find(cidxs == i))-1)];
                end
                pp = util_get_annopoint_by_id(rect.annopoints.point, jidx);
                d = sqrt(sum((repmat([pp.x pp.y],size(unPos(detidxs,:),1),1) - unPos(detidxs,:)).^2,2))./refDist;
                [val,id] = min(d);
                det = unPos(detidxs(id),:);
            end

            % DEBUG
%             [val,id] = max(unProb(detidxs,i));
%             det = unPos(detidxs(id),:);
%             figure(100); clf; imagesc(img); axis equal; hold on;
%             detidxsAll = unLab(:,2) == clusidxsUniq(id(1));
%             plot(unPos(detidxsAll,1),unPos(detidxsAll,2),'r+','Markersize',10);
%             plot(unPos(detidxs,1),unPos(detidxs,2),'b+','Markersize',10);
%             plot(det(:,1),det(:,2),'y+','Markersize',10);
        end

        % part is a joint
        keypointsAll(imgidx).det(jidx+1,:) = [det,-1];
    end
    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(keypointsAll));
    end
end
fprintf(' done\n');
fprintf('nmissing: %d\n',nmissing);

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [labels,classes] = ml_string2label(strings)
%ML_STRING2LABEL Converts a cell array of strings to a num array of labels.
%
%   input -----------------------------------------------------------------
%
%       o strings : cell(N x 1), cell array of string labels
%
%   output ----------------------------------------------------------------
%
%       o labels :  (N x 1), numerical array of class labels
%
%       o classes: cell(num_classes x 1), array of unique class names
%

classes = unique(strings);
labels  = zeros(size(strings,1),1);

for i=1:size(classes,1)

  idx    = ismember(strings,classes{i});
  labels(idx) = ones(sum(idx(:)),1) .* i;

end

end

author: epfl-lasa


text:
function y = pcone(z,x,y,alpha)
%PCONE Defines a power cone x^alpha y ^(1-alpha) > |z|
%
% Input
%    z,y,x   : sclar SDPVAR objects.
%    alpha   : scalar double 0<=alpha<=1
%
% Example
%    F = pcone(z,x,y,alpha)
%
% An alternative syntax with only one argument is also possible
%    F = pcone(z)
%
%
% See also CONE

if numel(z)>1
    error('x must be a scalar')
end
if numel(x)>1
    error('x must be a scalar')
end
if numel(y)>1
    error('y must be a scalar')
end

try
    y = [x;y;z;alpha];
    y.typeflag = 20;
    y = lmi(y);
catch
    rethrow(lasterror)
end

author: shamilmamedov


text:
function par_est
% parameter estimation with derivatives              Holzbecher September 2005
% for exponential fit of exponent lambda for c(0)=c0
global tfit cfit c0

% specify fitting data
tfit =  [0.25 1 2 4 8];
cfit = [0.7716 0.5791 0.4002 0.1860 0.1019];
c0 = .8;

lambda = fzero(@myfun,0.11);
normc = norm(cfit - c0*exp(-lambda*tfit));
display (['Best fit for lambda = ' num2str(lambda)]);
display (['Norm of residuals =' num2str(normc)]);
tmax = tfit(size(tfit,2));
t = [0:0.01*tmax:tmax];
figure; plot (tfit,cfit,'or',t,c0*exp(-lambda*t),'-');
legend ('given','modelled');
text(0.5*tmax,c0*0.7,['\lambda: ' num2str(lambda)]);
text(0.5*tmax,c0*0.6,['norm of residuals: ' num2str(normc)]);
xlabel ('time'); ylabel ('concentration')

function f = myfun(lambda)
global tfit cfit c0

t = tfit;
c = c0*exp(-lambda*t);    % solve linear decay equation for c with c(0)=c0
clambda = -c.*t;          % equation for dc/dlambda
f = (cfit-c)*clambda';    % specify function f to vanish

author: Sable

repo: mcbench-benchmarks


text:
classdef Tooltipable < handle
    % Mixin for component with Tooltip property

    % Copyright 2020-2021 The MathWorks Inc.

    %% Properties
    properties (AbortSet)

        % Tooltip of the component
        Tooltip (1,1) string

    end %properties

    %% Internal properties
    properties (AbortSet, Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % List of graphics controls to apply to
        TooltipableComponents (:,1) matlab.graphics.Graphics

    end %properties

    %% Accessors
    methods

        function set.Tooltip(obj,value)
            obj.Tooltip = value;
            obj.updateTooltipableComponents()
        end

        function set.TooltipableComponents(obj,value)
            obj.TooltipableComponents = value;
            obj.updateTooltipableComponents()
        end

    end %methods

    %% Methods
    methods (Access = protected)

        function updateTooltipableComponents(obj)

            wt.utility.fastSet(obj.TooltipableComponents,"Tooltip",obj.Tooltip);

        end %function

    end %methods

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%begin
% We will begin with the trophy from the CollisionModel demo, which we
% do not want to hit.
trophy = cmdl_trophy

% We will load a model of the SerialLink object with STL data included.
load('p560_col.mat');
p560

% The two are precariously close!
fig1 = figure;
qz = zeros(1,6);
p560.plot3d(qz);
hold on
trophy.plot;

% The point data from the STL data is not dense. We can only check for
% collisions where there are points. You can rectify this in a solid
% model program, such as SolidWorks, by adding a grid of split lines on
% the relevant surfaces.

% We can check for collisions
c = p560.collisions(qz, trophy)

% The false/zero value indicates the two are not colliding. If the 2nd
% joint turns a little to move the red link lower,
q = [0 -0.1 0 0 0 0]

c = p560.collisions(q, trophy)

% There is now a collision. Let's verify this
close(fig1); figure;
p560.plot3d(q);
hold on
trophy.plot;

% The robot does indeed collide with the trophy. It is important to
% note that the program sees the cup of the plot as a solid object,
% even though we have plotted is as an open, concave vessel. If the arm
% was to be in the void of the cup but not actually touching, it would
% still be considered colliding.

author: Allopart

repo: rbpf-gmapping


text:
function [ y, ierror ] = y_check_gregorian ( y )

%*****************************************************************************80
%
%% Y_CHECK_GREGORIAN checks a Gregorian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, the year, which must not be 0.
%
%    Output, integer IERROR, is 0 if Y is legal, and 1 otherwise.
%
  if ( y ~= 0 )
    ierror = 0;
  else
    ierror = 1;
  end

  return
end


text:
function ex5bvp
%EX5BVP  Example 5 of the BVP tutorial.
%   Falkner-Skan BVPs are discussed in T. Cebeci and H.B. Keller,
%   Shooting and parallel shooting methods for solving the Falkner-Skan
%   boundary-layer equation, J. Comp. Phy., 7 (1971) 289-300.  This is
%   the positive wall shear case for which the parameter beta is known
%   and the problem is to be solved for a range of the parameter.  This
%   is the hardest case of the table in the paper.
%
%   The problem is posed on [0 infinity).  As in the paper cited, the
%   boundary condition at infinity is imposed at a finite point, here
%   called 'infinity'.  It is best to start with a relatively small value
%   and increase it until consistent results are obtained.  A value of 6
%   appears to be satisfactory.  Starting with a "large" value is tempting,
%   but not a good tactic because the code will fail with the crude guess
%   and default tolerances used here.

% Copyright 1999, The MathWorks, Inc.

%  'infinity' is a variable to facilitate experimentation.
infinity = 6;

%  The constant guess for the solution satisfies the boundary conditions.
solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);

options = bvpset('stats','on');

sol = bvp4c(@ex5ode,@ex5bc,solinit,options);
eta = sol.x;
f = sol.y;

fprintf('\n');
fprintf('Cebeci & Keller report f''''(0) = 0.92768.\n')
fprintf('Value computed here is f''''(0) = %7.5f.\n',f(3,1))

clf reset
plot(eta,f(2,:));
axis([0 infinity 0 1.4]);
title('Falkner-Skan equation, positive wall shear, \beta = 0.5.')
xlabel('\eta')
ylabel('df/d\eta')
shg

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

function dfdeta = ex5ode(eta,f)
%EX5ODE  ODE function for Example 5 of the BVP tutorial.
beta = 0.5;
dfdeta = [ f(2)
           f(3)
          -f(1)*f(3) - beta*(1 - f(2)^2) ];

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

function res = ex5bc(f0,finf)
%EX5BC  Boundary conditions for Example 5 of the BVP tutorial.
res = [f0(1)
       f0(2)
       finf(2) - 1];


text:
%FINDNLAB Determine indices of specified classes (numeric)
%
%    J = FINDNLAB(A,NLAB)
%
% INPUT
%     A        Dataset
%     NLAB     vector with numerical indices of reqested classes
%
% OUTPUT
%     J        vector of sample indices
%
% DESCRIPTION
% J is a column vector of indices to the objects in A that have one of the
% numeric labels stored in NLAB. Note that NLAB should contain indices to
% the labellist of A, which may be retrieved by LABLIST = GETLABLIST(A).
% This command thereby retrieves the indices to all objects a specific class
% or set of classes.

author: marianux

repo: ecg-kit

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


text:
function [M0, T1] = mtv_compute_m0_t1(data, flipAngles, TR, b1Map, roi, fixT1, verbose)

% function function [M0 T1] = mtv_compute_m0_t1 (data, flipAngles, TR [, b1Map, roi, fixT1, verbose])
% -----------------------------------------------------------
% This function performs a weighted-least squares data fit
% on SPGR T1-weighted data set
% INPUTS:
% data: width x length x slices x flipAngles matrix
% flipAngles: a vector of flip angles (in degrees) that corresponds to
% data(:,:,:,4)
% TR: in s
% b1Map: a width x length x slices matrix that contains the relative flip
% angle (i.e. if nominal alpha is 60, and measured alpha is 61, then b1Map
% = 61/60
% roi: binary mask (Matrix)
% fixT1: scalar or Matrix. If 0 --> do not fix T1.
% verbose: logical
%

if ndims(data)<3, data = permute(data(:),[2 3 4 1]); end
dataSize = size(data);
if (nargin < 4) || isempty(b1Map)
    b1Map = ones([dataSize(1:end-1) 1]);
end

if nargin<5 || isempty(roi)
    roi = ones(dataSize(1:end-1));
end

if nargin<6
    fixT1 = 0;
end

if nargin<7
    verbose = 0;
end

if max(size(b1Map) ~= dataSize(1:length(size(b1Map)))), error('B1 size is different from data size'); end

dims = size(data);
T1 = zeros(dims(1:end-1));
M0 = zeros(dims(1:end-1));

warning('off');
%sprintf('%s\n\n\n','loop over voxels...')
for vox=1:dims(1)*dims(2)*dims(3)
    %disp([sprintf('\b\b\b\b\b%3i',floor(vox/(dims(1)*dims(2)*dims(3))*100)) '%'])
    if roi(vox) && b1Map(vox)~=0

    kk=floor((vox-1)/(dims(1)*dims(2))); jj=floor((vox-kk*dims(1)*dims(2)-1)/(dims(1))); ii=vox-kk*dims(1)*dims(2)-jj*dims(1);
    kk=kk+1; jj=jj+1;
    %ii
        if ~fixT1
            %% T1 Mapping Using Variable Flip Angle SPGR Data With
            % Flip Angle Correction - Liberman, et al.
            y = squeeze(squeeze(data(ii,jj,kk, :)))./sin(flipAngles/180*pi*b1Map(vox))';
            x = squeeze(squeeze(data(ii,jj,kk, :)))./tan(flipAngles/180*pi*b1Map(vox))';

            % fit data
            param = polyfit(x,y,1);
%             fitresult = LinearFit(x, y, verbose);
%             param = coeffvalues(fitresult); % slope and intercept of the fitting
            param(isnan(param))=0;

%             ci = confint(fitresult,0.682); % confidence interval of the fitting (returns the slope and intercept of the lines framing the fit) -- corresponding to 2*sigma
%             ci(isnan(ci)) = 0;
%
            % compute PD and T1
            [~,T1(vox)]=getT1(param,TR);
        else
            %% if T1 is known
            T1(vox)=fixT1(min(vox,end));
        end
        %% Get M0
        [FA,iFA]=min(flipAngles);

        M0(vox)=getM0fromT1(T1(vox),TR(1),data(ii,jj,kk,iFA),FA*b1Map(vox));
        % add length of the confidence interval at 68.2% in the fourth
        % dimension

        if ~fixT1
%             [~,t1_min]=getT1(ci(1,:),TR);
%             [~,t1_max]=getT1(ci(2,:),TR);
%             M02(vox)=abs(getM0fromT1(t1_max,TR,data(ii,jj,kk,ialpha),alpha)-getM0fromT1(t1_min,TR,data(ii,jj,kk,ialpha),alpha));
%             T12(vox)=abs(t1_max-t1_min);
        end

    end
end
if ~fixT1
%    T1(:,:,:,2)=T12;
%     M0(:,:,:,2)=M02;
end
%display('...done')

function [fitresult, gof] = LinearFit(x, y, verbose)
%CREATEFIT(X,Y)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input : x
%      Y Output: y
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 22-Jan-2015 15:51:35

% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( x, y );

% Set up fittype and options.
ft = fittype( 'poly1' );

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft );

% Plot fit with data.
if verbose && max(xData~=0) && max(yData~=0)

    figure(100)
    h = plot( fitresult, xData, yData,'+');
    set(h,'MarkerSize',30)
    legend( h, 'y vs. x', 'untitled fit 1', 'Location', 'NorthEast' );
    p11 = predint(fitresult,x,0.95,'observation','off');
    hold on
    plot(x,p11,'m--'); drawnow;
    hold off
    % Label axes
    xlabel( 'x' );
    ylabel( 'y' );
    grid on
    saveas(gcf,['temp.jpg']);
end

function [M0,T1]=getT1(param,TR)
a=param(1); % slope
b=param(2); % intercept
if a>0
    T1 = -TR/log(a);
else  % due to noise or bad fitting
    T1 = 0.000000000000001;
end
M0 = b/(1-exp(-TR/T1));

function M0=getM0fromT1(T1,TR,S,FA)
% Volz, S., N�th, U., Deichmann, R., 2012. Correction of systematic errors in quantitative proton density mapping. Magn. Reson. Med. 68, 74?85.
% Steady state
ST=(1-exp(-TR./T1))./(1-cos(FA*pi/180).*exp(-TR./T1)).*sin(FA*pi/180);
% M0=RP*PD (RP=Receiver Profile)
M0=S/ST;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
clear; clc;
disp('compile native libraries')
mex -O -largeArrayDims sparse_inp.c
mex -O -largeArrayDims sparse_update.c
disp('done compiling.')
test_consist1
test_consist2

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR