text:
function varargout = panel_inverse_2018(varargin)
% PANEL_INVERSE_2018: Inverse modeling GUI
%
% USAGE:  bstPanel = panel_inverse_2018('CreatePanel', sProcess, sFiles)                                                 : Called from the pipeline editor
%         bstPanel = panel_inverse_2018('CreatePanel', Modalities, isShared, HeadModelType, nSamplesNoise, nSamplesData) : Called from the interactive interface

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

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(Modalities, isShared, HeadModelType, nSamplesNoise, nSamplesData) %#ok<DEFNU>
    panelName = 'InverseOptions';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    % Initializations
    strDiagRecommNoise = '';
    strDiagRecommData = '';
    NoiseMethod = [];
    DataMethod = [];

    % GUI CALL:  panel_inverse_2018('CreatePanel', Modalities, isShared, HeadModelType, nSamplesNoise, nSamplesData)
    if iscell(Modalities)
        isProcess = 0;
        isFull = 0;
        % Get default inverse options
        OPTIONS = bst_inverse_linear_2018();

        % Recommendations for noise covariance
        NoiseMethod = OPTIONS.NoiseMethod;
        if ~isempty(nSamplesNoise)
            nChannels = nnz(any(~isnan(nSamplesNoise) & (nSamplesNoise > 2)));
            nSamples  = max(nSamplesNoise(~isnan(nSamplesNoise)));
            % Check if we should force the diagonal selection
            if (nSamples < nChannels*(nChannels+1)/2)
                strDiagRecommNoise = '  (recommended)';
                NoiseMethod = 'diag';
                OPTIONS.NoiseMethod = 'diag';
            end
        end
        % Recommendations for data covariance
        DataMethod = OPTIONS.NoiseMethod;
        if ~isempty(nSamplesData)
            nChannels = nnz(any(~isnan(nSamplesData) & (nSamplesData > 2)));
            nSamples  = max(nSamplesData(~isnan(nSamplesData)));
            % Check if we should force the diagonal selection
            if (nSamples < nChannels*(nChannels+1)/2)
                strDiagRecommData = '  (recommended)';
                DataMethod = 'diag';
            end
        end

    % PROCESS CALL:  panel_inverse_2018('CreatePanel', sProcess, sFiles)
    else
        isProcess = 1;
        % Get inputs
        sProcess = Modalities;
        sFiles   = isShared;
        % Get inverse options
        OPTIONS = sProcess.options.inverse.Value;
        % List of sensors
        Modalities = {'MEG MAG', 'MEG GRAD', 'MEG', 'EEG', 'ECOG', 'SEEG'};
        if ~isempty(sFiles(1).ChannelTypes)
            Modalities = intersect(sFiles(1).ChannelTypes, Modalities);
            if any(ismember({'MEG MAG','MEG GRAD'}, Modalities))
                Modalities = setdiff(Modalities, 'MEG');
            end
        end
        % Shared kernel
        isShared = (sProcess.options.output.Value == 1);
        isFull   = (sProcess.options.output.Value == 3);
        % Head model type: Get from head model
        sStudyChan = bst_get('ChannelFile', sFiles(1).ChannelFile);
        if ~isempty(sStudyChan) && ~isempty(sStudyChan.HeadModel)
            HeadModelFile = sStudyChan.HeadModel(sStudyChan.iHeadModel).FileName;
            HeadModelMat = in_bst_headmodel(HeadModelFile, 0, 'HeadModelType');
            HeadModelType = HeadModelMat.HeadModelType;
        else
            HeadModelType = 'surface';
        end
    end
    % Initializations
    isFirstCombinationWarning = 1;

    % ==== FRAME STRUCTURE ====
    % Create main panel: split in top (interface) / left(standard options) / right (details)
    jPanelNew = gui_component('panel');
    jPanelNew.setBorder(BorderFactory.createEmptyBorder(12,12,12,12));
    % Create top panel
    jPanelTop = gui_river([1,1], [0,6,6,6]);
    jPanelNew.add(jPanelTop, BorderLayout.NORTH);
    % Create left panel
    jPanelLeft = java_create('javax.swing.JPanel');
    jPanelLeft.setLayout(GridBagLayout());
    jPanelNew.add(jPanelLeft, BorderLayout.WEST);
    % Create right panel
    jPanelRight = java_create('javax.swing.JPanel');
    %jPanelRight.setLayout(BoxLayout(jPanelRight, BoxLayout.Y_AXIS));
    jPanelRight.setLayout(GridBagLayout());
    jPanelRight.setBorder(BorderFactory.createEmptyBorder(0,12,0,0));
    jPanelNew.add(jPanelRight, BorderLayout.EAST);
    % Default constrains
    c = GridBagConstraints();
    c.fill    = GridBagConstraints.HORIZONTAL;
    c.weightx = 1;
    c.weighty = 0;

    % ==== TOP PANEL ====
    % Comment
    gui_component('label', jPanelTop, 'br', 'Comment:  ', [], [], [], []);
    jTextComment = gui_component('text', jPanelTop, 'hfill', '', [], '', [], []);
    if isProcess && isfield(OPTIONS, 'Comment') && ~isempty(OPTIONS.Comment)
        jTextComment.setText(OPTIONS.Comment);
    end

    % ==== PANEL: METHOD ====
    jPanelMethod = gui_river([1,1], [0,6,6,6], 'Method');
        jGroupMethod  = ButtonGroup();
        jRadioMethodMn  = gui_component('radio', jPanelMethod, [],   'Minimum norm imaging', jGroupMethod, '', @Method_Callback, []);
        jRadioMethodBf  = gui_component('radio', jPanelMethod, 'br', 'LCMV beamformer',      jGroupMethod, '', @Method_Callback, []);
        jRadioMethodDip = gui_component('radio', jPanelMethod, 'br', 'Dipole modeling',      jGroupMethod, '', @Method_Callback, []);
        if ~isProcess
            jRadioMethodMem = gui_component('radio', jPanelMethod, 'br', 'MEM: Max entropy on the mean', jGroupMethod, '', @Method_Callback, []);
        else
            jRadioMethodMem = [];
        end
        % Default selection
        switch lower(OPTIONS.InverseMethod)
            case 'minnorm',  jRadioMethodMn.setSelected(1);
            case 'gls',      jRadioMethodDip.setSelected(1);
            case 'lcmv',     jRadioMethodBf.setSelected(1);
            case 'mem',      disp('BST> Warning: Running MEM from a script is not handled yet.');
        end
        % Disable Beamformer if no data covariance
        if ~isProcess && isempty(nSamplesData)
            jRadioMethodBf.setEnabled(0);
        end
        % Disable Dipoles/Beamformer if mixed head models
        if ~isProcess && strcmpi(HeadModelType, 'mixed')
            jRadioMethodBf.setEnabled(0);
            jRadioMethodDip.setEnabled(0);
        end
        % Disable MEM for shared/volume
        if ~isempty(jRadioMethodMem) && ~isProcess && (~strcmpi(HeadModelType, 'surface') || isShared)
            jRadioMethodMem.setEnabled(0);
        end

    c.gridy = 1;
    jPanelLeft.add(jPanelMethod, c);

    % ==== PANEL: MEASURE MIN NORM ====
    jPanelMeasureMN = gui_river([1,1], [0,6,6,6], 'Measure');
        jGroupMnMeasure = ButtonGroup();
        jRadioMnCurrent = gui_component('radio', jPanelMeasureMN, [],   'Current density map',  jGroupMnMeasure, '', @(h,ev)UpdatePanel(1), []);
        jRadioMnDspm    = gui_component('radio', jPanelMeasureMN, 'br', 'dSPM',                 jGroupMnMeasure, '', @(h,ev)UpdatePanel(1), []);
        jButtonDspmWarning = gui_component('label', jPanelMeasureMN, 'hfill', ' Warning: unscaled values', '', '', @(h,ev)WarningDspm(), []);
        jButtonDspmWarning.setHorizontalAlignment(jButtonDspmWarning.RIGHT);
        jRadioMnSloreta = gui_component('radio', jPanelMeasureMN, 'br', 'sLORETA',              jGroupMnMeasure, '', @(h,ev)UpdatePanel(1), []);
        % Default selection
        switch lower(OPTIONS.InverseMeasure)
            case 'amplitude',    jRadioMnCurrent.setSelected(1);
            case 'dspm2018',     jRadioMnDspm.setSelected(1);
            case 'sloreta',      jRadioMnSloreta.setSelected(1);
            otherwise,           jRadioMnCurrent.setSelected(1);
        end
    c.gridy = 2;
    jPanelLeft.add(jPanelMeasureMN, c);

    % ==== PANEL: MEASURE BEAMFORMER ====
    jPanelMeasureBf = gui_river([1,1], [0,6,6,6], 'Measure');
        jGroupBfMeasure = ButtonGroup();
        jRadioMethodBfNai = gui_component('radio', jPanelMeasureBf, 'br', 'Pseudo Neural Activity Index',  jGroupBfMeasure, '', @(h,ev)UpdatePanel(), []);
        % Default selection
        jRadioMethodBfNai.setSelected(1);
    c.gridy = 2;
    jPanelLeft.add(jPanelMeasureBf, c);

    % ==== PANEL: SOURCE MODEL ====
    jPanelModel = gui_river([1,1], [0,6,6,6], 'Source model: Dipole orientations');
        jGroupModel    = ButtonGroup();
        jRadioConstr   = gui_component('radio', jPanelModel, [],   'Constrained:  Normal to cortex',    jGroupModel, '', @(h,ev)UpdatePanel(), []);
        jRadioLoose    = gui_component('radio', jPanelModel, 'br', 'Loose constraints',                jGroupModel, '', @(h,ev)UpdatePanel(), []);
        jTextLoose     = gui_component('texttime', jPanelModel, [], '', [], '', [], []);
        gui_validate_text(jTextLoose, [], [], 0:0.1:1, '', 1, OPTIONS.Loose, []);
        jRadioUnconstr = gui_component('radio', jPanelModel, 'br', 'Unconstrained', jGroupModel, '', @(h,ev)UpdatePanel(), []);
        % Default selection
        if strcmpi(HeadModelType, 'surface')
            switch lower(OPTIONS.SourceOrient{1})
                case 'fixed',    jRadioConstr.setSelected(1);
                case 'loose',    jRadioLoose.setSelected(1);
                case 'free',     jRadioUnconstr.setSelected(1);
            end
        elseif strcmpi(HeadModelType, 'volume')
            jRadioConstr.setEnabled(0);
            jRadioLoose.setEnabled(0);
            jRadioUnconstr.setSelected(1);
        elseif strcmpi(HeadModelType, 'mixed')
            jRadioConstr.setEnabled(0);
            jRadioLoose.setEnabled(0);
            jRadioUnconstr.setEnabled(0);
        end
    c.gridy = 3;
    jPanelLeft.add(jPanelModel, c);

    % ==== PANEL: DATA TYPE ====
    jPanelSensors = gui_river([1,1], [0,6,6,6], 'Sensors');
        jCheckMeg = [];
        jCheckMegGrad = [];
        jCheckMegMag = [];
        jCheckEeg = [];
        jCheckEcog = [];
        jCheckSeeg = [];
        if ismember('MEG', Modalities)
            jCheckMeg = gui_component('checkbox', jPanelSensors, '', 'MEG', [], '', @Modality_Callback, []);
            jCheckMeg.setSelected(1);
        end
        if ismember('MEG GRAD', Modalities)
            jCheckMegGrad = gui_component('checkbox', jPanelSensors, '', 'MEG GRAD', [], '', @Modality_Callback, []);
            jCheckMegGrad.setSelected(1);
        end
        if ismember('MEG MAG', Modalities)
            jCheckMegMag = gui_component('checkbox', jPanelSensors, '', 'MEG MAG', [], '', @Modality_Callback, []);
            jCheckMegMag.setSelected(1);
        end
        if ismember('EEG', Modalities)
            jCheckEeg = gui_component('checkbox', jPanelSensors, '', 'EEG', [], '', @Modality_Callback, []);
            if (length(Modalities) == 1)
                jCheckEeg.setSelected(1);
            end
        end
        if ismember('ECOG', Modalities)
            jCheckEcog = gui_component('checkbox', jPanelSensors, '', 'ECOG', [], '', @Modality_Callback, []);
            if (length(Modalities) == 1)
                jCheckEcog.setSelected(1);
            end
        end
        if ismember('SEEG', Modalities)
            jCheckSeeg = gui_component('checkbox', jPanelSensors, '', 'SEEG', [], '', @Modality_Callback, []);
            if (length(Modalities) == 1)
                jCheckSeeg.setSelected(1);
            end
        end
    c.gridy = 4;
    jPanelLeft.add(jPanelSensors, c);

    % ==== PANEL: MEM INFO ====
    jPanelMemInfo = gui_river([1,1], [0,6,6,6], 'MEM');
        gui_component('label', jPanelMemInfo, '', 'Requires the BrainEntropy plugin.
Options defined in a separate panel.', [], '', [], []);
    c.gridy = 2;
    jPanelLeft.add(jPanelMemInfo, c);
    % ======================================================================================================

    % ==== DEPTH WEIGHTING =====
    jPanelDepth = gui_river([1,1], [0,6,6,6], 'Depth weighting');
        % Use depth weighting
        jCheckDepth = gui_component('checkbox', jPanelDepth, 'br', 'Use depth weighting', [], '', @(h,ev)UpdatePanel(), []);
        jCheckDepth.setSelected(OPTIONS.UseDepth);
        % Weightexp
        jLabelWeightExp = gui_component('label', jPanelDepth, 'br', '       Order [0,1]: ', [], '', [], []);
        jTextWeightExp  = gui_component('texttime', jPanelDepth, 'tab', num2str(OPTIONS.WeightExp), [], '', [], []);
        % Weightlimit
        jLabelWeightLimit = gui_component('label', jPanelDepth, 'br', '       Maximal amount: ', [], '', [], []);
        jTextWeightLimit  = gui_component('texttime', jPanelDepth, 'tab', num2str(OPTIONS.WeightLimit), [], '', [], []);
    c.gridy = 1;
    jPanelRight.add(jPanelDepth, c);

    % ==== PANEL: NOISE COVARIANCE ====
    jPanelNoiseCov = gui_river([1,1], [0,6,6,6], 'Title');
        jGroupReg    = ButtonGroup();
        jRadioReg    = gui_component('radio',    jPanelNoiseCov, [], 'Text:', jGroupReg, '', @(h,ev)UpdatePanel(), []);
            jTextReg = gui_component('texttime', jPanelNoiseCov, [], '', [], '', [], []);
            gui_validate_text(jTextReg, [], [], 0:0.1:1, '', 1, OPTIONS.NoiseReg, []);
        jRadioMedian = gui_component('radio',    jPanelNoiseCov, 'br', 'Median eigenvalue', jGroupReg, '', @(h,ev)UpdatePanel(), []);
        jRadioDiag   = gui_component('radio',    jPanelNoiseCov, 'br', 'Text', jGroupReg, '', @(h,ev)UpdatePanel(), []);
        jRadioNoReg  = gui_component('radio',    jPanelNoiseCov, 'br', 'No covariance regularization', jGroupReg, '', @(h,ev)UpdatePanel(), []);
        jRadioShrink = gui_component('radio',    jPanelNoiseCov, 'br', 'Automatic shrinkage', jGroupReg, '', @(h,ev)UpdatePanel(), []);
        % Default selection
        switch lower(OPTIONS.NoiseMethod)
            case 'reg',    jRadioReg.setSelected(1);
            case 'diag',   jRadioDiag.setSelected(1);
            case 'none',   jRadioNoReg.setSelected(1);
            case 'shrink', jRadioShrink.setSelected(1);
            case 'median', jRadioMedian.setSelected(1);
        end
    c.gridy = 2;
    jPanelRight.add(jPanelNoiseCov, c);

    % ==== PANEL: SNR ====
    jPanelSnr = gui_river([1,1], [0,6,6,6], 'Regularization parameter: 1 / λ');
        jGroupSnr  = ButtonGroup();
        % Fixed SNR
        jRadioSnrFix = gui_component('radio', jPanelSnr, [], 'Signal-to-noise ratio: ', jGroupSnr, '', @(h,ev)UpdatePanel(), []);
        jTextSnrFix  = gui_component('texttime', jPanelSnr, [], '', [], '', [], []);
        gui_validate_text(jTextSnrFix, [], [], {0, 10000, 100}, '', 2, OPTIONS.SnrFixed, []);
        % Maximum source amplitude
        jRadioSnrRms = gui_component('radio', jPanelSnr, 'br', 'RMS source amplitude: ', jGroupSnr, '', @(h,ev)UpdatePanel(), []);
        jTextSnrRms  = gui_component('texttime', jPanelSnr, [], '', [], '', [], []);
        gui_validate_text(jTextSnrRms, [], [], {0, 1000000, 100}, 'scalar', 2, OPTIONS.SnrRms, []);
        gui_component('label', jPanelSnr, [], 'nAm', [], '', [], []);
        % Default selection
        switch lower(OPTIONS.SnrMethod)
            case 'rms',    jRadioSnrRms.setSelected(1);
            case 'fixed',  jRadioSnrFix.setSelected(1);
        end
    c.gridy = 3;
    jPanelRight.add(jPanelSnr, c);

    % ===== PANEL: OUTPUT MODE =====
    jPanelOutput = gui_river([1,1], [0,6,6,6], 'Output mode');
        jGroupOutput = ButtonGroup();
        jRadioKernel = gui_component('radio', jPanelOutput, [],   'Inverse kernel only', jGroupOutput, 'Time independant computation.
To get the sources estimations for a time frame,
 the kernel is applied to the recordings (matrix product).', @(h,ev)UpdatePanel(), []);
        jRadioFull   = gui_component('radio', jPanelOutput, 'br', 'Full results (Kernel*Recordings)', jGroupOutput, 'Compute sources for all the time samples.', @(h,ev)UpdatePanel(), []);
        if isShared
            jRadioFull.setEnabled(0);
            jRadioKernel.setSelected(1);
        elseif isFull
            jRadioFull.setSelected(1);
        else
            jRadioKernel.setSelected(1);
        end
    c.gridy = 4;
    jPanelRight.add(jPanelOutput, c);

    % ===== VALIDATION BUTTONS =====
    jPanelValid = gui_river([1,1], [0,6,6,6]);
    % Add a glue at the bottom of the right panel (for appropriate scaling with the left)
    c.gridy   = 5;
    c.weighty = 1;
    jPanelRight.add(Box.createVerticalGlue(), c);
    % Add a glue at the bottom of the left panel (for appropriate scaling with the right)
    c.gridy   = 5;
    c.weighty = 1;
    jPanelLeft.add(Box.createVerticalGlue(), c);
    % Expert/normal mode
    jButtonExpert = gui_component('button', jPanelValid, [], 'Show details', [], [], @SwitchExpertMode_Callback, []);
    gui_component('label', jPanelValid, 'hfill', ' ');
    % Ok/Cancel
    gui_component('Button', jPanelValid, 'right', 'Cancel', [], [], @ButtonCancel_Callback, []);
    gui_component('Button', jPanelValid, [],      'OK',     [], [], @ButtonOk_Callback,     []);
    c.gridy = 6;
    c.weighty = 0;
    jPanelLeft.add(jPanelValid, c);

    % ===== PANEL CREATION =====
    % Return a mutex to wait for panel close
    bst_mutex('create', panelName);
    % Create the BstPanel object that is returned by the function
    ctrl = struct(...
            'HeadModelType',  HeadModelType, ...
            'jTextComment',   jTextComment, ...
            ... % ==== PANEL: METHOD ====
            'jRadioMethodMn',  jRadioMethodMn, ...
            'jRadioMethodBf',  jRadioMethodBf, ...
            'jRadioMethodDip', jRadioMethodDip, ...
            ... % ==== PANEL: MEASURE ====
            'jRadioMnCurrent',   jRadioMnCurrent, ...
            'jRadioMnDspm',      jRadioMnDspm, ...
            'jRadioMnSloreta',   jRadioMnSloreta, ...
            'jRadioMethodBfNai', jRadioMethodBfNai, ...
            ... % ==== PANEL: SOURCE MODEL ====
            'jRadioConstr',   jRadioConstr, ...
            'jRadioUnconstr', jRadioUnconstr, ...
            'jRadioLoose',    jRadioLoose, ...
            'jTextLoose',     jTextLoose, ...
            ... % ==== PANEL: DEPTH WEIGHTING ====
            'jCheckDepth',      jCheckDepth, ...
            'jTextWeightExp',   jTextWeightExp, ...
            'jTextWeightLimit', jTextWeightLimit, ...
            ... % ==== PANEL: SNR ====
            'jRadioSnrRms',  jRadioSnrRms, ...
            'jTextSnrRms',   jTextSnrRms, ...
            'jRadioSnrFix',  jRadioSnrFix, ...
            'jTextSnrFix',   jTextSnrFix, ...
            ... % ==== PANEL: NOISE COVARIANCE ====
            'jRadioShrink',  jRadioShrink, ...
            'jRadioMedian',  jRadioMedian, ...
            'jRadioReg',     jRadioReg, ...
            'jTextReg',      jTextReg, ...
            'jRadioDiag',    jRadioDiag, ...
            'jRadioNoReg',   jRadioNoReg, ...
            ... % ==== PANEL: NON-LINEAR ====
            'jRadioMethodMem',     jRadioMethodMem, ...
            ... % ==== PANEL: DATA TYPE ====
            'jCheckMeg',     jCheckMeg, ...
            'jCheckMegGrad', jCheckMegGrad, ...
            'jCheckMegMag',  jCheckMegMag, ...
            'jCheckEeg',     jCheckEeg, ...
            'jCheckEcog',    jCheckEcog, ...
            'jCheckSeeg',    jCheckSeeg, ...
            ... % ===== PANEL: OUTPUT MODE =====
            'jRadioFull',    jRadioFull, ...
            'jRadioKernel',  jRadioKernel);
    % Create the BstPanel object that is returned by the function
    bstPanelNew = BstPanel(panelName, jPanelNew, ctrl);
    % Update comments
    UpdatePanel(1, 1);

%% =================================================================================
%  === LOCAL CALLBACKS  ============================================================
%  =================================================================================
    %% ===== BUTTON: CANCEL =====
    function ButtonCancel_Callback(varargin)
        % Close panel
        gui_hide(panelName);
    end

    %% ===== BUTTON: OK =====
    function ButtonOk_Callback(varargin)
        % Release mutex and keep the panel opened
        bst_mutex('release', panelName);
    end

    %% ===== METHOD CALLBACK =====
    function Method_Callback(hObject, event)
        % Change the default for source orientation
        if jRadioMethodMn.isSelected() && jRadioConstr.isEnabled()
            jRadioConstr.setSelected(1);
        elseif jRadioMethodDip.isSelected()
            jRadioUnconstr.setSelected(1);
        end
        % Update the panel
        UpdatePanel(1);
    end

    %% ===== MODALITY CALLBACK =====
    function Modality_Callback(hObject, event)
        % If only one checkbox: can't deselect it
        if (length(Modalities) == 1)
            event.getSource().setSelected(1);
        % Warning if both MEG and EEG are selected
        elseif isFirstCombinationWarning && ~isempty(jCheckEeg) && jCheckEeg.isSelected() && (...
                (~isempty(jCheckMeg) && jCheckMeg.isSelected()) || ...
                (~isempty(jCheckMegGrad) && jCheckMegGrad.isSelected()) || ...
                (~isempty(jCheckMegMag) && jCheckMegMag.isSelected()))
            java_dialog('warning', ['Warning: Brainstorm inverse models do not properly handle the combination of MEG and EEG yet.' 10 10 ...
                                    'For now, we recommend to compute separatly the sources for MEG and EEG.'], 'EEG/MEG combination');
            isFirstCombinationWarning = 0;
        end
        % Update comment
        UpdatePanel();
    end

    %% ===== SWITCH EXPERT MODE =====
    function SwitchExpertMode_Callback(varargin)
        % Toggle expert mode
        ExpertMode = bst_get('ExpertMode');
        bst_set('ExpertMode', ~ExpertMode);
        % Update comment
        UpdatePanel(1);
    end

    %% ===== UPDATE PANEL ======
    % USAGE:  UpdatePanel(isForced = 0)
    function UpdatePanel(isForced, isFirstCall)
        % Default values
        if (nargin < 2) || isempty(isFirstCall)
            isFirstCall = 0;
        end
        if (nargin < 1) || isempty(isForced)
            isForced = 0;
        end
        % Get the main categories of options
        isLinear = isempty(jRadioMethodMem) || ~jRadioMethodMem.isSelected();
        % Expert mode / Normal mode
        if isForced
            ExpertMode = bst_get('ExpertMode');
            % Left panels
            jPanelModel.setVisible(isLinear);
            jPanelMeasureMN.setVisible(isLinear && jRadioMethodMn.isSelected());
            jPanelMeasureBf.setVisible(isLinear && jRadioMethodBf.isSelected());
            jPanelMemInfo.setVisible(~isLinear);
            % Right panels (expert)
            jPanelRight.setVisible(ExpertMode);
            jPanelNoiseCov.setVisible(isLinear);
            jPanelSnr.setVisible(isLinear && ~jRadioMethodDip.isSelected() && ~jRadioMethodBf.isSelected());
            jPanelDepth.setVisible(isLinear && jRadioMethodMn.isSelected() && ~jRadioMnSloreta.isSelected());
            jPanelOutput.setVisible(isLinear && ~isProcess);
            % Update expert button
            if ExpertMode
                jButtonExpert.setText('Hide details');
            else
                jButtonExpert.setText('Show details');
            end
            % Matrix to regularize
            regMethod = '';
            if jRadioMethodBf.isSelected()
                jPanelNoiseCov.getBorder().setTitle('Data covariance regularization');
                jRadioReg.setText('Regularize data covariance:');
                jRadioDiag.setText(['Diagonal data covariance' strDiagRecommData]);
                if ~isempty(DataMethod)
                    regMethod = DataMethod;
                else
                    regMethod = 'median';
                end
            else
                jPanelNoiseCov.getBorder().setTitle('Noise covariance regularization');
                jRadioReg.setText('Regularize noise covariance:');
                jRadioDiag.setText(['Diagonal noise covariance' strDiagRecommNoise]);
                if jRadioMethodDip.isSelected() && ~isequal(NoiseMethod, 'diag')
                    regMethod = 'median';
                elseif ~isempty(NoiseMethod)
                    regMethod = NoiseMethod;
                end
            end
            jPanelNoiseCov.repaint();
            % Force selection for the covariance reg
            switch lower(regMethod)
                case 'reg',    jRadioReg.setSelected(1);
                case 'diag',   jRadioDiag.setSelected(1);
                case 'none',   jRadioNoReg.setSelected(1);
                case 'shrink', jRadioShrink.setSelected(1);
                case 'median', jRadioMedian.setSelected(1);
            end
            % Select default regularization method
            if jRadioMethodMn.isSelected()
                jRadioSnrFix.setSelected(1);
                jRadioSnrRms.setEnabled(0);
                if strcmpi(HeadModelType, 'surface')
                    jRadioLoose.setEnabled(1);
                end
            else
                jRadioSnrRms.setSelected(1);
                jRadioSnrRms.setEnabled(1);
                jRadioLoose.setEnabled(0);
                if jRadioLoose.isSelected()
                    jRadioUnconstr.setSelected(1);
                end
            end

            % Get old panel
            [bstPanelOld, iPanel] = bst_get('Panel', 'InverseOptions');
            container = get(bstPanelOld, 'container');
            % Re-pack frame
            if ~isempty(container)
                jFrame = container.handle{1};
                if ~isempty(jFrame)
                    jFrame.pack();
                end
            end
        end
        % Enable/disable text boxes
        jTextLoose.setEnabled(jRadioLoose.isSelected());
        jLabelWeightExp.setEnabled(jCheckDepth.isSelected());
        jTextWeightExp.setEnabled(jCheckDepth.isSelected());
        jLabelWeightLimit.setEnabled(jCheckDepth.isSelected());
        jTextWeightLimit.setEnabled(jCheckDepth.isSelected());
        jTextSnrRms.setEnabled(jRadioSnrRms.isSelected());
        jTextSnrFix.setEnabled(jRadioSnrFix.isSelected());
        jTextReg.setEnabled(jRadioReg.isSelected());

        % Selected modalities
        selModalities = {};
        if ~isempty(jCheckMeg) && jCheckMeg.isSelected()
            selModalities{end+1} = 'MEG';
        end
        if ~isempty(jCheckMegGrad) && jCheckMegGrad.isSelected()
            selModalities{end+1} = 'MEG GRAD';
        end
        if ~isempty(jCheckMegMag) && jCheckMegMag.isSelected()
            selModalities{end+1} = 'MEG MAG';
        end
        if ~isempty(jCheckEeg) && jCheckEeg.isSelected()
            selModalities{end+1} = 'EEG';
        end
        if ~isempty(jCheckEcog) && jCheckEcog.isSelected()
            selModalities{end+1} = 'ECOG';
        end
        if ~isempty(jCheckSeeg) && jCheckSeeg.isSelected()
            selModalities{end+1} = 'SEEG';
        end

        % Linear methods
        if isLinear
            % Get selected method
            [InverseMethod, InverseMeasure] = GetSelectedMethod(ctrl);
            % Get comment for this method
            Comment = GetMethodComment(InverseMethod, InverseMeasure);
        else
            if ~isempty(jRadioMethodMem) && jRadioMethodMem.isSelected()
                Comment = 'MEM: ';
            end
        end
        % Add modality comment
        Comment = [Comment, ': ', process_inverse_2018('GetModalityComment', selModalities)];
        % Update comment field
        if ~isProcess || ~isFirstCall || ~isfield(OPTIONS, 'Comment') || isempty(OPTIONS.Comment)
            jTextComment.setText(Comment);
        end
    end
end

%% =================================================================================
%  === EXTERNAL CALLBACKS  =========================================================
%  =================================================================================
%% ===== GET PANEL CONTENTS =====
function s = GetPanelContents() %#ok<DEFNU>
    % Get panel controls handles
    ctrl = bst_get('PanelControls', 'InverseOptions');
    if isempty(ctrl)
        s = [];
        return;
    end
    % Comment
    s.Comment = char(ctrl.jTextComment.getText());
    % Linear models
    isLinear = isempty(ctrl.jRadioMethodMem) || ~ctrl.jRadioMethodMem.isSelected();
    if isLinear
        % Get selected method
        [s.InverseMethod, s.InverseMeasure] = GetSelectedMethod(ctrl);
        % Source model
        if strcmpi(ctrl.HeadModelType, 'mixed')
            s.SourceOrient = [];
            s.Loose = [];
        else
            if ctrl.jRadioConstr.isSelected()
                s.SourceOrient = {'fixed'};
            elseif ctrl.jRadioUnconstr.isSelected()
                s.SourceOrient = {'free'};
            elseif ctrl.jRadioLoose.isSelected()
                s.SourceOrient = {'loose'};
            end
            s.Loose = str2num(char(ctrl.jTextLoose.getText()));
        end
        % Depth weighting
        s.UseDepth    = ctrl.jCheckDepth.isSelected() && ~strcmpi(s.InverseMeasure, 'sloreta');
        s.WeightExp   = str2num(char(ctrl.jTextWeightExp.getText()));
        s.WeightLimit = str2num(char(ctrl.jTextWeightLimit.getText()));
        % Noise covariance
        if ctrl.jRadioShrink.isSelected()
            s.NoiseMethod = 'shrink';
        elseif ctrl.jRadioMedian.isSelected()
            s.NoiseMethod = 'median';
        elseif ctrl.jRadioReg.isSelected()
            s.NoiseMethod = 'reg';
        elseif ctrl.jRadioDiag.isSelected()
            s.NoiseMethod = 'diag';
        elseif ctrl.jRadioNoReg.isSelected()
            s.NoiseMethod = 'none';
        end
        s.NoiseReg = str2num(char(ctrl.jTextReg.getText()));
        % Signal to noise
        if ctrl.jRadioSnrRms.isSelected()
            s.SnrMethod = 'rms';
        elseif ctrl.jRadioSnrFix.isSelected()
            s.SnrMethod = 'fixed';
        end
        s.SnrRms   = str2num(char(ctrl.jTextSnrRms.getText())) * 1e-9;  % Convert to Amper.mter
        s.SnrFixed = str2num(char(ctrl.jTextSnrFix.getText()));
        % Output mode
        if ctrl.jRadioFull.isSelected()
            s.ComputeKernel = 0;
        else
            s.ComputeKernel = 1;
        end
    % Non-linear models
    else
        % Get selected method
        if ctrl.jRadioMethodMem.isSelected()
            s.InverseMethod = 'mem';
        end
        % Other fields that are not defined
        s.InverseMeasure = [];
        s.SourceOrient   = [];
        s.Loose          = [];
        s.UseDepth       = [];
        s.WeightExp      = [];
        s.WeightLimit    = [];
        s.NoiseMethod    = [];
        s.NoiseReg       = [];
        s.SnrMethod      = [];
        s.SnrRms         = [];
        s.SnrFixed       = [];
        s.ComputeKernel  = 0;
    end
    % Selected modalities
    s.DataTypes = {};
    if ~isempty(ctrl.jCheckMeg) && ctrl.jCheckMeg.isSelected()
        s.DataTypes{end+1} = 'MEG';
    end
    if ~isempty(ctrl.jCheckMegGrad) && ctrl.jCheckMegGrad.isSelected()
        s.DataTypes{end+1} = 'MEG GRAD';
    end
    if ~isempty(ctrl.jCheckMegMag) && ctrl.jCheckMegMag.isSelected()
        s.DataTypes{end+1} = 'MEG MAG';
    end
    if ~isempty(ctrl.jCheckEeg) && ctrl.jCheckEeg.isSelected()
        s.DataTypes{end+1} = 'EEG';
    end
    if ~isempty(ctrl.jCheckEcog) && ctrl.jCheckEcog.isSelected()
        s.DataTypes{end+1} = 'ECOG';
    end
    if ~isempty(ctrl.jCheckSeeg) && ctrl.jCheckSeeg.isSelected()
        s.DataTypes{end+1} = 'SEEG';
    end
end

%% ===== GET SELECTED METHOD =====
function [Method, Measure] = GetSelectedMethod(ctrl)
    if ctrl.jRadioMethodMn.isSelected()
        Method = 'minnorm';
        if ctrl.jRadioMnCurrent.isSelected()
            Measure = 'amplitude';
        elseif ctrl.jRadioMnDspm.isSelected()
            Measure = 'dspm2018';
        elseif ctrl.jRadioMnSloreta.isSelected()
            Measure = 'sloreta';
        end
    elseif ctrl.jRadioMethodDip.isSelected()
        Method = 'gls';
        Measure = 'performance';
    elseif ctrl.jRadioMethodBf.isSelected()
        Method = 'lcmv';
        if ctrl.jRadioMethodBfNai.isSelected()
            Measure = 'nai';
        end
    end
end

%% ===== GET METHOD COMMENT =====
function Comment = GetMethodComment(Method, Measure)
    Comment = 'UKNOWN';
    switch (lower(Method))
        case 'minnorm'
            switch (lower(Measure))
                case 'amplitude', Comment = 'MN';
                case 'dspm2018',  Comment = 'dSPM-unscaled';
                case 'sloreta',   Comment = 'sLORETA';
            end
        case 'gls'
            Comment = 'Dipoles';
        case 'lcmv'
            switch (lower(Measure))
                case 'nai',       Comment = 'PNAI';
            end
        case 'mem'
            Comment = 'MEM';
    end
end

%% ===== WARNING DSPM 2018 =====
function WarningDspm()
    java_dialog('msgbox', [...
        'The dSPM implementation in "Compute sources [2018]" changed in July 2018:' 10 ...
        'The values are not scaled by the effective number of trials any more.' 10 ...
        'To get proper dSPM values for averages, run process "Sources > Scale averaged dSPM".' 10 10 ...
        'You will be directed to the Brainstorm website for additional information.' 10 10], 'Warning: dSPM update.');
    web('https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Averaging_normalized_values', '-browser');
end

author: brainstorm-tools


text:
function [TR, TT, ER, maxD, t] = icp(q,p,varargin)
% this is modified version of the original version
%
% Perform the Iterative Closest Point algorithm on three dimensional point
% clouds.
%
% [TR, TT] = icp(q,p)   returns the rotation matrix TR and translation
% vector TT that minimizes the distances from (TR * p + TT) to q.
% p is a 3xm matrix and q is a 3xn matrix.
%
% [TR, TT] = icp(q,p,k)   forces the algorithm to make k iterations
% exactly. The default is 10 iterations.
%
% [TR, TT, ER] = icp(q,p,k)   also returns the RMS of errors for k
% iterations in a (k+1)x1 vector. ER(0) is the initial error.
%
% [TR, TT, ER, t] = icp(q,p,k)   also returns the calculation times per
% iteration in a (k+1)x1 vector. t(0) is the time consumed for preprocessing.
%
% Additional settings may be provided in a parameter list:
%
% Boundary
%       {[]} | 1x? vector
%       If EdgeRejection is set, a vector can be provided that indexes into
%       q and specifies which points of q are on the boundary.
%
% EdgeRejection
%       {false} | true
%       If EdgeRejection is true, point matches to edge vertices of q are
%       ignored. Requires that boundary points of q are specified using
%       Boundary or that a triangulation matrix for q is provided.
%
% Extrapolation
%       {false} | true
%       If Extrapolation is true, the iteration direction will be evaluated
%       and extrapolated if possible using the method outlined by
%       Besl and McKay 1992.
%
% Matching
%       bruteForce | Delaunay | {kDtree}
%       Specifies how point matching should be done.
%       bruteForce is usually the slowest and kDtree is the fastest.
%       Note that the kDtree option is depends on the Statistics Toolbox
%       v. 7.3 or higher.
%
% Minimize
%       {point} | plane | lmaPoint
%       Defines whether point to point or point to plane minimization
%       should be performed. point is based on the SVD approach and is
%       usually the fastest. plane will often yield higher accuracy. It
%       uses linearized angles and requires surface normals for all points
%       in q. Calculation of surface normals requires substantial pre
%       proccessing.
%       The option lmaPoint does point to point minimization using the non
%       linear least squares Levenberg Marquardt algorithm. Results are
%       generally the same as in points, but computation time may differ.
%
% Normals
%       {[]} | n x 3 matrix
%       A matrix of normals for the n points in q might be provided.
%       Normals of q are used for point to plane minimization.
%       Else normals will be found through a PCA of the 4 nearest
%       neighbors.
%
% ReturnAll
%       {false} | true
%       Determines whether R and T should be returned for all iterations
%       or only for the last one. If this option is set to true, R will be
%       a 3x3x(k+1) matrix and T will be a 3x1x(k+1) matrix.
%
% Triangulation
%       {[]} | ? x 3 matrix
%       A triangulation matrix for the points in q can be provided,
%       enabling EdgeRejection. The elements should index into q, defining
%       point triples that act together as triangles.
%
% Verbose
%       {false} | true
%       Enables extrapolation output in the Command Window.
%
% Weight
%       {@(match)ones(1,m)} | Function handle
%       For point or plane minimization, a function handle to a weighting
%       function can be provided. The weighting function will be called
%       with one argument, a 1xm vector that specifies point pairs by
%       indexing into q. The weighting function should return a 1xm vector
%       of weights for every point pair.
%
% WorstRejection
%       {0} | scalar in ]0; 1[
%       Reject a given percentage of the worst point pairs, based on their
%       Euclidean distance.
%
% Martin Kjer and Jakob Wilm, Technical University of Denmark, 2012

% Use the inputParser class to validate input arguments.
inp = inputParser;

inp.addRequired('q', @(x)isreal(x) && size(x,1) == 3);
inp.addRequired('p', @(x)isreal(x) && size(x,1) == 3);

inp.addOptional('iter', 10, @(x)x > 0 && x < 10^5);

inp.addParamValue('Boundary', [], @(x)size(x,1) == 1);

inp.addParamValue('EdgeRejection', false, @(x)islogical(x));

inp.addParamValue('Extrapolation', false, @(x)islogical(x));

validMatching = {'bruteForce','Delaunay','kDtree'};
inp.addParamValue('Matching', 'kDtree', @(x)any(strcmpi(x,validMatching)));

validMinimize = {'point','plane','lmapoint'};
inp.addParamValue('Minimize', 'point', @(x)any(strcmpi(x,validMinimize)));

inp.addParamValue('Normals', [], @(x)isreal(x) && size(x,1) == 3);

inp.addParamValue('NormalsData', [], @(x)isreal(x) && size(x,1) == 3);

inp.addParamValue('ReturnAll', false, @(x)islogical(x));

inp.addParamValue('Triangulation', [], @(x)isreal(x) && size(x,2) == 3);

inp.addParamValue('Verbose', false, @(x)islogical(x));

inp.addParamValue('Weight', @(x)ones(1,length(x)), @(x)isa(x,'function_handle'));

inp.addParamValue('WorstRejection', 0, @(x)isscalar(x) && x > 0 && x < 1);

inp.addParamValue('SmartRejection', 0, @(x)isscalar(x) && x > 0);

inp.parse(q,p,varargin{:});
arg = inp.Results;
clear('inp');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Actual implementation

% Allocate vector for RMS of errors in every iteration.
t = zeros(arg.iter+1,1);

% Start timer
tic;

Np = size(p,2);

% Transformed data point cloud
pt = p;

% Allocate vector for RMS of errors in every iteration.
ER = zeros(arg.iter+1,1);
maxD = zeros(arg.iter,1);

% Initialize temporary transform vector and matrix.
T = zeros(3,1);
R = eye(3,3);

% Initialize total transform vector(s) and rotation matric(es).
TT = zeros(3,1, arg.iter+1);
TR = repmat(eye(3,3), [1,1, arg.iter+1]);

% If Minimize == 'plane', normals are needed
if (strcmp(arg.Minimize, 'plane') && isempty(arg.Normals))
    arg.Normals = lsqnormest(q,4);
end

% If Matching == 'Delaunay', a triangulation is needed
if strcmp(arg.Matching, 'Delaunay')
    DT = DelaunayTri(transpose(q));
end

% If Matching == 'kDtree', a kD tree should be built (req. Stat. TB >= 7.3)
if strcmp(arg.Matching, 'kDtree')
    kdOBJ = KDTreeSearcher(transpose(q));
end

% If edge vertices should be rejected, find edge vertices
if arg.EdgeRejection
    if isempty(arg.Boundary)
        bdr = find_bound(q, arg.Triangulation);
    else
        bdr = arg.Boundary;
    end
end

if arg.Extrapolation
    % Initialize total transform vector (quaternion ; translation vec.)
    qq = [ones(1,arg.iter+1);zeros(6,arg.iter+1)];
    % Allocate vector for direction change and change angle.
    dq = zeros(7,arg.iter+1);
    theta = zeros(1,arg.iter+1);
end

t(1) = toc;

% Go into main iteration loop
for k=1:arg.iter

    % Do matching
    switch arg.Matching
        case 'bruteForce'
            [match mindist] = match_bruteForce(q,pt);
        case 'Delaunay'
            [match mindist] = match_Delaunay(q,pt,DT);
        case 'kDtree'
            [match mindist] = match_kDtree(q,pt,kdOBJ);
    end

    % If matches to edge vertices should be rejected
    if arg.EdgeRejection
        p_idx = not(ismember(match, bdr));
        q_idx = match(p_idx);
        mindist = mindist(p_idx);
    else
        p_idx = true(1, Np);
        q_idx = match;
    end

    if k==1 && arg.SmartRejection
        arg.WorstRejection = sum(mindist > (median(mindist)* arg.SmartRejection))/ length(mindist);
        fprintf('ICP: Median=%f Threshold=%f WorstRejection=%f\n', median(mindist),median(mindist)* arg.SmartRejection,arg.WorstRejection);
        % vis
        %{
        figure
        hist(mindist,0:0.01/10:max(mindist)); hold on;
        plot(median(mindist),0,'g*');
        plot(median(mindist)* arg.SmartRejection,0,'r*');
        %}
    end

    % If worst matches should be rejected
    if arg.WorstRejection
        edge = round((1-arg.WorstRejection)*sum(p_idx));
        pairs = find(p_idx);
        [~, idx] = sort(mindist);
        p_idx(pairs(idx(edge:end))) = false;
        q_idx = match(p_idx);
        mindist = mindist(p_idx);
    end

    maxD(k) = max(mindist);

    if k == 1
        ER(k) = sqrt(sum(mindist.^2)/length(mindist));
    end

    switch arg.Minimize
        case 'point'
            % Determine weight vector
            weights = arg.Weight(match);
            [R,T] = eq_point(q(:,q_idx),pt(:,p_idx), weights(p_idx));
        case 'plane'
            weights = arg.Weight(match);
            [R,T] = eq_plane(q(:,q_idx),pt(:,p_idx),arg.Normals(:,q_idx),weights(p_idx));
        case 'lmaPoint'
            [R,T] = eq_lmaPoint(q(:,q_idx),pt(:,p_idx));
    end

    % Add to the total transformation
    TR(:,:,k+1) = R*TR(:,:,k);
    TT(:,:,k+1) = R*TT(:,:,k)+T;

    % Apply last transformation
    pt = TR(:,:,k+1) * p + repmat(TT(:,:,k+1), 1, Np);

    % Root mean of objective function
    ER(k+1) = rms_error(q(:,q_idx), pt(:,p_idx));

    % If Extrapolation, we might be able to move quicker
    if arg.Extrapolation
        qq(:,k+1) = [rmat2quat(TR(:,:,k+1));TT(:,:,k+1)];
        dq(:,k+1) = qq(:,k+1) - qq(:,k);
        theta(k+1) = (180/pi)*acos(dot(dq(:,k),dq(:,k+1))/(norm(dq(:,k))*norm(dq(:,k+1))));
        if arg.Verbose
            disp(['Direction change ' num2str(theta(k+1)) ' degree in iteration ' num2str(k)]);
        end
        if k>2 && theta(k+1) < 10 && theta(k) < 10
            d = [ER(k+1), ER(k), ER(k-1)];
            v = [0, -norm(dq(:,k+1)), -norm(dq(:,k))-norm(dq(:,k+1))];
            vmax = 25 * norm(dq(:,k+1));
            dv = extrapolate(v,d,vmax);
            if dv ~= 0
                q_mark = qq(:,k+1) + dv * dq(:,k+1)/norm(dq(:,k+1));
                q_mark(1:4) = q_mark(1:4)/norm(q_mark(1:4));
                qq(:,k+1) = q_mark;
                TR(:,:,k+1) = quat2rmat(qq(1:4,k+1));
                TT(:,:,k+1) = qq(5:7,k+1);
                % Reapply total transformation
                pt = TR(:,:,k+1) * p + repmat(TT(:,:,k+1), 1, Np);
                % Recalculate root mean of objective function
                % Note this is costly and only for fun!
                switch arg.Matching
                    case 'bruteForce'
                        [~, mindist] = match_bruteForce(q,pt);
                    case 'Delaunay'
                        [~, mindist] = match_Delaunay(q,pt,DT);
                    case 'kDtree'
                        [~, mindist] = match_kDtree(q,pt,kdOBJ);
                end
                ER(k+1) = sqrt(sum(mindist.^2)/length(mindist));
            end
        end
    end
    t(k+1) = toc;
end

if not(arg.ReturnAll)
    TR = TR(:,:,end);
    TT = TT(:,:,end);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [match mindist] = match_bruteForce(q, p)
    m = size(p,2);
    n = size(q,2);
    match = zeros(1,m);
    mindist = zeros(1,m);
    for ki=1:m
        d=zeros(1,n);
        for ti=1:3
            d=d+(q(ti,:)-p(ti,ki)).^2;
        end
        [mindist(ki),match(ki)]=min(d);
    end

    mindist = sqrt(mindist);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [match mindist] = match_Delaunay(q, p, DT)
	match = transpose(nearestNeighbor(DT, transpose(p)));
	mindist = sqrt(sum((p-q(:,match)).^2,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [match mindist] = match_kDtree(~, p, kdOBJ)
	[match mindist] = knnsearch(kdOBJ,transpose(p));
    match = transpose(match);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [R,T] = eq_point(q,p,weights)

m = size(p,2);
n = size(q,2);

% normalize weights
weights = weights ./ sum(weights);

% find data centroid and deviations from centroid
q_bar = q * transpose(weights);
q_mark = q - repmat(q_bar, 1, n);
% Apply weights
q_mark = q_mark .* repmat(weights, 3, 1);

% find data centroid and deviations from centroid
p_bar = p * transpose(weights);
p_mark = p - repmat(p_bar, 1, m);
% Apply weights
%p_mark = p_mark .* repmat(weights, 3, 1);

N = p_mark*transpose(q_mark); % taking points of q in matched order

[U,~,V] = svd(N); % singular value decomposition

R = V*diag([1 1 det(U*V')])*transpose(U);

T = q_bar - R*p_bar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [R,T] = eq_plane(q,p,n,weights)

n = n .* repmat(weights,3,1);

c = cross(p,n);

cn = vertcat(c,n);

C = cn*transpose(cn);

b = - [sum(sum((p-q).*repmat(cn(1,:),3,1).*n));
       sum(sum((p-q).*repmat(cn(2,:),3,1).*n));
       sum(sum((p-q).*repmat(cn(3,:),3,1).*n));
       sum(sum((p-q).*repmat(cn(4,:),3,1).*n));
       sum(sum((p-q).*repmat(cn(5,:),3,1).*n));
       sum(sum((p-q).*repmat(cn(6,:),3,1).*n))];

X = C\b;

cx = cos(X(1)); cy = cos(X(2)); cz = cos(X(3));
sx = sin(X(1)); sy = sin(X(2)); sz = sin(X(3));

R = [cy*cz cz*sx*sy-cx*sz cx*cz*sy+sx*sz;
     cy*sz cx*cz+sx*sy*sz cx*sy*sz-cz*sx;
     -sy cy*sx cx*cy];

T = X(4:6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [R,T] = eq_lmaPoint(q,p)

Rx = @(a)[1     0       0;
          0     cos(a)  -sin(a);
          0     sin(a)  cos(a)];

Ry = @(b)[cos(b)    0   sin(b);
          0         1   0;
          -sin(b)   0   cos(b)];

Rz = @(g)[cos(g)    -sin(g) 0;
          sin(g)    cos(g)  0;
          0         0       1];

Rot = @(x)Rx(x(1))*Ry(x(2))*Rz(x(3));

myfun = @(x,xdata)Rot(x(1:3))*xdata+repmat(x(4:6),1,length(xdata));

options = optimset('Algorithm', 'levenberg-marquardt');
x = lsqcurvefit(myfun, zeros(6,1), p, q, [], [], options);

R = Rot(x(1:3));
T = x(4:6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extrapolation in quaternion space. Details are found in:
%
% Besl, P., & McKay, N. (1992). A method for registration of 3-D shapes.
% IEEE Transactions on pattern analysis and machine intelligence, 239?256.

function [dv] = extrapolate(v,d,vmax)

p1 = polyfit(v,d,1); % linear fit
p2 = polyfit(v,d,2); % parabolic fit
v1 = -p1(2)/p1(1); % linear zero crossing
v2 = -p2(2)/(2*p2(1)); % polynomial top point

if issorted([0 v2 v1 vmax]) || issorted([0 v2 vmax v1])
    disp('Parabolic update!');
    dv = v2;
elseif issorted([0 v1 v2 vmax]) || issorted([0 v1 vmax v2])...
        || (v2 < 0 && issorted([0 v1 vmax]))
    disp('Line based update!');
    dv = v1;
elseif v1 > vmax && v2 > vmax
    disp('Maximum update!');
    dv = vmax;
else
    disp('No extrapolation!');
    dv = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Determine the RMS error between two point equally sized point clouds with
% point correspondance.
% ER = rms_error(p1,p2) where p1 and p2 are 3xn matrices.

function ER = rms_error(p1,p2)
dsq = sum(power(p1 - p2, 2),1);
ER = sqrt(mean(dsq));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Converts (orthogonal) rotation matrices R to (unit) quaternion
% representations
%
% Input: A 3x3xn matrix of rotation matrices
% Output: A 4xn matrix of n corresponding quaternions
%
% http://en.wikipedia.org/wiki/Rotation_matrix#Quaternion

function quaternion = rmat2quat(R)

Qxx = R(1,1,:);
Qxy = R(1,2,:);
Qxz = R(1,3,:);
Qyx = R(2,1,:);
Qyy = R(2,2,:);
Qyz = R(2,3,:);
Qzx = R(3,1,:);
Qzy = R(3,2,:);
Qzz = R(3,3,:);

w = 0.5 * sqrt(1+Qxx+Qyy+Qzz);
x = 0.5 * sign(Qzy-Qyz) .* sqrt(1+Qxx-Qyy-Qzz);
y = 0.5 * sign(Qxz-Qzx) .* sqrt(1-Qxx+Qyy-Qzz);
z = 0.5 * sign(Qyx-Qxy) .* sqrt(1-Qxx-Qyy+Qzz);

quaternion = reshape([w;x;y;z],4,[]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Converts (unit) quaternion representations to (orthogonal) rotation matrices R
%
% Input: A 4xn matrix of n quaternions
% Output: A 3x3xn matrix of corresponding rotation matrices
%
% http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#From_a_quaternion_to_an_orthogonal_matrix

function R = quat2rmat(quaternion)
q0(1,1,:) = quaternion(1,:);
qx(1,1,:) = quaternion(2,:);
qy(1,1,:) = quaternion(3,:);
qz(1,1,:) = quaternion(4,:);

R = [q0.^2+qx.^2-qy.^2-qz.^2 2*qx.*qy-2*q0.*qz 2*qx.*qz+2*q0.*qy;
     2*qx.*qy+2*q0.*qz q0.^2-qx.^2+qy.^2-qz.^2 2*qy.*qz-2*q0.*qx;
     2*qx.*qz-2*q0.*qy 2*qy.*qz+2*q0.*qx q0.^2-qx.^2-qy.^2+qz.^2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Least squares normal estimation from point clouds using PCA
%
% H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle.
% Surface reconstruction from unorganized points.
% In Proceedings of ACM Siggraph, pages 71:78, 1992.
%
% p should be a matrix containing the horizontally concatenated column
% vectors with points. k is a scalar indicating how many neighbors the
% normal estimation is based upon.
%
% Note that for large point sets, the function performs significantly
% faster if Statistics Toolbox >= v. 7.3 is installed.
%
% Jakob Wilm 2010

function n = lsqnormest(p, k)
m = size(p,2);
n = zeros(3,m);

v = ver('stats');
if str2double(v.Version) >= 7.5
    neighbors = transpose(knnsearch(transpose(p), transpose(p), 'k', k+1));
else
    neighbors = k_nearest_neighbors(p, p, k+1);
end

for i = 1:m
    x = p(:,neighbors(2:end, i));
    p_bar = 1/k * sum(x,2);

    P = (x - repmat(p_bar,1,k)) * transpose(x - repmat(p_bar,1,k)); %spd matrix P
    %P = 2*cov(x);

    [V,D] = eig(P);

    [~, idx] = min(diag(D)); % choses the smallest eigenvalue

    n(:,i) = V(:,idx);   % returns the corresponding eigenvector
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Program to find the k - nearest neighbors (kNN) within a set of points.
% Distance metric used: Euclidean distance
%
% Note that this function makes repetitive use of min(), which seems to be
% more efficient than sort() for k < 30.

function [neighborIds neighborDistances] = k_nearest_neighbors(dataMatrix, queryMatrix, k)

numDataPoints = size(dataMatrix,2);
numQueryPoints = size(queryMatrix,2);

neighborIds = zeros(k,numQueryPoints);
neighborDistances = zeros(k,numQueryPoints);

D = size(dataMatrix, 1); %dimensionality of points

for i=1:numQueryPoints
    d=zeros(1,numDataPoints);
    for t=1:D % this is to avoid slow repmat()
        d=d+(dataMatrix(t,:)-queryMatrix(t,i)).^2;
    end
    for j=1:k
        [s,t] = min(d);
        neighborIds(j,i)=t;
        neighborDistances(j,i)=sqrt(s);
        d(t) = NaN; % remove found number from d
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Boundary point determination. Given a set of 3D points and a
% corresponding triangle representation, returns those point indices that
% define the border/edge of the surface.

function bound = find_bound(pts, poly)

%Correcting polygon indices and converting datatype
poly = double(poly);
pts = double(pts);

%Calculating freeboundary points:
TR = TriRep(poly, pts(1,:)', pts(2,:)', pts(3,:)');
FF = freeBoundary(TR);

%Output
bound = FF(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: jianxiongxiao

repo: ProfXkit


text:
function tests = test_metaImageRead
% Test suite for the file metaImageRead.
%
%   Test suite for the file metaImageRead
%
%   Example
%   test_metaImageRead
%
%   See also
%     metaImageRead

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2021-05-14,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2021 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_UInt8_2D(testCase)

fileName = 'img_UInt8_rampXY_20x10.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 2);
assertEqual(testCase, size(res), [10 20]);

function test_UInt8_2D_NoExtension(testCase)

fileName = 'img_UInt8_rampXY_20x10';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 2);
assertEqual(testCase, size(res), [10 20]);

function test_RGB8_2D_peppers(testCase)

fileName = 'img_RGB8_peppers.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [384 512 3]);

function test_UInt8_3D_RampXYZ(testCase)

fileName = 'img_UInt8_rampXYZ_6x5x4.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1,1,1), uint8(0));
assertEqual(testCase, res(1,end,1), uint8(5));
assertEqual(testCase, res(end,1,1), uint8(40));
assertEqual(testCase, res(end,end,1), uint8(45));
assertEqual(testCase, res(end,end,end), uint8(255));

function test_read_slice_list(testCase)

filename = 'BRNOR39e5p1List.mhd';
img = metaImageRead(fullfile('ratBrainMriSlices', filename));

exp = [96 96 96];
assertEqual(testCase, size(img), exp, 'image does not have the right size');

assertTrue(max(img(:))>0);

function test_read_slices_pattern(testCase)

filename = 'BRNOR39e5p1Pattern.mhd';
img = metaImageRead(fullfile('ratBrainMriSlices', filename));

exp = [96 96 96];
assertEqual(testCase, size(img), exp, 'image does not have the right size');

assertTrue(max(img(:))>0);

function test_UInt16_3D_RampXYZ_msb(testCase)

fileName = 'img_UInt16_rampXYZ_6x5x4_msb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), uint16(0));
assertEqual(testCase, res(1, 6, 1), uint16(5));
assertEqual(testCase, res(5, 1, 1), uint16(40));
assertEqual(testCase, res(5, 6, 1), uint16(45));
assertEqual(testCase, res(1, 1, 4), uint16(300));
assertEqual(testCase, res(1, 6, 4), uint16(305));
assertEqual(testCase, res(5, 1, 4), uint16(340));
assertEqual(testCase, res(5, 6, 4), uint16(345));

function test_UInt16_3D_RampXYZ_lsb(testCase)

fileName = 'img_UInt16_rampXYZ_6x5x4_lsb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), uint16(0));
assertEqual(testCase, res(1, 6, 1), uint16(5));
assertEqual(testCase, res(5, 1, 1), uint16(40));
assertEqual(testCase, res(5, 6, 1), uint16(45));
assertEqual(testCase, res(1, 1, 4), uint16(300));
assertEqual(testCase, res(1, 6, 4), uint16(305));
assertEqual(testCase, res(5, 1, 4), uint16(340));
assertEqual(testCase, res(5, 6, 4), uint16(345));

function test_Int16_3D_RampXYZ_msb(testCase)

fileName = 'img_Int16_rampXYZ_6x5x4_msb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), int16(0));
assertEqual(testCase, res(1, 6, 1), int16(5));
assertEqual(testCase, res(5, 1, 1), int16(40));
assertEqual(testCase, res(5, 6, 1), int16(45));
assertEqual(testCase, res(1, 1, 4), int16(300));
assertEqual(testCase, res(1, 6, 4), int16(305));
assertEqual(testCase, res(5, 1, 4), int16(340));
assertEqual(testCase, res(5, 6, 4), int16(345));

function test_Int16_3D_RampXYZ_lsb(testCase)

fileName = 'img_Int16_rampXYZ_6x5x4_lsb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), int16(0));
assertEqual(testCase, res(1, 6, 1), int16(5));
assertEqual(testCase, res(5, 1, 1), int16(40));
assertEqual(testCase, res(5, 6, 1), int16(45));
assertEqual(testCase, res(1, 1, 4), int16(300));
assertEqual(testCase, res(1, 6, 4), int16(305));
assertEqual(testCase, res(5, 1, 4), int16(340));
assertEqual(testCase, res(5, 6, 4), int16(345));

function test_Int32_3D_RampXYZ_msb(testCase)

fileName = 'img_Int32_rampXYZ_6x5x4_msb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), int32(0));
assertEqual(testCase, res(1, 6, 1), int32(5));
assertEqual(testCase, res(5, 1, 1), int32(40));
assertEqual(testCase, res(5, 6, 1), int32(45));
assertEqual(testCase, res(1, 1, 4), int32(300));
assertEqual(testCase, res(1, 6, 4), int32(305));
assertEqual(testCase, res(5, 1, 4), int32(340));
assertEqual(testCase, res(5, 6, 4), int32(345));

function test_Int32_3D_RampXYZ_lsb(testCase)

fileName = 'img_Int32_rampXYZ_6x5x4_lsb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), int32(0));
assertEqual(testCase, res(1, 6, 1), int32(5));
assertEqual(testCase, res(5, 1, 1), int32(40));
assertEqual(testCase, res(5, 6, 1), int32(45));
assertEqual(testCase, res(1, 1, 4), int32(300));
assertEqual(testCase, res(1, 6, 4), int32(305));
assertEqual(testCase, res(5, 1, 4), int32(340));
assertEqual(testCase, res(5, 6, 4), int32(345));

function test_Float32_3D_RampXYZ_msb(testCase)

fileName = 'img_Float32_rampXYZ_6x5x4_msb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), single(0));
assertEqual(testCase, res(1, 6, 1), single(5));
assertEqual(testCase, res(5, 1, 1), single(40));
assertEqual(testCase, res(5, 6, 1), single(45));
assertEqual(testCase, res(1, 1, 4), single(300));
assertEqual(testCase, res(1, 6, 4), single(305));
assertEqual(testCase, res(5, 1, 4), single(340));
assertEqual(testCase, res(5, 6, 4), single(345));

function test_Float32_3D_RampXYZ_lsb(testCase)

fileName = 'img_Float32_rampXYZ_6x5x4_lsb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), single(0));
assertEqual(testCase, res(1, 6, 1), single(5));
assertEqual(testCase, res(5, 1, 1), single(40));
assertEqual(testCase, res(5, 6, 1), single(45));
assertEqual(testCase, res(1, 1, 4), single(300));
assertEqual(testCase, res(1, 6, 4), single(305));
assertEqual(testCase, res(5, 1, 4), single(340));
assertEqual(testCase, res(5, 6, 4), single(345));

function test_Float64_3D_RampXYZ_msb(testCase)

fileName = 'img_Float64_rampXYZ_6x5x4_msb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), 0);
assertEqual(testCase, res(1, 6, 1), 5);
assertEqual(testCase, res(5, 1, 1), 40);
assertEqual(testCase, res(5, 6, 1), 45);
assertEqual(testCase, res(1, 1, 4), 300);
assertEqual(testCase, res(1, 6, 4), 305);
assertEqual(testCase, res(5, 1, 4), 340);
assertEqual(testCase, res(5, 6, 4), 345);

function test_Float64_3D_RampXYZ_lsb(testCase)

fileName = 'img_Float64_rampXYZ_6x5x4_lsb.mhd';

res = metaImageRead(fullfile('mhd_images', fileName));

assertEqual(testCase, ndims(res), 3);
assertEqual(testCase, size(res), [5 6 4]);

assertEqual(testCase, res(1, 1, 1), 0);
assertEqual(testCase, res(1, 6, 1), 5);
assertEqual(testCase, res(5, 1, 1), 40);
assertEqual(testCase, res(5, 6, 1), 45);
assertEqual(testCase, res(1, 1, 4), 300);
assertEqual(testCase, res(1, 6, 4), 305);
assertEqual(testCase, res(5, 1, 4), 340);
assertEqual(testCase, res(5, 6, 4), 345);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function[]=makefigs_mspec
%MAKEFIGS_MSPEC  Makes some sample figures for MSPEC.

%First figure
load bravo94
use bravo94
use bravo94.rcm

figure
cv=cv(:,5);

clear psi f spp snn spn
psi{1}=ones(size(cv))./sqrt(length(cv));
psi{2}=sleptap(length(cv),4);
psi{3}=sleptap(length(cv),32);

for i=1:length(psi)
    [f(:,i),spp(:,i),snn(:,i),spn(:,i)]=mspec(cv,conj(cv),psi{i});
end

h=twospecplot(f,spp,snn);
axes(h(1)),vlines(abs(corfreq(lat))), linestyle D b r k:
ax=axis;axis([10^-2.95 ax(2) 10^-3 10^5]),xtick(10.^[-3 -2 -1 0])
xlabel('Frequency (rad/hour)'),ylabel('Power Spectral Density')
axes(h(2)),vlines(abs(corfreq(lat))), linestyle D b r k:
ax=axis;axis([10^-2.95 ax(2) 10^-3 10^5]),xtick(10.^[-3 -2 -1 0])
xlabel('Frequency (rad/hour)')

%To print
if 0
    set(gcf,'paperposition',[1 1 10 5.5])
    currentdir=pwd;
    cd([whichdir('jlab_license') '/figures'])
    print -dpng mspec
    crop mspec.png
    cd(currentdir)
end

%Former example
load bravo94
x=bravo94.rcm.cv;
vswap(x,nan,0);
[psi,lambda]=sleptap(length(x),16);
[f,sp,sn,spn]=mspec(x,psi);
[f,su,sv,suv]=mspec(real(x),imag(x),psi);

figure,plot(f,[sp sn]),xlog,ylog,axis tight
title('Counterclockwise (blue) and clockwise (red) spectra'),
linestyle b b b b b b r r r r r r

load bravo94
x=bravo94.rcm.cv;
vswap(x,nan,0);
[psi,lambda]=sleptap(length(x),8);

for i=1:size(x,2)
    [f,Suu,Suu3,Cuu]=mspec(real(x(:,i)),real(x(:,3)),psi);
    gammauu(:,i)=Cuu./sqrt(Suu.*Suu3);
end

figure,
plot(f,abs(gammauu)),xlog,yoffset 1,axis tight
title('Coherence of u(t) at each depth vs. u(t) at \#3','interpreter','latex')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function menu_fieldtrip(handle, cfg, allowsubplot)

% MENU_FIELDTRIP adds a FieldTrip-specific menu to a figure.
%
% See also MENU_VIEWPOINT

if nargin<2
  cfg = [];
end

if nargin<3
  allowsubplot = true;
end

if ~allowsubplot && numel(findobj(gcf, 'type', 'axes')) > 1
  % do not add a menu if there are multiple axes, i.e. subplots
  % these may contain different data structures
  return
end

if ft_platform_supports('uimenu')
  % delete any possibly existing previous menu, this is safe because delete([]) does nothing
  delete(findobj(handle, 'type', 'uimenu', 'label', 'FieldTrip'));

  % add the menu with callbacks
  ftmenu = uimenu(handle, 'Label', 'FieldTrip');
  if ~isempty(cfg)
    uimenu(ftmenu, 'Label', 'Show pipeline', 'Callback', {@cb_pipeline, cfg});
  end
  uimenu(ftmenu, 'Label', 'About', 'Callback', @cb_about);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_pipeline(handle, eventdata, varargin)
ft_analysispipeline([], varargin{1});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_about(handle, eventdata, varargin)
msgbox( {
  'FieldTrip is the MATLAB toolbox for MEG and EEG analysis that is being developed at the Centre for Cognitive Neuroimaging of the Donders Institute for Brain, Cognition and Behaviour together with collaborating institutes. The FieldTrip software is released as open source under the GNU general public license.'
  ''
  sprintf('This is FieldTrip, version %s', ft_version)
  } );

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypefun {[@var{a}, @var{b}, @dots{}] =} pycall_sympy__ (@var{cmd}, @var{x}, @var{y}, @dots{})
%% Run some Python command on some objects and return other objects.
%%
%% This function is not really intended for end users.
%%
%% Here @var{cmd} is a string of Python code.
%% Inputs @var{x}, @var{y}, @dots{} can be a variety of objects
%% (possible types listed below). Outputs @var{a}, @var{b}, @dots{} are
%% converted from Python objects: not all types are possible, see
%% below.
%%
%% Example:
%% @example
%% @group
%% x = 10; y = 2;
%% cmd = '(x, y) = _ins; return (x+y, x-y)';
%% [a, b] = pycall_sympy__ (cmd, x, y)
%%   @result{} a =  12
%%   @result{} b =  8
%% @end group
%% @end example
%%
%% The inputs will be in a list called @code{_ins}.
%% Instead of @code{return}, you can append to the Python list
%% @code{_outs}:
%% @example
%% @group
%% cmd = '(x, y) = _ins; _outs.append(x**y)';
%% a = pycall_sympy__ (cmd, x, y)
%%   @result{} a =  100
%% @end group
%% @end example
%%
%% If you want to return a list, one way is to append a comma
%% to the return command.  Compare these two examples:
%% @example
%% @group
%% L = pycall_sympy__ ('return [1, -3.4, "python"],')
%%   @result{} L =
%%     @{
%%       [1,1] =  1
%%       [1,2] = -3.4000
%%       [1,3] = python
%%     @}
%% [a, b, c] = pycall_sympy__ ('return [1, -3.4, "python"]')
%%   @result{} a =  1
%%   @result{} b = -3.4000
%%   @result{} c = python
%% @end group
%% @end example
%%
%% You can also pass a cell-array of lines of code.  But be careful
%% with whitespace: its Python!
%% @example
%% @group
%% cmd = @{ '(x,) = _ins'
%%         'if x.is_Matrix:'
%%         '    return x.T'
%%         'else:'
%%         '    return x' @};
%% @end group
%% @end example
%% The cell array can be either a row or a column vector.
%% Each of these strings probably should not have any newlines
%% (other than escaped ones e.g., inside strings).  An exception
%% might be Python triple-quoted """ multiline strings """.
%% FIXME: test this.
%% It might be a good idea to avoid blank lines as they can cause
%% problems with some of the IPC mechanisms.
%%
%% Possible input types:
%% @itemize
%% @item sym objects
%% @item strings (char)
%% @item scalar doubles
%% @item structs
%% @end itemize
%% They can also be cell arrays of these items.  Multi-D cell
%% arrays may not work properly.
%%
%% Possible output types:
%% @itemize
%% @item SymPy objects
%% @item int
%% @item float
%% @item string
%% @item unicode strings
%% @item bool
%% @item dict (converted to structs)
%% @item lists/tuples (converted to cell vectors)
%% @end itemize
%% FIXME: add a py_config to change the header?  The python
%% environment is defined in python_header.py.  Changing it is
%% currently harder than it should be.
%%
%% Note: if you don't pass in any syms, this shouldn't need SymPy.
%% But it still imports it in that case.  If  you want to run this
%% w/o having the SymPy package, you'd need to hack a bit.
%%
%% @end deftypefun

function varargout = pycall_sympy__ (cmd, varargin)

  if (~iscell(cmd))
    if (isempty(cmd))
      cmd = {};
    else
      cmd = {cmd};
    end
  end

  % empty command will cause empty try: except: block
  if isempty(cmd)
    cmd = {'pass'};
  end

  %% IPC interface
  % the ipc mechanism shall put the input variables in the tuple
  % '_ins' and it will return to us whatever we put in the tuple
  % '_outs'.  There is no particular reason this needs to define
  % a function, I just thought it isolates local variables a bit.

  % Careful: fix this constant if you change the code below.
  % Test with "pycall_sympy__ ('raise')" which should say "line 1".
  LinesBeforeCmdBlock = 3;

  % replace blank lines w/ empty comments (unnec. b/c of try:?)
  %I = cellfun(@isempty, cmd);
  %cmd(I) = repmat({'#'}, 1, nnz(I));

  cmd = indent_lines(cmd, 8);
  cmd = { 'def _fcn(_ins):' ...
          '    _outs = []' ...
          '    try:' ...
          cmd{:} ...
          '    except Exception as e:' ...
          '        ers = type(e).__name__ + ": " + str(e) if str(e) else type(e).__name__' ...
          '        _outs = ("COMMAND_ERROR_PYTHON", ers, sys.exc_info()[-1].tb_lineno)' ...
          '    return _outs' ...
        };

  [A, db] = python_ipc_driver('run', cmd, varargin{:});

  if (~iscell(A))
    A={A};
  end

  %% Error reporting
  % ipc drivers are supposed to give back these specially formatting error strings
  if (~isempty(A) && ischar(A{1}) && strcmp(A{1}, 'COMMAND_ERROR_PYTHON'))
    errcmdlineno = A{3} - db.prelines;
    errlineno = errcmdlineno - LinesBeforeCmdBlock;
    if (errcmdlineno <= 0 || errcmdlineno > length (cmd))
      error ('Python exception: %s\n    occurred at unexpected line number, maybe near line %d?', ...
              A{2}, errlineno);
    else
      error ('Python exception: %s\n    occurred at line %d of the Python code block:\n    %s', ...
              A{2}, errlineno, strtrim (cmd{errcmdlineno}));
    end
  elseif (~isempty(A) && ischar(A{1}) && strcmp(A{1}, 'INTERNAL_PYTHON_ERROR'))
    % Here A{3} is the error msg and A{2} is more info about where it happened
    if (strfind (A{3}, 'KeyboardInterrupt'))
      what_to_do = ['Probably something was interrupted by "Ctrl-C".\n' ...
                    '    Do "sympref reset" and repeat your command.'];
    else
      what_to_do = ['Try "sympref reset" and repeat your command?\n' ...
                    '    (consider filing an issue at ' ...
                    'https://github.com/cbm755/octsympy/issues)'];
    end
    error ('Python exception: %s\n    occurred %s.\n    %s', ...
           A{3}, A{2}, sprintf (what_to_do))
  end

  M = length(A);
  varargout = cell(1,M);
  for i=1:M
    varargout{i} = A{i};
  end

  if nargout ~= M
    warning('number of outputs don''t match, was this intentional?')
  end
end

%!test
%! % general test
%! x = 10; y = 6;
%! cmd = '(x,y) = _ins; return (x+y,x-y)';
%! [a,b] = pycall_sympy__ (cmd, x, y);
%! assert (a == x + y && b == x - y)

%!test
%! % bool
%! assert (pycall_sympy__ ('return True,'))
%! assert (~pycall_sympy__ ('return False,'))

%!test
%! % float
%! assert (abs(pycall_sympy__ ('return 1.0/3,') - 1/3) < 1e-15)

%!test
%! % int
%! r = pycall_sympy__ ('return 123456');
%! assert (r == 123456)
%! assert (isinteger (r))

%!test
%! % string
%! x = 'octave';
%! cmd = 's = _ins[0]; return s.capitalize(),';
%! y = pycall_sympy__ (cmd, x);
%! assert (strcmp(y, 'Octave'))

%!test
%! % string with escaped newlines, comes back as escaped newlines
%! x = 'a string\nbroke off\nmy guitar\n';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % string with actual newlines, comes back as actual newlines
%! x = sprintf('a string\nbroke off\nmy guitar\n');
%! y = pycall_sympy__ ('return _ins', x);
%! y2 = strrep(y, sprintf('\n'), sprintf('\r\n'));  % windows
%! assert (strcmp(x, y) || strcmp(x, y2))

%!test
%! % cmd string with newlines, works with cell
%! y = pycall_sympy__ ('return "string\nbroke",');
%! y2 = sprintf('string\nbroke');
%! y3 = strrep(y2, sprintf('\n'), sprintf('\r\n'));  % windows
%! assert (strcmp(y, y2) || strcmp(y, y3))

%!test
%! % string with XML escapes
%! x = '<> >< <<>>';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))
%! x = '&';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % strings with double quotes
%! x = 'a\"b\"c';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))
%! x = '\"';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % cmd has double quotes, these must be escaped by user
%! % (of course: she is writing python code)
%! expy = 'a"b"c';
%! y = pycall_sympy__ ('return "a\"b\"c",');
%! assert (strcmp(y, expy))

%!test
%! % strings with quotes
%! x = 'a''b';  % this is a single quote
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % strings with quotes
%! x = '\"a''b\"c''\"d';
%! y = pycall_sympy__ ('return _ins[0]', x);
%! assert (strcmp(y, x))

%!test
%! % strings with quotes
%! expy = '"a''b"c''"d';
%! y = pycall_sympy__ ('s = "\"a''b\"c''\"d"; return s');
%! assert (strcmp(y, expy))

%!test
%! % strings with printf escapes
%! x = '% %% %%% %%%% %s %g %%s';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % cmd with printf escapes
%! x = '% %% %%% %%%% %s %g %%s';
%! y = pycall_sympy__ (['return "' x '",']);
%! assert (strcmp(y, x))

%!test
%! % cmd w/ backslash and \n must be escaped by user
%! expy = 'a\b\\c\nd\';
%! y = pycall_sympy__ ('return "a\\b\\\\c\\nd\\",');
%! assert (strcmp(y, expy))

%!test
%! % slashes
%! x = '/\\ // \\\\ \\/\\/\\';
%! z = '/\ // \\ \/\/\';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % slashes
%! z = '/\ // \\ \/\/\';
%! y = pycall_sympy__ ('return "/\\ // \\\\ \\/\\/\\"');
%! assert (strcmp(y, z))

%!test
%! % strings with special chars
%! x = '!@#$^&* you!';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))
%! x = '~-_=+[{]}|;:,.?';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!xtest
%! % string with backtick trouble for system -c (sysoneline)
%! x = '`';
%! y = pycall_sympy__ ('return _ins', x);
%! assert (strcmp(y, x))

%!test
%! % unicode
%! s1 = '我爱你';
%! cmd = 'return u"\u6211\u7231\u4f60",';
%! s2 = pycall_sympy__ (cmd);
%! assert (strcmp (s1, s2))

%!test
%! % unicode with \x
%! s1 = '我';
%! cmd = 'return b"\xe6\x88\x91".decode("utf-8")';
%! s2 = pycall_sympy__ (cmd);
%! assert (strcmp (s1, s2))

%!test
%! % unicode with \x and some escaped backslashes
%! s1 = '\我\';
%! cmd = 'return b"\\\xe6\x88\x91\\".decode("utf-8")';
%! s2 = pycall_sympy__ (cmd);
%! assert (strcmp (s1, s2))

%!xtest
%! % unicode passthru
%! s = '我爱你';
%! s2 = pycall_sympy__ ('return _ins', s);
%! assert (strcmp (s, s2))
%! s = '我爱你<>\&//\#%% %\我';
%! s2 = pycall_sympy__ ('return _ins', s);
%! assert (strcmp (s, s2))

%!xtest
%! % unicode w/ slashes, escapes
%! s = '我<>\&//\#%% %\我';
%! s2 = pycall_sympy__ ('return "我<>\\&//\\#%% %\\我"');
%! assert (strcmp (s, s2))

%!test
%! % list, tuple
%! assert (isequal (pycall_sympy__ ('return [1,2,3],'), {1, 2, 3}))
%! assert (isequal (pycall_sympy__ ('return (4,5),'), {4, 5}))
%! assert (isequal (pycall_sympy__ ('return (6,),'), {6,}))
%! assert (isequal (pycall_sympy__ ('return [],'), {}))

%!test
%! % dict
%! cmd = 'd = dict(); d["a"] = 6; d["b"] = 10; return d,';
%! d = pycall_sympy__ (cmd);
%! assert (d.a == 6 && d.b == 10)

%!test
%! r = pycall_sympy__ ('return 6');
%! assert (isequal (r, 6))

%!test
%! r = pycall_sympy__ ('return "Hi"');
%! assert (strcmp (r, 'Hi'))

%!test
%! % blank lines, lines with spaces
%! a = pycall_sympy__ ({ '', '', '     ', 'return 6', '   ', ''});
%! assert (isequal (a, 6))

%!test
%! % blank lines, strange comment lines
%! cmd = {'a = 1', '', '#', '', '#   ', '     #', 'a = a + 2', '  #', 'return a'};
%! a = pycall_sympy__ (cmd);
%! assert (isequal (a, 3))

%!test
%! % return empty string (was https://bugs.python.org/issue25270)
%! assert (isempty (pycall_sympy__ ('return ""')))

%!test
%! % return nothing (via an empty list)
%! % note distinct from 'return [],'
%! pycall_sympy__ ('return []')

%!test
%! % return nothing (because no return command)
%! pycall_sympy__ ('dummy = 1')

%!test
%! % return nothing (because no command)
%! pycall_sympy__ ('')

%!test
%! % return nothing (because no command)
%! pycall_sympy__ ({})

%!error <AttributeError>
%! % python exception while passing variables to python
%! % This tests the "INTERNAL_PYTHON_ERROR" path.
%! % FIXME: this is a very specialized test, relies on internal octsympy
%! % implementation details, and may need to be adjusted for changes.
%! disp('')
%! disp('**********************************************************************')
%! disp('')
%! disp('  Some deliberate AttributeError and ValueError may appear next')
%! disp('')
%! disp('v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v  v')
%! b = sym([], 'S.This_is_supposed_to_make_an_exception', [1 1], 'Test', 'Test', 'Test');
%! c = b + 1;
%!test
%! % ...and after the above test, the pipe should still work
%! a = pycall_sympy__ ('return _ins[0]*2', 3);
%! assert (isequal (a, 6))

%!test
%! % This command does not fail with native interface and '@pyobject'
%! s = warning ('off', 'OctSymPy:pythonic_no_convert');
%! try
%!   q = pycall_sympy__ ({'return type(int)'});
%! catch
%!   msg = lasterror.message;
%!   assert (~ isempty (regexp (msg, '.*does not know how to.*')))
%! end
%! warning (s)
%! disp('^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^')
%! disp('                  End of deliberate errors!')
%! disp('  You should see "46/46" or "Passes 46 out of 46 tests" below')
%! disp('**********************************************************************')
%! % ...and after the above test, the pipe should still work
%! a = pycall_sympy__ ('return _ins[0]*2', 3);
%! assert (isequal (a, 6))

%!test
%! % complex input
%! [A, B] = pycall_sympy__ ('z = 2*_ins[0]; return (z.real,z.imag)', 3+4i);
%! assert (A, 6)
%! assert (B, 8)

%!test
%! % complex output
%! z = pycall_sympy__ ('return 3+2j');
%! assert (z, 3+2i)

%!error <multirow char array>
%! s = char ('abc', 'defgh', '12345');
%! r = pycall_sympy__ ('return _ins[0]', s);

%!test
%! r = pycall_sympy__ ('return len(_ins[0])', '');
%! assert (r == 0)

%!test
%! % ensure the 0x0 matrix in octave gets mapped to the empty list in python
%! % @sym/subsasgn currently replies on this behaviour
%! assert (pycall_sympy__ ('return _ins[0] == []', []));

%!test
%! % ensure the empty string gets mapped to the empty string in python
%! assert (pycall_sympy__ ('return _ins[0] == ""', ''));

repo: octsympy