text:
function spm_eeg_prep_ui(callback)
% User interface for spm_eeg_prep function performing several tasks
% for preparation of converted MEEG data for further analysis
% FORMAT spm_eeg_prep_ui(callback)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_prep_ui.m 7169 2017-09-19 10:42:27Z vladimir $

spm('Pointer','Watch');

if ~nargin, callback = 'CreateMenu'; end

eval(callback);

spm('Pointer','Arrow');

end
%==========================================================================
% function CreateMenu
%==========================================================================
function CreateMenu

SVNrev = '$Rev: 7169 $';
spm('FnBanner', 'spm_eeg_prep_ui', SVNrev);
Finter = spm('FnUIsetup', 'M/EEG prepare', 0);

%-Draw top level menu
% ====== File ===================================
FileMenu = uimenu(Finter,'Label','File',...
    'Tag','EEGprepUI',...
    'HandleVisibility','on');

FileOpenMenu = uimenu(FileMenu, ...
    'Label','Open',...
    'Separator','off',...
    'Tag','EEGprepUI',...
    'HandleVisibility', 'on',...
    'Accelerator','O',...
    'Callback', 'spm_eeg_prep_ui(''FileOpenCB'')');

FileHeaderMenu = uimenu(FileMenu, ...
    'Label','Load header',...
    'Separator','off',...
    'Tag','EEGprepUI',...
    'HandleVisibility', 'on',...
    'Accelerator','H',...
    'Callback', 'spm_eeg_prep_ui(''FileHeaderCB'')');

FileImportMenu = uimenu(FileMenu, ...
    'Label','Import from workspace',...
    'Separator','off',...
    'Tag','EEGprepUI',...
    'HandleVisibility', 'on',...
    'Accelerator','I',...
    'Callback', 'spm_eeg_prep_ui(''FileImportCB'')');

FileSaveMenu = uimenu(FileMenu, ...
    'Label','Save',...
    'Separator','off',...
    'Tag','EEGprepUI',...
    'Enable','off',...
    'HandleVisibility', 'on',...
    'Accelerator','S',...
    'Callback', 'spm_eeg_prep_ui(''FileSaveCB'')');

FileExitMenu = uimenu(FileMenu, ...
    'Label','Quit',...
    'Separator','on',...
    'Tag','EEGprepUI',...
    'HandleVisibility', 'on',...
    'Accelerator','Q',...
    'Callback', 'spm_eeg_prep_ui(''FileExitCB'')');

% ====== Batch inputs ======================================

BatchInputsMenu = uimenu(Finter,'Label','Batch inputs',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on');

BInputsChannelsMenu = uimenu(BatchInputsMenu, 'Label', 'Channel selection',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''ChannelsCB'')');

BInputsTrialsMenu = uimenu(BatchInputsMenu, 'Label', 'Trial definition',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''TrialsCB'')');

BInputsEventsBIDSMenu = uimenu(BatchInputsMenu, ...
    'Label','Load events from BIDS',...
    'Separator','off',...
    'Enable','off',...
    'Tag','EEGprepUI',...
    'HandleVisibility', 'on',...
    'Callback', 'spm_eeg_prep_ui(''EventsBIDSCB'')');

BInputsEventsMenu = uimenu(BatchInputsMenu, 'Label', 'Event list',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''EventsCB'')');

BInputsMontageMenu = uimenu(BatchInputsMenu, 'Label', 'Montage',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on');

BInputsSortCondMenu = uimenu(BatchInputsMenu, 'Label', 'Sort conditions',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''SortCondCB'')');

MontageCustomMenu = uimenu(BInputsMontageMenu, 'Label', 'Custom Montage',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''MontageCB'')');

MontageRereferenceMenu = uimenu(BInputsMontageMenu, 'Label', 'Re-reference',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''RereferenceCB'')');

MontageROIMenu = uimenu(BInputsMontageMenu, 'Label', 'ROI',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''roiCB'')');

% ====== Channel types ===============================

ChanTypeMenu = uimenu(Finter,'Label','Channel types',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on');

chanTypes = {'EEG', 'EOG', 'ECG', 'EMG', 'LFP', 'PHYS', 'ILAM', 'SRC', 'Other'};

for i = 1:length(chanTypes)
    CTypesMenu(i) = uimenu(ChanTypeMenu, 'Label', chanTypes{i},...
        'Tag','EEGprepUI',...
        'Enable', 'on', ...
        'HandleVisibility','on',...
        'Callback', 'spm_eeg_prep_ui(''ChanTypeCB'')');
end

CTypesRef2MEGMenu = uimenu(ChanTypeMenu, 'Label', 'MEGREF=>MEG',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on',...
    'Separator', 'on',...
    'Callback', 'spm_eeg_prep_ui(''MEGChanTypeCB'')');

CTypesDefaultMenu = uimenu(ChanTypeMenu, 'Label', 'Default',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on',...
    'Callback', 'spm_eeg_prep_ui(''ChanTypeDefaultCB'')');

CTypesBIDSMenu = uimenu(ChanTypeMenu, 'Label', 'From BIDS',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''ChanTypeBIDSCB'')');

CTypesReviewMenu = uimenu(ChanTypeMenu, 'Label', 'Review',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''ChanTypeCB'')');

% ====== Sensors ===================================

Coor3DMenu = uimenu(Finter,'Label','Sensors',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on');

LoadMEGSensMenu = uimenu(Coor3DMenu, 'Label', 'Load MEG sensors',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''LoadMEGSensCB'')');

HeadshapeMenu = uimenu(Coor3DMenu, 'Label', 'Load MEG Fiducials/Headshape',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''HeadshapeCB'')');

LoadEEGSensMenu = uimenu(Coor3DMenu, 'Label', 'Load EEG sensors',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'Separator', 'on', ...
    'HandleVisibility','on');

LoadEEGSensTemplateMenu = uimenu(LoadEEGSensMenu, 'Label', 'Assign default',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''LoadEEGSensTemplateCB'')');

LoadEEGSensMatMenu = uimenu(LoadEEGSensMenu, 'Label', 'From *.mat file',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''LoadEEGSensCB'')');

LoadEEGSensOtherMenu = uimenu(LoadEEGSensMenu, 'Label', 'Convert locations file',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''LoadEEGSensCB'')');

ReferenceMenu = uimenu(Coor3DMenu, 'Label', 'Define EEG referencing',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on');

ReferenceSelectMenu = uimenu(ReferenceMenu, 'Label', 'Select reference sensors',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''ReferenceCB'')');

ReferenceMontageMenu = uimenu(ReferenceMenu, 'Label', 'Specify montage matrix',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''ReferenceCB'')');

CoregisterMenu = uimenu(Coor3DMenu, 'Label', 'Coregister',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on', ...
    'Callback', 'spm_eeg_prep_ui(''CoregisterCB'')');

% ====== 2D projection ===================================

Coor2DMenu = uimenu(Finter, 'Label','2D projection',...
    'Tag','EEGprepUI',...
    'Enable', 'off', ...
    'HandleVisibility','on');

EditMEGMenu = uimenu(Coor2DMenu, 'Label', 'Edit existing MEG',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''EditExistingCoor2DCB'')');

EditEEGMenu = uimenu(Coor2DMenu, 'Label', 'Edit existing EEG',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''EditExistingCoor2DCB'')');

LoadTemplateMenu = uimenu(Coor2DMenu, 'Label', 'Load template',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on', ...
    'Callback', 'spm_eeg_prep_ui(''LoadTemplateCB'')');

SaveTemplateMenu = uimenu(Coor2DMenu, 'Label', 'Save template',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''SaveTemplateCB'')');

Project3DEEGMenu = uimenu(Coor2DMenu, 'Label', 'Project 3D (EEG)',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on', ...
    'Callback', 'spm_eeg_prep_ui(''Project3DCB'')');

Project3DMEGMenu = uimenu(Coor2DMenu, 'Label', 'Project 3D (MEG)',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''Project3DCB'')');

AddCoor2DMenu = uimenu(Coor2DMenu, 'Label', 'Add sensor',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on', ...
    'Callback', 'spm_eeg_prep_ui(''AddCoor2DCB'')');

DeleteCoor2DMenu = uimenu(Coor2DMenu, 'Label', 'Delete sensor',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''DeleteCoor2DCB'')');

UndoMoveCoor2DMenu = uimenu(Coor2DMenu, 'Label', 'Undo move',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''UndoMoveCoor2DCB'')');

ApplyCoor2DMenu = uimenu(Coor2DMenu, 'Label', 'Apply',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Separator', 'on', ...
    'Callback', 'spm_eeg_prep_ui(''ApplyCoor2DCB'')');

Clear2DMenu = uimenu(Coor2DMenu, 'Label', 'Clear',...
    'Tag','EEGprepUI',...
    'Enable', 'on', ...
    'HandleVisibility','on',...
    'Callback', 'spm_eeg_prep_ui(''Clear2DCB'')');
end

%==========================================================================
% function FileOpenCB
%==========================================================================
function FileOpenCB

D = spm_eeg_load;
setD(D);

update_menu;

end

%==========================================================================
% function FileHeaderCB
%==========================================================================
function FileHeaderCB

S = [];
S.mode = 'header';
[S.dataset, sts] = spm_select(1, '.*', 'Select M/EEG data file');
if ~sts, return; end
D = spm_eeg_convert(S);
setD(D);

update_menu;

end

%==========================================================================
% function FileImportCB
%==========================================================================
function FileImportCB

var = evalin('base', 'whos');
var = {var(ismember({var.class}, {'single', 'double'})).name};

[selection, ok]= listdlg('ListString', var, 'SelectionMode', 'single' ,'Name', 'Select variable' , 'ListSize', [400 300]);

if ok
    data = full(double(evalin('base', var{selection})));
    dim = size(data);
    ndim = length(dim);

    if length(dim)>3
        error('Data of up to 3 dimensions are supported.');
    end

    prompt = sprintf('%d|', dim);
    prompt(end) = [];

    chandim = spm_input('Number of channels:', 1, prompt, 1:length(dim));

    nchan = dim(chandim);

    otherdim = setdiff(1:ndim, chandim);

    if ndim == 3
        prompt = sprintf('%d|', dim(otherdim));
        prompt(end) = [];
        trialdim = spm_input('Number of trials:','+1', prompt, otherdim);

        ntrials = dim(trialdim);
    else
        trialdim = [];
        ntrials  = 1;
    end

    sampledim = setdiff(1:ndim, [trialdim, chandim]);
    nsamples  = dim(sampledim);

    data    = permute(data, [chandim, sampledim, trialdim]);

    fs      = spm_input('Sampling rate (Hz):', '+1','r');
    t_onset = spm_input('Time axis onset (ms):', '+1','r');

    fname   = spm_input('Output file name:', '+1','s');

    time    = (0:(nsamples-1))/fs - 1e-3*t_onset;

    for i = 1:ntrials
        ftdata.trial{i} = data(:, :, i);
        ftdata.time{i}  = time;
    end

    for i = 1:nchan
        ftdata.label{i, 1} = sprintf('Ch%d', i);
    end

    D = spm_eeg_ft2spm(ftdata, fname);

    spm_eeg_review(D);
    spm_eeg_review_callbacks('edit', 'prep');
end
end
%==========================================================================
% function FileSaveCB
%==========================================================================
function FileSaveCB

D = getD;
if ~isempty(D)
    D.save;
end

update_menu;

end

%==========================================================================
% function FileExitCB
%==========================================================================
function FileExitCB

spm_figure('Clear','Interactive');
spm('FigName','M/EEG prepare: done');

end

%==========================================================================
% function ChannelsCB
%==========================================================================
function ChannelsCB

D = getD;

[selection, ok]= listdlg('ListString', D.chanlabels, 'SelectionMode', 'multiple' ,'Name', 'Select channels' , 'ListSize', [400 300]);
if ~ok, return; end

[chanfilename, chanpathname] = uiputfile( ...
    {'*.mat', 'MATLAB File (*.mat)'}, 'Save channel selection as');

label = D.chanlabels(selection);

save(fullfile(chanpathname, chanfilename), 'label', spm_get_defaults('mat.format'));

end

%==========================================================================
% function TrialsCB
%==========================================================================
function TrialsCB

D = getD;

S = [];
S.D = D;
S.save = 1;
spm_eeg_definetrial(S);

CreateMenu;

setD(D);

update_menu;

end

%==========================================================================
% function EventsBIDSCB
%==========================================================================
function EventsBIDSCB

D = getD;

S = [];
S.task = 'loadbidsevents';
S.D = D;
[S.filename, sts] = spm_select(1, '.*_events.tsv$', 'Select BIDS tsv file');

if sts
    S.replace = spm_input('Replace existing','1','replace|add',[1 0], 1);
end

D = spm_eeg_prep(S);

setD(D);

update_menu;

end

%==========================================================================
% function EventsCB
%==========================================================================
function EventsCB

D = getD;

ev = D.events;

if iscell(ev)
    sev = ev{1};
    for i = 2:numel(ev)
        sev = spm_cat_struct(sev, ev{i});
    end
else
    sev = ev;
end

selected = spm_eeg_select_event_ui(sev);

events = [];
for i = 1:size(selected, 1)
    events(i).type  = selected{i, 1};
    events(i).value = selected{i, 2};
end

if ~isempty(events)
    [filename, pathname] = uiputfile('*.mat', 'Save event list as');

    if ~isequal(filename, 0)
        save(fullfile(pathname, filename), 'events', spm_get_defaults('mat.format'));
    end
end
end

%==========================================================================
% function MontageCB
%==========================================================================
function MontageCB

D = getD;

label = D.chanlabels(D.indchantype({'EEG', 'EOG', 'ECG', 'EMG'}));

montage = [];
montage.labelorg = label;
montage.labelnew = label;
montage.tra = eye(numel(label));

spm_eeg_montage_ui(montage);

end

%==========================================================================
% function RereferenceCB
%==========================================================================
function RereferenceCB

D = getD;

% Get indices for just EEG channels and remove any bad channels
%--------------------------------------------------------------------------
eegchan  = D.indchantype('EEG');
goodind  = D.indchantype('EEG', 'GOOD');

if isempty(goodind)
    error('No good EEG channels.')
end

%-Get reference channel indices
%--------------------------------------------------------------------------
[selection, ok]= listdlg('ListString', D.chanlabels(goodind), 'SelectionMode', 'multiple' ,...
    'Name', 'Select reference channels' , 'ListSize', [400 300]);
if ~ok
    return;
end

refind  = goodind(selection);
badind  = D.indchantype('EEG', 'BAD');

goodind = find(ismember(eegchan, goodind));
badind  = find(ismember(eegchan,  badind));
refind  = find(ismember(eegchan,  refind));

tra                 = eye(length(eegchan));
tra(goodind,refind) = tra(goodind,refind) - 1/length(refind);
tra(badind, refind) = tra(badind,refind)  - 1/length(refind);

montage          = [];
montage.labelorg = D.chanlabels(eegchan);
montage.labelnew = D.chanlabels(eegchan);
montage.tra      = tra;

[filename, pathname] = uiputfile('*.mat', 'Save montage as');

if ~isequal(filename, 0)
    save(fullfile(pathname, filename), 'montage', spm_get_defaults('mat.format'));
end

end

%==========================================================================
% function roiCB
%==========================================================================
function roiCB

D = getD;

[modality, chanind]  = spm_eeg_modality_ui(D, 0, 1);

chanind = setdiff(chanind, D.badchannels);

if isempty(chanind)
    error(['No good ' modality ' channels.']);
end

montage = [];
montage.labelorg = D.chanlabels(chanind)';
montage.labelnew = {};
montage.tra      = [];
while 1
    roilabel = spm_input('ROI Label:', '+1', 's');

    [selection, ok]= listdlg('ListString', D.chanlabels(chanind), 'SelectionMode', 'multiple' ,...
        'Name', 'Select reference channels' , 'ListSize', [400 300]);
    if ok
       montage.labelnew{end+1, 1} = roilabel;
       tra = zeros(1, length(chanind));
       tra(selection) = 1;
       tra = tra./sum(tra);
       montage.tra = [montage.tra; tra];
    end

    if spm_input('Add another?','+1','yes|no',[0, 1], 0);
        break;
    end
end

[filename, pathname] = uiputfile('*.mat', 'Save montage as');

if ~isequal(filename, 0)
    save(fullfile(pathname, filename), 'montage', spm_get_defaults('mat.format'));
end

update_menu;

end

%==========================================================================
% function SortCondCB
%==========================================================================
function SortCondCB

D = getD;

S   = [];
S.task = 'sortconditions';
S.D = D;
oldcondlist = D.condlist;
S.condlist = cell(size(oldcondlist));
for i = 1:D.nconditions
    str = sprintf('%s|', oldcondlist{:});
    str = str(1:(end-1));

    ind = spm_input(['Select condition ' num2str(i)], 1, 'm', str, 1:numel(oldcondlist));
    S.condlist(i) = oldcondlist(ind);
    oldcondlist(ind) = [];
end

[filename, pathname] = uiputfile('*.mat', 'Save conditions list as');

if ~isequal(filename, 0)
    condlist = S.condlist;

    save(fullfile(pathname, filename), 'condlist', spm_get_defaults('mat.format'));
end

D = spm_eeg_prep(S);

setD(D);
update_menu;

end

%==========================================================================
% function ChanTypeCB
%==========================================================================
function ChanTypeCB

type = get(gcbo, 'Label');

D = getD;

if ~isempty(D)
    chanlist ={};
    for i = 1:D.nchannels
        if strncmp(D.chantype(i), 'MEG', 3) || strncmp(D.chantype(i), 'REF', 3)
            chanlist{i} = [num2str(i) '    Label:    ' D.chanlabels(i) '    Type:    ' D.chantype(i) , ' (nonmodifiable)'];
        else
            chanlist{i} = [num2str(i) '    Label:    ' D.chanlabels(i) '    Type:    ' D.chantype(i)];
        end

        chanlist{i} = [chanlist{i}{:}];
    end

    if strcmpi(type, 'review')
        listdlg('ListString', chanlist, 'SelectionMode', 'single', 'Name', 'Review channels', 'ListSize', [400 300]);
        return
    else

        [selection,ok]= listdlg('ListString', chanlist, 'SelectionMode', 'multiple',...
            'InitialValue', strmatch(type, D.chantype) ,'Name', ['Set type to ' type], 'ListSize', [400 300]);

        selection(strmatch('MEG', chantype(D, selection))) = [];

        if ok && ~isempty(selection)
            S.task = 'settype';
            S.D = D;
            S.ind = selection;
            S.type = type;
            D = spm_eeg_prep(S);
            setD(D);
        end
    end
end

update_menu;

end
%==========================================================================
% function MEGChanTypeCB
%==========================================================================
function MEGChanTypeCB

S = [];
S.D = getD;
S.task = 'settype';

switch get(gcbo, 'Label')
    case 'MEGREF=>MEG'
        dictionary = {
            'REFMAG',   'MEGMAG';
            'REFGRAD',  'MEGGRAD';
            };

        ind = spm_match_str(S.D.chantype, dictionary(:,1));

        grad = S.D.sensors('meg');
        if ~isempty(grad)
            % Under some montages only subset of the reference sensors are
            % in the grad
            [junk, sel] = intersect(S.D.chanlabels(ind), grad.label);
            ind = ind(sel);
        end

        S.ind = ind;

        [sel1, sel2] = spm_match_str(S.D.chantype(S.ind), dictionary(:, 1));

        S.type = dictionary(sel2, 2);

        D = spm_eeg_prep(S);
end

setD(D);
update_menu;

end
%==========================================================================
% function ChanTypeDefaultCB
%==========================================================================
function ChanTypeDefaultCB

S.D    = getD;
S.task = 'defaulttype';
D      = spm_eeg_prep(S);

setD(D);
update_menu;

end

%==========================================================================
% function ChanTypeBIDSCB
%==========================================================================
function ChanTypeBIDSCB

S = [];
S.D    = getD;
S.task = 'bidschantype';
[S.filename, sts] = spm_select(1, '.*_channels.tsv$', 'Select BIDS tsv file');

if sts
    D      = spm_eeg_prep(S);

    setD(D);
    update_menu;
end

end

%==========================================================================
% function LoadEEGSensTemplateCB
%==========================================================================
function LoadEEGSensTemplateCB

S.D    = getD;
S.task = 'defaulteegsens';

if strcmp(S.D.modality(1, 0), 'Multimodal')
    fid = fiducials(S.D);
    if ~isempty(fid)
        lblfid = fid.fid.label;

        S.regfid = match_fiducials({'nas'; 'lpa'; 'rpa'}, lblfid);
        S.regfid(:, 2) = {'spmnas'; 'spmlpa'; 'spmrpa'};
    else
        warndlg(strvcat('Could not match EEG fiducials for multimodal dataset.', ...
            '           EEG coregistration might fail.'));
    end
end

D = spm_eeg_prep(S);
setD(D);
update_menu;

end
%==========================================================================
% function LoadEEGSensCB
%==========================================================================
function LoadMEGSensCB

D = getD;

S = [];
S.task = 'loadmegsens';
S.D = D;
[S.source,sts] = spm_select(1, '\.*', 'Select a raw MEG data file');
if ~sts, return; end

D = spm_eeg_prep(S);

setD(D);

update_menu;

end
%==========================================================================
% function LoadEEGSensCB
%==========================================================================
function LoadEEGSensCB

D = getD;

switch get(gcbo, 'Label')
    case 'From *.mat file'
        [S.sensfile, sts] = spm_select(1,'mat','Select EEG sensors file');
        if ~sts, return, end
        S.source = 'mat';
        [S.headshapefile, sts] = spm_select(1,'mat','Select EEG fiducials file');
        if ~sts, return, end
        S.fidlabel = spm_input('Fiducial labels:', '+1', 's', 'nas lpa rpa');
    case 'Convert locations file'
        [S.sensfile, sts] = spm_select(1, '.*', 'Select locations file');
        if ~sts, return, end
        S.source = 'locfile';
end

if strcmp(D.modality(1, 0), 'Multimodal')
    if ~isempty(D.fiducials)
        S.regfid = {};
        if strcmp(S.source, 'mat')
            fidlabel = S.fidlabel;
            lblshape = {};
            fidnum = 0;
            while ~all(isspace(fidlabel))
                fidnum = fidnum+1;
                [lblshape{fidnum},fidlabel] = strtok(fidlabel);
            end
            if (fidnum < 3)
                error('At least 3 labeled fiducials are necessary');
            end
        else
            shape = ft_read_headshape(S.sensfile);
            lblshape = shape.fid.label;
        end

        fid = fiducials(D);
        lblfid = fid.fid.label;

        S.regfid = match_fiducials(lblshape, lblfid);
    else
        warndlg(strvcat('Could not match EEG fiducials for multimodal dataset.', ...
            '           EEG coregistration might fail.'));
    end
end

S.D    = D;
S.task = 'loadeegsens';
D      = spm_eeg_prep(S);

setD(D);

update_menu;

end
%==========================================================================
% function HeadshapeCB
%==========================================================================
function HeadshapeCB

S = [];
S.D = getD;
S.task = 'headshape';

[S.headshapefile, sts] = spm_select(1, '.*', 'Select fiducials/headshape file');
if ~sts, return, end
S.source = 'convert';

shape = ft_read_headshape(S.headshapefile);
lblshape = shape.fid.label;

fid = fiducials(S.D);

if ~isempty(fid)
    lblfid = fid.fid.label;
    S.regfid = match_fiducials(lblshape, lblfid);
end

D = spm_eeg_prep(S);

setD(D);

update_menu;

end
%==========================================================================
% function CoregisterCB
%==========================================================================
function ReferenceCB

D = getD;

S = [];
S.D = D;
S.task = 'sens2chan';
switch get(gcbo, 'Label')
    case 'Select reference sensors'
        elec = D.sensors('EEG');

        [refind, ok]= listdlg('ListString', elec.label, 'SelectionMode', 'multiple' ,'Name', 'Select reference electrodes' , 'ListSize', [400 300]);
        if ~ok
            return;
        end

        S.refelec = elec.label(refind);
    case 'Specify montage matrix'

        sens  = D.sensors('EEG');
        label = D.chanlabels(D.indchantype('EEG'));

        [sel1, sel2] = spm_match_str(label, sens.label);

        montage = [];
        montage.labelorg = sens.label;
        montage.labelnew = label;
        montage.tra = zeros(numel(label), numel(sens.label));
        montage.tra(sub2ind(size(montage.tra), sel1, sel2)) = 1;

        S.montage = spm_eeg_montage_ui(montage);
end

D = spm_eeg_prep(S);

setD(D);

update_menu;

end
%==========================================================================
% function CoregisterCB
%==========================================================================
function CoregisterCB

S = [];
S.D = getD;
S.task = 'coregister';

D = spm_eeg_prep(S);

% Bring the menu back
spm_eeg_prep_ui;

setD(D);

update_menu;

end
%==========================================================================
% function EditExistingCoor2DCB
%==========================================================================
function EditExistingCoor2DCB

D = getD;

switch get(gcbo, 'Label')
    case 'Edit existing MEG'
        xy = D.coor2D('MEG');
        label = D.chanlabels(strmatch('MEG', D.chantype));
    case 'Edit existing EEG'
        xy = D.coor2D('EEG');
        label = D.chanlabels(strmatch('EEG', D.chantype, 'exact'));
end

plot_sensors2D(xy, label);

update_menu;

end
%==========================================================================
% function LoadTemplateCB
%==========================================================================
function LoadTemplateCB

[sensorfile, sts] = spm_select(1, 'mat', 'Select sensor template file', ...
    [], fullfile(spm('dir'), 'EEGtemplates'));
if ~sts, return, end

template = load(sensorfile);

if isfield(template, 'Cnames') && isfield(template, 'Cpos')
    plot_sensors2D(template.Cpos, template.Cnames);
end

update_menu;

end
%==========================================================================
% function SaveTemplateCB
%==========================================================================
function SaveTemplateCB

handles=getHandles;

Cnames = handles.label;
Cpos = handles.xy;
Rxy = 1.5;
Nchannels = length(Cnames);

[filename, pathname] = uiputfile('*.mat', 'Save channel template as');

save(fullfile(pathname, filename), 'Cnames', 'Cpos', 'Rxy', 'Nchannels', spm_get_defaults('mat.format'));

end
%==========================================================================
% function Project3DCB
%==========================================================================
function Project3DCB

D = getD;

switch get(gcbo, 'Label')
    case 'Project 3D (EEG)'
        modality = 'EEG';
        chanind  = D.indchantype('EEG');
    case 'Project 3D (MEG)'
        modality = 'MEG';
        chanind  = D.indchantype('MEGANY');
end

if ~isfield(D, 'val')
    D.val = 1;
end

if isfield(D, 'inv') && isfield(D.inv{D.val}, 'datareg')
    datareg = D.inv{D.val}.datareg;
    ind     = strmatch(modality, {datareg(:).modality}, 'exact');
    sens    = datareg(ind).sensors;
else
    sens    = D.sensors(modality);
end

[xy, label] = spm_eeg_project3D(sens, modality);

[sel1, sel2] = spm_match_str(D.chanlabels(chanind), label);

plot_sensors2D(xy(:, sel2), label(sel2));

update_menu;

end
%==========================================================================
% function AddCoor2DCB
%==========================================================================
function AddCoor2DCB

newlabel = spm_input('Label?', '+1', 's');

if isempty(newlabel)
    return;
end

coord = spm_input('Coordinates [x y]', '+1', 'r', '0.5 0.5', 2);

handles = getHandles;

if ~isfield(handles, 'xy')
    handles.xy = [];
end

if ~isfield(handles, 'xy')
    handles.xy = [];
end

if ~isfield(handles, 'label')
    handles.label = {};
end

plot_sensors2D([handles.xy coord(:)], ...
    [handles.label newlabel]);

update_menu;

end
%==========================================================================
% function ApplyCoor2DCB
%==========================================================================
function ApplyCoor2DCB

handles = getHandles;
D = getD;

S = [];
S.task = 'setcoor2d';
S.D = D;
S.xy = handles.xy;
S.label = handles.label;

D = spm_eeg_prep(S);

setD(D);

update_menu;

end
%==========================================================================
% function update_menu
%==========================================================================
function update_menu

Finter = spm_figure('GetWin','Interactive');
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'File'), 'Enable', 'on');

IsEEG = 'off';
IsMEG = 'off';
IsEpochable = 'off';
HasEvents   = 'off';
IsEpoched   = 'off';
HasPlanar   = 'off';
HasSensors  = 'off';
HasSensorsEEG = 'off';
ReferenceSelectable = 'off';
HasSensorsMEG = 'off';
HasChannelsMEGREF = 'off';
HasFiducials = 'off';
HasDefaultLocs = 'off';
HasHistory = 'off';
if isa(get(Finter, 'UserData'), 'meeg')
    Dloaded = 'on';

    D = getD;

    if ~isempty(D.indchantype('EEG'))
        IsEEG = 'on';
    end

    if ~isempty(D.indchantype('MEG'));
        IsMEG = 'on';
    end

    if ~isempty(D.events)
        HasEvents = 'on';
    end

    if isequal(D.type, 'continuous') && ~isempty(D.events);
        IsEpochable = 'on';
    end

    if ~isequal(D.type, 'continuous')
        IsEpoched = 'on';
    end

    [res, list] = modality(D, 1, 1);
    if ismember('MEGPLANAR', list);
        HasPlanar = 'on';
    end

    if ~isempty(D.indchantype('REF'));
        HasChannelsMEGREF = 'on';
    end

    if ~isempty(D.sensors('EEG')) || ~isempty(D.sensors('MEG'))
        HasSensors = 'on';
    end

    if ~isempty(D.sensors('EEG'))
        HasSensorsEEG = 'on';

        elec = D.sensors('EEG');
        if ~isfield(elec, 'tra') &&  (size(elec.elecpos, 1) == numel(elec.label))
            ReferenceSelectable = 'on';
        end
    end

    if  ~isempty(D.sensors('MEG'))
        HasSensorsMEG = 'on';
    end

    if  ~isempty(D.fiducials)
        HasFiducials = 'on';
    end

    template_sfp = dir(fullfile(spm('dir'), 'EEGtemplates', '*.sfp'));
    template_sfp = {template_sfp.name};
    ind = strmatch([ft_senstype(D.chanlabels(D.indchantype('EEG'))) '.sfp'], template_sfp, 'exact');

    if ~isempty(ind) || ft_senstype(D.chanlabels(D.indchantype('EEG')), 'ext1020')
        HasDefaultLocs = 'on';
    end

    if  ~isempty(D.history)
        HasHistory = 'on';
    end

else
    Dloaded = 'off';
end

handles = getHandles;

IsTemplate = 'off';
IsSelected = 'off';
IsMoved = 'off';
if ~isempty(handles)
    if isfield(handles, 'xy') && size(handles.xy, 1)>0
        IsTemplate = 'on';
    end

    if isfield(handles, 'labelSelected') && ~isempty(handles.labelSelected)
        IsSelected = 'on';
    end

    if isfield(handles, 'lastMoved')
        isMoved = 'on';
    end
end

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Save'), 'Enable', 'on');
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Load events from BIDS'), 'Enable', Dloaded);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Batch inputs'), 'Enable', Dloaded);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Trial definition'), 'Enable', IsEpochable);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Event list'), 'Enable', HasEvents);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Montage'), 'Enable', Dloaded);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Re-reference'), 'Enable', IsEEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'ROI'), 'Enable', Dloaded);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Sort conditions'), 'Enable', IsEpoched);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Channel types'), 'Enable', Dloaded);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Sensors'), 'Enable', Dloaded);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', '2D projection'), 'Enable', Dloaded);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'MEGREF=>MEG'), 'Enable', HasChannelsMEGREF);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Assign default'), 'Enable', HasDefaultLocs);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Load EEG sensors'), 'Enable', IsEEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Load MEG sensors'), 'Enable', IsMEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Load MEG Fiducials/Headshape'), 'Enable', HasSensorsMEG);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Headshape'), 'Enable', HasSensorsMEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Define EEG referencing'), 'Enable', HasSensorsEEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Select reference sensors'), 'Enable', ReferenceSelectable);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Coregister'), 'Enable', HasSensors);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Edit existing EEG'), 'Enable', IsEEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Edit existing MEG'), 'Enable', IsMEG);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Project 3D (EEG)'), 'Enable', HasSensorsEEG);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Project 3D (MEG)'), 'Enable', HasSensorsMEG);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Delete sensor'), 'Enable', IsSelected);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Undo move'), 'Enable', IsMoved);

set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Apply'), 'Enable', IsTemplate);
set(findobj(Finter,'Tag','EEGprepUI', 'Label', 'Clear'), 'Enable', IsTemplate);

delete(setdiff(findobj(Finter), [Finter; findobj(Finter,'Tag','EEGprepUI')]));

if strcmp(Dloaded, 'on') && isfield(D,'PSD') && D.PSD == 1
    try
        hc = get(Finter,'children');
        hc = findobj(hc,'flat','type','uimenu');
        hc = findobj(hc,'flat','label','File');
        delete(hc)
    end
    uicontrol(Finter,...
        'style','pushbutton','string','OK',...
        'callback','spm_eeg_review_callbacks(''get'',''prep'')',...
        'tooltipstring','Send changes to ''SPM Graphics'' window',...
        'BusyAction','cancel',...
        'Interruptible','off',...
        'Tag','EEGprepUI');
end

figure(Finter);

end
%==========================================================================
% function getD
%==========================================================================
function D = getD

Finter = spm_figure('GetWin','Interactive');
D = get(Finter, 'UserData');
if ~isa(D, 'meeg')
    D = [];
end

end
%==========================================================================
% function setD
%==========================================================================
function setD(D)
Finter = spm_figure('GetWin','Interactive');
set(Finter, 'UserData', D);

end
%==========================================================================
% function getHandles
%==========================================================================
function handles = getHandles
Fgraph = spm_figure('GetWin','Graphics');
handles = get(Fgraph, 'UserData');

end
%==========================================================================
% function setHandles
%==========================================================================
function setHandles(handles)
Fgraph = spm_figure('GetWin','Graphics');
set(Fgraph, 'UserData', handles);

end
%==========================================================================
% function plot_sensors2D
%==========================================================================
function plot_sensors2D(xy, label)

Fgraph = spm_figure('GetWin','Graphics');

spm_clf(Fgraph);

handles = [];

if ~isempty(xy)
    if size(xy, 1) ~= 2
        xy = xy';
    end

    xy(xy < 0.05) = 0.05;
    xy(xy > 0.95) = 0.95;

    handles.h_lbl=text(xy(1,:), xy(2, :),strvcat(label),...
        'FontSize', 9,...
        'Color','r',...
        'FontWeight','bold');

    set(handles.h_lbl, 'ButtonDownFcn', 'spm_eeg_prep_ui(''LabelClickCB'')');

    hold on

    handles.h_el = [];
    for i=1:size(xy, 2)
        handles.h_el(i) = plot(xy(1,i), xy(2,i), 'or');
    end

    set(handles.h_el,'MarkerFaceColor','r','MarkerSize', 2,'MarkerEdgeColor','k');

    handles.TemplateFrame = ...
        plot([0.05 0.05 0.95 0.95 0.05], [0.05 0.95 0.95 0.05 0.05], 'k-');

    axis off;

end

handles.xy = xy;
handles.label = label(:)';

setHandles(handles);

update_menu;

end
%==========================================================================
% function DeleteCoor2DCB
%==========================================================================
function DeleteCoor2DCB

handles = getHandles;
graph = spm_figure('GetWin','Graphics');

if isfield(handles, 'labelSelected') && ~isempty(handles.labelSelected)
    set(graph, 'WindowButtonDownFcn', '');

    label=get(handles.labelSelected, 'String');
    ind=strmatch(label, handles.label, 'exact');

    delete([handles.labelSelected handles.pointSelected]);

    handles.xy(:, ind)=[];
    handles.label(ind) = [];

    plot_sensors2D(handles.xy, handles.label)
end

end
%==========================================================================
% function UndoMoveCoor2DCB
%==========================================================================
function UndoMoveCoor2DCB

handles = getHandles;

if isfield(handles, 'lastMoved')
    label = get(handles.lastMoved(end).label, 'String');
    ind = strmatch(label, handles.label, 'exact');
    handles.xy(:, ind) = handles.lastMoved(end).coords(:);

    set(handles.lastMoved(end).point, 'XData', handles.lastMoved(end).coords(1));
    set(handles.lastMoved(end).point, 'YData', handles.lastMoved(end).coords(2));
    set(handles.lastMoved(end).label, 'Position', handles.lastMoved(end).coords);

    if length(handles.lastMoved)>1
        handles.lastMoved = handles.lastMoved(1:(end-1));
    else
        handles = rmfield(handles, 'lastMoved');
    end

    setHandles(handles);
    update_menu;
end

end
%==========================================================================
% function LabelClickCB
%==========================================================================
function LabelClickCB

handles=getHandles;
Fgraph = spm_figure('GetWin','Graphics');

if isfield(handles, 'labelSelected') && ~isempty(handles.labelSelected)
    if handles.labelSelected == gcbo
        set(handles.labelSelected, 'Color', 'r');
        set(handles.pointSelected,'MarkerFaceColor', 'r');
        set(Fgraph, 'WindowButtonDownFcn', '');
    else
        handles.pointSelected=[];
        handles.labelSelected=[];
    end
else
    set(Fgraph, 'WindowButtonDownFcn', 'spm_eeg_prep_ui(''LabelMoveCB'')');

    coords = get(gcbo, 'Position');

    handles.labelSelected=gcbo;
    handles.pointSelected=findobj(gca, 'Type', 'line',...
        'XData', coords(1), 'YData', coords(2));

    set(handles.labelSelected, 'Color', 'g');
    set(handles.pointSelected,'MarkerFaceColor', 'g');
end

setHandles(handles);

update_menu;

end
%==========================================================================
% function LabelMoveCB
%==========================================================================
function LabelMoveCB

handles = getHandles;
Fgraph = spm_figure('GetWin','Graphics');

coords=mean(get(gca, 'CurrentPoint'));

coords(coords < 0.05) = 0.05;
coords(coords > 0.95) = 0.95;

set(handles.pointSelected, 'XData', coords(1));
set(handles.pointSelected, 'YData', coords(2));

set(handles.labelSelected, 'Position', coords);
set(handles.labelSelected, 'Color', 'r');
set(handles.pointSelected,'MarkerFaceColor','r','MarkerSize',2,'MarkerEdgeColor','k');
set(Fgraph, 'WindowButtonDownFcn', '');
set(Fgraph, 'WindowButtonMotionFcn', 'spm_eeg_prep_ui(''CancelMoveCB'')');

labelind=strmatch(get(handles.labelSelected, 'String'), handles.label);

if isfield(handles, 'lastMoved')
    handles.lastMoved(end+1).point = handles.pointSelected;
    handles.lastMoved(end).label = handles.labelSelected;
    handles.lastMoved(end).coords = handles.xy(:, labelind);
else
    handles.lastMoved.point = handles.pointSelected;
    handles.lastMoved.label = handles.labelSelected;
    handles.lastMoved.coords = handles.xy(:, labelind);
end

handles.xy(:, labelind) = coords(1:2)';

setHandles(handles);

update_menu;

end
%==========================================================================
% function CancelMoveCB
%==========================================================================
function CancelMoveCB

Fgraph = spm_figure('GetWin','Graphics');
handles = getHandles;

handles.pointSelected=[];
handles.labelSelected=[];
set(Fgraph, 'WindowButtonMotionFcn', '');

setHandles(handles);

update_menu;

end
%==========================================================================
% function Clear2DCB
%==========================================================================
function Clear2DCB

plot_sensors2D([], {});

update_menu;

end
%==========================================================================
% function match_fiducials
%==========================================================================
function regfid = match_fiducials(lblshape, lblfid)

if numel(intersect(upper(lblshape), upper(lblfid))) < 3
    if numel(lblshape)<3 || numel(lblfid)<3
        warndlg('3 fiducials are required');
        return;
    else
        regfid = {};
        for i = 1:length(lblfid)
            [selection,ok]= listdlg('ListString',lblshape, 'SelectionMode', 'single',...
                'InitialValue', strmatch(upper(lblfid{i}), upper(lblshape)), ...
                'Name', ['Select matching fiducial for ' lblfid{i}], 'ListSize', [400 300]);

            if ~ok
                continue
            end

            regfid = [regfid; [lblfid(i) lblshape(selection)]];
        end

        if size(regfid, 1) < 3
            warndlg('3 fiducials are required to load headshape');
            return;
        end
    end
else
    [sel1, sel2] = spm_match_str(upper(lblfid), upper(lblshape));
    lblfid = lblfid(sel1);
    lblshape = lblshape(sel2);
    regfid = [lblfid(:) lblshape(:)];
end
end

author: spm


text:
function uscita = deg2rad(ingresso)
uscita = 2*pi*ingresso/360;

author: Sable

repo: mcbench-benchmarks


text:
function y = variablereplace(y,oldVar,newVar)

var = y.lmi_variables;

[~,pos] = ismember(oldVar,var);
if any(pos)
    index = find(pos);
    pos = pos(index);
    y.lmi_variables(pos) = newVar(index);

    [var,pos] = sort(y.lmi_variables);
    y.lmi_variables = var;
    y.basis = [y.basis(:,1) y.basis(:,1+pos)];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function mObject = setAnchorFromRawAnchor(mObject,rawAnchor);
%   Set anchor points using raw anchor information
%   mObject = setAnchorFromRawAnchor(mObject,rawAnchor);

%   Designed and coded by Hideki Kawahara
%   27/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara

TIMINGmARGIN = 10; % threshould for merging location
[dm1,indsrt] = sort(rawAnchor(:,1));
sortedAnchor = rawAnchor(indsrt,1);
sortedFrequency = rawAnchor(indsrt,2);
indexNumber = 1:length(sortedAnchor);

%anchorCandidate = sortedAnchor(diff([-100;sortedAnchor])>TIMINGmARGIN);
anchorIndex = indexNumber(diff([-100;sortedAnchor])>TIMINGmARGIN);
anchorCandidate = sortedAnchor(anchorIndex);
mObject.anchorTimeLocation = anchorCandidate;
nFrequency = mObject.maximumFrequencyPoints;

nAnchor = length(anchorCandidate);
sortedAnchor(end+1) = sortedAnchor(end)+1;
anchorIndex(end+1) = anchorIndex(end)+1; % Terminator
frequencyAnchor = zeros(nAnchor,nFrequency);
for ii=1:nAnchor
    iFrequency = 0;
    anchorLocation = 0;
    for jj=1:min(nFrequency,anchorIndex(ii+1)-anchorIndex(ii)+1)
        if sortedAnchor((jj-1)+anchorIndex(ii)) < sortedAnchor(anchorIndex(ii+1))
            frequencyAnchor(ii,jj) = sortedFrequency((jj-1)+anchorIndex(ii));
            anchorLocation = anchorLocation+sortedAnchor((jj-1)+anchorIndex(ii));
            iFrequency = iFrequency+1;
        end;
    end;
    if iFrequency>1
        [dmy1,indsrt] = sort(frequencyAnchor(ii,1:iFrequency));
        frequencyAnchor(ii,1:iFrequency) = frequencyAnchor(ii,indsrt);
        mObject.anchorTimeLocation(ii) = anchorLocation/iFrequency;
    end;
end;
mObject.anchorFrequency = frequencyAnchor;

author: HidekiKawahara


text:
% MATCORR - Find matching rows in two matrices and their corrs.
%             Uses the Hungarian (default), VAM, or maxcorr assignment methods.
%             (Follow with MATPERM to permute and sign x -> y).
%
% Usage: >> [corr,indx,indy,corrs] = matcorr(x,y,rmmean,method,weighting);
%
% Inputs:
%   x     = first input matrix
%   y     = matrix with same number of columns as x
%
% Optional inputs:
%   rmmean    = When present and non-zero, remove row means prior to correlation
%               {default: 0}
%   method    = Method used to find assignments.
%               0= Hungarian Method - maximize sum of abs corrs {default: 2}
%               1= Vogel's Assignment Method -find pairs in order of max contrast
%               2= Max Abs Corr Method - find pairs in order of max abs corr
%               Note that the methods 0 and 1 require matrices to be square.
%   weighting = An optional weighting matrix size(weighting) = size(corrs) that
%               weights the corrs matrix before pair assignment {def: 0/[] -> ones}
% Outputs:
%   corr  = a column vector of correlation coefficients between
%           best-correlating rows of matrice x and y
%   indx  = a column vector containing the index of the maximum
%           abs-correlated x row in descending order of abs corr
%           (no duplications)
%   indy  = a column vector containing the index of the maximum
%           abs-correlated row of y in descending order of abs corr
%           (no duplications)
%   corrs = an optional square matrix of row-correlation coefficients
%           between matrices x and y
%
% Note: outputs are sorted by abs(corr)
%
% Authors: Scott Makeig & Sigurd Enghoff, SCCN/INC/UCSD, La Jolla, 11-30-96

% Copyright (C) 11-30-96 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 04-22-99 Re-written using VAM by Sigurd Enghoff, CNL/Salk
% 04-30-99 Added revision of algorithm loop by SE -sm
% 05-25-99 Added Hungarian method assignment by SE
% 06-15-99 Maximum correlation method reinstated by SE
% 08-02-99 Made order of outputs match help msg -sm
% 02-16-00 Fixed order of corr output under VAM added method explanations,
%          and returned corr signs in abs max method -sm
% 01-25-02 reformated help & license, added links -ad

% Uses function hungarian.m

function [corr,indx,indy,corrs] = matcorr(x,y,rmmean,method,weighting)
%
if nargin < 2 || nargin > 5
  help matcorr
  return
end

if nargin < 4
	method = 2; % default: Max Abs Corr - select successive best abs(corr) pairs
end

[m,n] = size(x);
[p,q] = size(y);
m = min(m,p);

if m~=n || p~=q
   if nargin>3 && method~=2
     fprintf('matcorr(): Matrices are not square: using max abs corr method (2).\n');
   end
   method = 2; % Can accept non-square matrices
end

if n~=q
  error('Rows in the two input matrices must be the same length.');
end

if nargin < 3 || isempty(rmmean)
  rmmean = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rmmean
  x = x - mean(x')'*ones(1,n); % optionally remove means
  y = y - mean(y')'*ones(1,n);
end
dx = sum(x'.^2);
dy = sum(y'.^2);
dx(find(dx==0)) = 1;
dy(find(dy==0)) = 1;
corrs = x*y'./sqrt(dx'*dy);

if nargin > 4 && ~isempty(weighting) && norm(weighting) > 0,
  if any(size(corrs) ~= size(weighting))
    fprintf('matcorr(): weighting matrix size must match that of corrs\n.')
    return
  else
	corrs = corrs.*weighting;
  end
end

cc = abs(corrs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch method
case 0
	ass = hungarian(-cc);   % Performs Hungarian algorithm matching

	idx1 = sub2ind(size(cc),ass,1:m);
	[dummy idx2] = sort(-cc(idx1));
	corr = corrs(idx1);
	corr = corr(idx2)';
	indy = [1:m]';
	indx = ass(idx2)';
	indy = indy(idx2);

case 1                      % Implements the VAM assignment method
	indx = zeros(m,1);
	indy = zeros(m,1);
	corr = zeros(m,1);

	for i=1:m,
		[sx ix] = sort(cc);  % Looks for maximum salience along a row/column
		[sy iy] = sort(cc'); % rather than maximum correlation.
		[sxx ixx] = max(sx(end,:)-sx(end-1,:));
		[syy iyy] = max(sy(end,:)-sy(end-1,:));

		if sxx == syy
			if sxx == 0 && syy == 0
        	   [sxx ixx] = max((sx(end,:)-sx(end-1,:)) .* sx(end,:));
    	       [syy iyy] = max((sy(end,:)-sy(end-1,:)) .* sy(end,:));
			else
				sxx = sx(end,ixx); % takes care of identical vectors
				syy = sy(end,iyy); % and zero vectors
			end
		end

		if sxx > syy
			indx(i) = ix(end,ixx);
			indy(i) = ixx;
		else
			indx(i) = iyy;
			indy(i) = iy(end,iyy);
		end
		cc(indx(i),:) = -1;
		cc(:,indy(i)) = -1;
	end

	i = sub2ind(size(corrs),indx,indy);
	corr = corrs(i);

	[tmp j] = sort(-abs(corr)); % re-sort by abs(correlation)
	corr = corr(j);
	indx = indx(j);
	indy = indy(j);

case 2                       % match successive max(abs(corr)) pairs
	indx = zeros(size(cc,1),1);
	indy = zeros(size(cc,1),1);
	corr = zeros(size(cc,1),1);

	for i = 1:size(cc,1)
		[tmp j] = max(cc(:));
		% [corr(i) j] = max(cc(:));
		[indx(i) indy(i)] = ind2sub(size(cc),j);
        corr(i) = corrs(indx(i),indy(i));
		cc(indx(i),:) = -1; % remove from contention
		cc(:,indy(i)) = -1;
	end

otherwise
	error('Unknown method');
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function angle = i4_to_angle ( i )

%*****************************************************************************80
%
%% I4_TO_ANGLE maps integers to points on a circle.
%
%  Discussion:
%
%    The angles are intended to be used to select colors on a color
%    hexagon whose 6 vertices are red, yellow, green, cyan, blue,
%    magenta.
%
%  Example:
%
%     I   X      ANGLE
%
%     0   0/3      0
%     1   1/3    120
%     2   2/3    240
%
%     3   1/6     60
%     4   3/6    180
%     5   5/6    300
%
%     6   1/12    30
%     7   3/12    90
%     8   5/12   150
%     9   7/12   210
%    10   9/12   270
%    11  11/12   330
%
%    12   1/24    15
%    13   3/24    45
%    14   5/24    75
%    etc
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the desired color.
%
%    Output, real ANGLE, an angle, measured in degrees,
%    between 0 and 360.
%
  if ( 0 <= abs ( i ) && abs ( i ) <= 2 )

    angle = 120.0 * abs ( i );

  else

    i1 = i4_log_2 ( floor ( abs ( i ) / 3 ) );
    i2 = abs ( i ) + 1 - 3 * 2^i1;
    i3 = 2 * ( i2 - 1 ) + 1;
    i4 = 3 * 2^( i1 + 1 );

    angle = 360.0 * i3 / i4;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

syms q1 q2

A01 = dh_sym(q1, L3, 0, pi/2);
A12 = dh_sym(q2, 0, 0, -pi/2);
A02 = A01*A12;

z0 = [0 0 1]';
z1 = A01(1:3,3);
z2 = A02(1:3,3);

Jw = [z0 z1 z2];

singularities = det(Jw)
singularities = simplify(singularities)
q=solve(singularities == 0)

function A = dh_sym(theta, d, a, alpha)
syms q1 q2
ca = cos(alpha);
sa = sin(alpha);
if abs(ca) < 1e-6
    ca = 0;
end
if abs(sa) < 1e-6
    sa = 0;
end

A=[cos(theta)  -ca*sin(theta)   sa*sin(theta)   a*cos(theta);
   sin(theta)   ca*cos(theta)  -sa*cos(theta)   a*sin(theta);
            0              sa             ca             d;
            0         0                     0              1];

repo: ARTE


text:
% CONVERT ROTATION MATRIX TO QUATERNION
function [q] = GetQuaternionFromRotationMatrix(R)
% This function is designed to convert from a rotation matrix
% to an equivalent quaternion. This function is also parallel
% to "rotm2quat.m".

q = zeros(4,1);
% Assemble the quaternion elements
q(1) = 0.5*sqrt(1 + R(1,1) + R(2,2) + R(3,3));
q(2) = (R(3,2) - R(2,3))/(4*q(1));
q(3) = (R(1,3) - R(3,1))/(4*q(1));
q(4) = (R(2,1) - R(1,2))/(4*q(1));

assert(any(isnan(q)) == 0,'quaternion calculation failed');

% Normalise the quaternion
q = unit(q);
end

repo: OpenMAS


text:
function [GrainBoundaryX,GrainBoundaryY]=IdentifyGrainBoundary(x,y,state)

count=1;
for i=1:size(x,1)-1
    for j=1:(size(x,1)-1)
        if state(i,j)~=state(i,j+1)
            GrainBoundaryX(1,count)=x(i,j);
            GrainBoundaryY(1,count)=y(i,j);
            count=count+1;
        end
        if state(i,j)~=state(i+1,j)
            GrainBoundaryX(1,count)=x(i,j);
            GrainBoundaryY(1,count)=y(i,j);
            count=count+1;
        end
    end
end

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34784-monte-carlo-simulation-of-two-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 2D square-lattice - microstructure/IdentifyGrainBoundary.m


text:
 function xs = de_pl_osps(x, Gb, ymi, Im, rmi, R, ftab, Gtab, masseff, niter, pixmax, gi, denom)
%function xs = de_pl_osps(x, Gb, ymi, Im, rmi, R, ftab, Gtab, masseff, niter, pixmax, gi, denom)
% The Dual-Energy PL-OS-SPS algorithm for transmission Poisson problem
% (ordered subsets, separable paraboloidal surrogates)
% uses fast precomputed curvatures described in paper fessler::
% in:
%	x		[np,2] initial guess
%	Gb		Gblock object
%	ymi		[nb,na,2] raw polyenergetic measurements
%	Im		[2] source intensities
%	rmi		[nb,na,2] (or scalar) background/scatter
%	R		penalty object (see Robject.m)
%	ftab		F table
%	Gtab		G table
%	masseff		[2,2] effective mass atten coef's
%			2 energies by 2 materials
%	niter		# iterations (including initial guess)
%			fix: allow relaxation parameters here!
%	pixmax		upper constraint for pixel values
%		can be scalar (e.g. 'inf') or an array the size of x
%	gi		sum(Gt)'
%		denom		[np,2] (if precomputed denominator)
% out:
%	x [np,2,niter]	updated image vectors each iteration
%
% Copyright 2001-04-28, Jeff Fessler, The University of Michigan

if nargin < 6, ir_usage, end

if ~isvar('Im') | isempty(Im)
	Im = [1 1];
end
if ~isvar('rmi') | isempty(rmi)
	rmi = zeros(size(ymi));
end

trl_check(ymi(:,:,1), [], rmi(:,:,1));
trl_check(ymi(:,:,2), [], rmi(:,:,2));

if ~isvar('R') | isempty(R), error 'R required', end

nblock = block_ob(Gb, 'n');
if ~isvar('niter')	| isempty(niter),	niter = 1;	end
if ~isvar('pixmax')	| isempty(pixmax),	pixmax = inf;	end
if ~isvar('chat')	| isempty(chat),	chat = 0;	end

starts = subset_start(nblock);
[nb, na, nm] = size(ymi);

%
% precompute denominator
%
if ~isvar('denom'), error 'use de_pl_denom', end

%
% loop over iterations
%
xs = zeros([size(x) niter]);	% [np,2,niter]
x = max(x,0);
xs(:,:,1) = x;
for it=2:niter
	printf('DE-PL-OSPS iteration %d', it)

	%
	% loop over subsets
	%
	for iset=1:nblock
		iblock = starts(iset);
		ia = iblock:nblock:na;

		% s=G*x "mass integrals"
		s1 = Gb{iblock} * x(:,1);
		s2 = Gb{iblock} * x(:,2);

		% predicted meas. means
		s1 = reshape(s1, nb, length(ia));
		s2 = reshape(s2, nb, length(ia));
		fh = ftab.feval(ftab, s1, s2);
		yb1 = Im(1) * exp(-fh(:,:,1)) + rmi(:,ia,1);
		yb2 = Im(2) * exp(-fh(:,:,2)) + rmi(:,ia,2);

		dh1 = ymi(:,ia,1) ./ yb1 - 1;	% m=1
		dh2 = ymi(:,ia,2) ./ yb2 - 1;	% m=2

	if 0
		G11 = interp2(Gtab.s1, Gtab.s2, Gtab.Gm1(:,:,1)', s1, s2);
		G21 = interp2(Gtab.s1, Gtab.s2, Gtab.Gm1(:,:,2)', s1, s2);
		G12 = interp2(Gtab.s1, Gtab.s2, Gtab.Gm2(:,:,1)', s1, s2);
		G22 = interp2(Gtab.s1, Gtab.s2, Gtab.Gm2(:,:,2)', s1, s2);

		% undo nonlinearity associated with table of Gml
		G11 = -exp(-G11);
		G12 = -exp(-G12);
		G21 = -exp(-G21);
		G22 = -exp(-G22);

	else
		% Gml using polynomial model!
		% note: 'coef' is indexed by m, but d_l
if it*iset == 2, warning 'not tested', end
		d11 = reshape(ftab.basis_d1(s1(:), s2(:)) * ftab.coef(:,1), ...
				size(s1));
		d21 = reshape(ftab.basis_d1(s1(:), s2(:)) * ftab.coef(:,2), ...
				size(s1));
		d12 = reshape(ftab.basis_d2(s1(:), s2(:)) * ftab.coef(:,1), ...
				size(s2));
		d22 = reshape(ftab.basis_d2(s1(:), s2(:)) * ftab.coef(:,2), ...
				size(s2));
		G11 = -(yb1-rmi(:,ia,1)) .* d11;
		G21 = -(yb2-rmi(:,ia,2)) .* d21;
		G12 = -(yb1-rmi(:,ia,1)) .* d12;
		G22 = -(yb2-rmi(:,ia,2)) .* d22;
	end

		% g_l is sum over m
		g1 = G11 .* dh1 + G21 .* dh2;
		g2 = G12 .* dh1 + G22 .* dh2;

		grad = Gb{iblock}' * [g1(:) g2(:)];	% [np,2]

%%		for isub=1:1
			for ll=1:2
				pgrad(:,ll) = R{ll}.cgrad(R{ll}, x(:,ll));
				Rdenom(:,ll) = R{ll}.denom(R{ll}, x(:,ll));
			end

			num = nblock * grad - pgrad;
			den = denom + Rdenom;	% semi-constant denom

			x = x + num ./ den;		% the update!
%%			x = x + ((x-xold) .* Lden + num) ./ den; % the update!
			x = max(x,0);		% enforce nonnegativity
			x = min(x,pixmax);	% enforce upper bound constraint
%%		end
	end

	if chat, printf('Range %g %g', min(x), max(x)), end
	xs(:,:,it) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [hidden, output] = ForwardPropogation(wConvol, wHidden, wBias, xTrain)
    % Dimension: (numImages, numFeatures)
    xTrain = xTrain(:, 1: 1024);
    numImages = size(xTrain,1);

    % wConvol is (filterDim, filterDim, numFilters) dimension
    numFilters = size(wConvol,3);
    filterDim = size(wConvol,2);

    imageDim = 32;
    convDim = imageDim - filterDim + 1;
    poolSize = 2;
    poolDim = convDim / poolSize;

    convolutions = zeros(convDim, convDim, numFilters, numImages);
    pooling = zeros(poolDim, poolDim, numFilters, numImages);
    images = reshape(xTrain, imageDim, imageDim,[]);

    for imageNum = 1: numImages
        for filterNum = 1: numFilters
            % Convolution Layer
            convolvedImage = zeros(convDim, convDim);
            filter = wConvol(:,:,filterNum);
            filter = rot90(squeeze(filter),2);
            im = squeeze(images(:,:,imageNum));
            for i = 1 : convDim
                for j = 1 : convDim
                    temp = double(im(i:i+filterDim-1,j:j+filterDim-1));
                    temp = temp .* filter;
                    convolvedImage(i,j) = sum(temp(:));
                    bias = wBias(filterNum);
                    convolvedImage(i,j) = sigmoid(convolvedImage(i,j) + bias);
                end
            end
            convolutions(:,:,filterNum, imageNum) = convolvedImage;

            % Pooling Layer
            pooledImage = zeros(poolDim, poolDim);
            for i = 1 : poolDim
                for j = 1 : poolDim
                    x = ((i-1) * poolSize) + 1;
                    y = ((j-1) * poolSize) + 1;
                    temp = convolvedImage(x:x+poolSize-1, y:y+poolSize-1);
                    pooledImage(i,j) = max(temp(:));
                end
            end
            pooling(:,:,filterNum, imageNum) = pooledImage;
        end
    end

    % Hidden Layer
    % Dimension is (hiddenSize, numImages)
    hidden = reshape(pooling,[],numImages);

    % Output Layer
    % Calculate the probability of each output unit/Label
    output = sigmoid(wHidden * hidden);
    output = output';
    hidden = hidden';

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% demo for image registration

close all
clear variables

% 1: create original 2D image
[X,Y] = meshgrid(-2:.2:2, -2:.2:2);
Z = X .* exp(-X.^2 - Y.^2);
inG.Y = Z;
[nx,ny] = size(inG.Y);

% 2: dis-align original image with rotation & translation
phi = [2;-2;-1];
sigma = 1e6;
[gx] = g_rigid2D([],phi,[],inG);
y = gx +  (1/sqrt(sigma))*randn(size(gx));

% 3: display original and disaligned images
hf = figure;
subplot(2,2,1),imagesc(inG.Y),colormap(bone),colorbar,axis equal, axis tight,title('original image')
subplot(2,2,2),imagesc(reshape(y,nx,ny)),colormap(bone),colorbar,axis equal, axis tight,title('disaligned image (target)')

% 4: undo dis-alignement:
%   - data consist of original image
%   - alignement parameters are not known
g_fname = @g_rigid2D;
dim.n = 0;
dim.n_theta = 0;
dim.n_phi = 3;
options.inG = inG;
[posterior,out] = VBA_NLStateSpaceModel(...
    y,...
    [],...
    [],...
    g_fname,...
    dim,...
    options);
hg = out.suffStat.gx;

% 5: display realigned and residuals images
figure(hf)
subplot(2,2,3),imagesc(reshape(hg,nx,ny)),colormap(bone),colorbar,axis equal, axis tight,title('realigned image')
subplot(2,2,4),imagesc(reshape(hg-y,nx,ny)),colormap(bone),colorbar,axis equal, axis tight,title('realigned - disaligned')

displayResults(posterior,out,y,[],[],[],phi,[],sigma);

author: MBB-team

repo: VBA-toolbox


text:
function [t]=vec2strIntDouble(varargin)

%%
switch nargin
    case 1
        n=varargin{1};
        formatDouble='%6.7e';
        rowWrapLength=[];
    case 2
        n=varargin{1};
        formatDouble=varargin{2};
        rowWrapLength=[];
    case 3
        n=varargin{1};
        formatDouble=varargin{2};
        rowWrapLength=varargin{3};
end

%%

optionStruct.formatDouble=formatDouble;
optionStruct.rowWrapLength=rowWrapLength;
[t]=mat2strIntDouble(n,optionStruct);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef LatitudeTermCondition < AbstractEventTerminationCondition
    %LatitudeTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lat(1,1) double = 0; %rad
        bodyInfo KSPTOT_BodyInfo
        direction(1,1) double = 0;
    end

    methods
        function obj = LatitudeTermCondition(lat)
            obj.lat = lat;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Latitude (%.3f deg)', rad2deg(obj.lat));
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = true;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Latitude';
            params.paramUnit = 'deg';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = obj.lat;
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = LatitudeOptimizationVariable(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = LatitudeTermCondition(paramValue);
        end
    end

    methods(Access=protected)
        function [value,isterminal,direction] = eventTermCond(obj, t,y)
            rVect = y(1:3);
            vVect = y(4:6);
            cartElem = CartesianElementSet(t, rVect(:), vVect(:), obj.bodyInfo.getBodyCenteredInertialFrame());
            geoElem = cartElem.convertToFrame(obj.frame).convertToGeographicElementSet();

            value = geoElem.lat - obj.lat;
            isterminal = 1;
            direction = obj.direction;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef (Abstract) mme_branch_opf_ac < mp.mme_branch_opf

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function obj = add_constraints(obj, mm, nm, dm, mpopt)
            %% find branches with flow limits
            dme = obj.data_model_element(dm);
            nme = obj.network_model_element(nm);
            ibr = find(dme.rate_a ~= 0 & dme.rate_a < 1e10);
            nl2 = length(ibr);      %% number of constrained branches
            mm.userdata.flow_constrained_branch_idx = ibr;

            if nl2
                %% port indexes
                nl = nme.nk;
                idx = [ibr; nl+ibr];

                %% limits
                flow_max = dme.rate_a(ibr); %% RATE_A

                %% branch flow constraints
                lim_type = upper(mpopt.opf.flow_lim(1));
                if lim_type == 'S'
                    fcn_flow = @(x)port_apparent_power_lim_fcn(nme, ...
                        mm.convert_x_m2n(x, nm), nm, idx, ...
                        [flow_max; flow_max] .^ 2);
                    hess_flow = @(x, lam)port_apparent_power_lim_hess(nme, ...
                        mm.convert_x_m2n(x, nm), lam, nm, idx);
                elseif lim_type == 'P'
                    fcn_flow = @(x)port_active_power_lim_fcn(nme, ...
                        mm.convert_x_m2n(x, nm), nm, idx, [flow_max; flow_max]);
                    hess_flow = @(x, lam)port_active_power_lim_hess(nme, ...
                        mm.convert_x_m2n(x, nm), lam, nm, idx);
                elseif lim_type == '2' || lim_type == 'P'
                    fcn_flow = @(x)port_active_power2_lim_fcn(nme, ...
                        mm.convert_x_m2n(x, nm), nm, idx, ...
                        [flow_max; flow_max] .^ 2);
                    hess_flow = @(x, lam)port_active_power2_lim_hess(nme, ...
                        mm.convert_x_m2n(x, nm), lam, nm, idx);
                elseif lim_type == 'I'
                    fcn_flow = @(x)port_current_lim_fcn(nme, ...
                        mm.convert_x_m2n(x, nm), nm, idx, ...
                        [flow_max; flow_max] .^ 2);
                    hess_flow = @(x, lam)port_current_lim_hess(nme, ...
                        mm.convert_x_m2n(x, nm), lam, nm, idx);
                else
                    error('mp.mme_branch_opf_ac/add_constraints: MPOPT.opf.flow_lim = ''%s'' not yet implemented.', mpopt.opf.flow_lim);
                end

                mm.add_nln_constraint({'Sf', 'St'}, [nl2;nl2], 0, fcn_flow, hess_flow);
            end
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            nme = obj.network_model_element(nm);

            %% branch active power flow
            pp = nm.get_idx('port');
            S_fr = nm.soln.gs_(pp.i1.branch(1):pp.iN.branch(1)) * dm.base_mva;
            S_to = nm.soln.gs_(pp.i1.branch(2):pp.iN.branch(2)) * dm.base_mva;

            %% shadow prices on branch flow constraints
            ibr = mm.userdata.flow_constrained_branch_idx;
            mu_flow_fr_ub = zeros(nme.nk, 1);
            mu_flow_to_ub = mu_flow_fr_ub;
            if length(ibr)
                lim_type = upper(mpopt.opf.flow_lim(1));
                nni = mm.get_idx('nli');
                lambda = mm.soln.lambda;
                if lim_type == 'P'
                    mu_flow_fr_ub(ibr) = lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf);
                    mu_flow_to_ub(ibr) = lambda.ineqnonlin(nni.i1.St:nni.iN.St);
                else
                    rate_a = dme.rate_a(ibr);
                    mu_flow_fr_ub(ibr) = 2 * lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf) .* rate_a;
                    mu_flow_to_ub(ibr) = 2 * lambda.ineqnonlin(nni.i1.St:nni.iN.St) .* rate_a;
                end
            end

            %% shadow prices on angle difference limits
            [mu_vad_lb, mu_vad_ub] = obj.ang_diff_prices(mm, nme);

            %% update in the data model
            dme.tab.pl_fr(dme.on) = real(S_fr);
            dme.tab.ql_fr(dme.on) = imag(S_fr);
            dme.tab.pl_to(dme.on) = real(S_to);
            dme.tab.ql_to(dme.on) = imag(S_to);
            dme.tab.mu_flow_fr_ub(dme.on) = mu_flow_fr_ub / dm.base_mva;
            dme.tab.mu_flow_to_ub(dme.on) = mu_flow_to_ub / dm.base_mva;
            dme.tab.mu_vad_lb(dme.on) = mu_vad_lb * pi/180;
            dme.tab.mu_vad_ub(dme.on) = mu_vad_ub * pi/180;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function X = imag(X)
% REAL (overloaded)

X.basis = imag(X.basis);
X.conicinfo = [0 0];
X = clean(X);


text:
function s4_resave_annot_after_temporal()
%image chnaged
% stb = load('annot_test.mat');
% base_annot = stb.videos;
% clear stb;clear all
close all
stb = load('annotV5.mat');
%stb = load('annot.mat');
base_annot = stb.videos;
clear stb;

stp = load('annot_full_Philippe.mat'); % train+test
phill_annot = stp.videos;
clear stp;

stp = load('annot_apt.mat'); % train+test
apt_annot = stp.videos;
clear stp;

baseList = {};
for i=1:length(base_annot)
    baseList{i,1} = base_annot(i).name; % store all the cor vid in a cellarray
end

philList = {};
for i=1:length(phill_annot)
    philList{i,1} = phill_annot(i).name; % store all the cor vid in a cellarray
end

aptList = {};
for i=1:length(apt_annot)
    aptList{i,1} = apt_annot(i).name; % store all the cor vid in a cellarray
end

images_dir = '/mnt/sun-alpha/datasets/UCF101/images';
save_path = './saved/';
%for v = 403 % : length(base_annot)
% rvid = randperm(length(base_annot));

% FigHandle = figure(1);
% FigHandle = figure('Position', [100, 100, 1049, 895]);
% hold on
count = 0;
load('tdata.mat')
% merged_annot = struct();
load('merged_annot');
load('checkpoint_improved.mat')
v = 1;vc =0;
while v <= length(baseList)

    videoName = baseList{v,1};

    base_tubes = base_annot(v).tubes;

    index = find(strcmp(philList,videoName));
    phil_tubes = phill_annot(index).tubes;

    index = find(strcmp(aptList,videoName));
    apt_tubes = apt_annot(index).tubes;

    num_imgs = base_annot(v).num_images;
    class = base_annot(v).tubes(1).class;

    merged_annot(v).num_imgs = num_imgs;
    merged_annot(v).name = base_annot(v).name;
    base_annot(v).name;

    if ~isgood(v) && strcmp(inputs{v},'t')
        vc = vc +1;
        newcell = tdata{v};
        tline = newcell{1};
        stl = strsplit(tline,' ');

        if length(stl{1})>1
            inputs{v} = 'ra';
            isgood(v) = 0;
        else
            merged_annot(v).tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class);
            inputs{v} = 'y';
            isgood(v) = 1;
        end

    end
    v=v+1;

end
% vc
fprintf('%d\n',sum(isgood));
save('checkpoint_improved_temporal.mat','isgood','inputs','v');
save('merged_annot_temporal.mat','merged_annot')

function tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class)
tubes = struct();
alltubes = {base_tubes,phil_tubes,apt_tubes};
for t = 1 : length(newcell)
    tline = newcell{t};
    stl = strsplit(tline,' ');

    if length(stl)<6
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,0);
    else
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,1);
    end
    tubes(t).sf = sf;
    tubes(t).ef = ef;
    tubes(t).class = class;
    tubes(t).boxes = boxes;
end

function [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,isone)
tube  = [];
annotType = {'b','p','a'};
ind = find(strcmp(annotType,stl{1}));

annot = alltubes{ind};

tid = str2num(stl{2});

sf = max(str2num(stl{3}),annot(tid).sf);

ef = min(str2num(stl{4}),annot(tid).ef);
st = sf-annot(tid).sf+1;
et = ef-sf+st; ct = 1;
boxes = annot(tid).boxes(st:et,:);
if isone
    ind = find(strcmp(annotType,stl{5}));
    annot = alltubes{ind};
    tid = str2num(stl{6});
    sfa = str2num(stl{7});
    if (sfa-ef)>1
        for kk = 1:(sfa-ef-1)
            boxes = [boxes;boxes(end,:)];
        end
    end
    ef = str2num(stl{8});
    st = sfa-annot(tid).sf+1;
    et = ef-sfa+st; ct = 1;
    boxes = [boxes;annot(tid).boxes(st:et,:)];
end

author: gurkirt


text:
function data_dir = GetNestedDataDir()
data_dir = 'C:\Janelia2015\datasets';
end

author: xiuyechen

repo: FishExplorer


text:
function value = r8vec_any_nonzero ( n, a )

%*****************************************************************************80
%
%% R8VEC_ANY_NONZERO: ( any A nonzero ) for R8VEC's.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, real A(N), the vector.
%
%    Output, logical R8VEC_ANY_NONZERO is TRUE if any entry is nonzero.
%
  value = any ( a(1:n) ~= 0.0 );

  return
end


text:
function laguerre_polynomial_test02 ( )

%*****************************************************************************80
%
%% LAGUERRE_POLYNOMIAL_TEST02 tests L_POLYNOMIAL_COEFFICIENTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2012
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLYNOMIAL_TEST02\n' );
  fprintf ( 1, '  L_POLYNOMIAL_COEFFICIENTS determines polynomial coefficients of L(n,x).\n' );

  c = l_polynomial_coefficients ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L(%d) = \n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( c(i+1,j+1) == 0.0 )

      elseif ( j == 0 )
        fprintf ( 1, '  %g\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %g * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %g * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end


text:
% feldkamp_example.m
% example of how to use feldkamp.m for cone-beam CT reconstruction
% Copyright 2004-8-28, Nicole Caparanis, Patty Laskowsky, Taka Masuda,
% and Jeff Fessler, The University of Michigan
% modified version by Ajay Paidi

if ~isvar('proj')
	down = 1/2;
	nv = 240*down;%240*down;
	nh = 256*down;%256*down;
	na = 224*down;%224*down;
	ds = 1024/nh;
	dt = ds;
	dis_src_det = 949;
	dis_iso_det = 408;
	dis_src_iso = dis_src_det - dis_iso_det;
	dis_foc_src = inf; % flat detector panel
	offset_det_h = 0.25; % quarter detector
	offset_det_v = 0.0;
	horiz = ([-(nh-1)/2:(nh-1)/2]' - offset_det_h) * ds;
	verti = ([-(nv-1)/2:(nv-1)/2]' - offset_det_v) * dt;
	printf('rmax=%g', dis_src_iso*sin(atan(max(abs(horiz)) / dis_src_det)))

ell = [ ...
		[0 0 -72 150 200 15 0 5];
		[0 0 -32 150 200 20 0 10];
        [0 0 12 150 200 15 0 10];
		[0 0 52 150 200 20 0 7];
        [0 0 82 150 200 15 0 10];
        ];

	nx = 240*down;%256*down;
	ny = 256*down;%240*down;
	nz = 112*down;%184*down;
	dx = 2/down; dy = dx; dz = dx;
	x = ellipsoids(nx, ny, nz, ell, dx, dy, dz);

	% cone-beam system geometry, generalized from fan-beam geometry.
	% see ASPIRE users guide under tech. reports on web page for details.
	args = arg_pair('system', NaN, 'nx', nx, 'ny', ny, 'nz', nz, ...
		'nv', nv, 'nh', nh, 'na', na, 'support', 'all', ...
		'orbit', 360, 'orbit_start', 0, ...
		'pixel_size', dx, 'ray_spacing', ds, 'strip_width', 0, ...
		'dis_src_det', dis_src_det, ...
		'dis_iso_det', dis_iso_det, ...
		'dis_foc_src', dis_foc_src, ...
		'offset_source', 0, ...
		'offset_det_h', offset_det_h, ...
		'offset_det_v', offset_det_v);

	pl=230;
	im(pl+1, x, 'x'), cbar
	im(pl+2, proj, 'proj'), cbar
	drawnow
prompt
end

% cone-beam reconstruction
mask = true([nx ny nz]);
recon = feldkamp(proj, 'ramp', mask, args);

% show results (off-center slices worse than central slice)
im(pl+4, recon, 'recon'), cbar
im(pl+5, recon - x, 'error'), cbar
ix = 1:nx; iy = ny/2; iz = nz/2;
subplot(236)
plot(ix, x(ix,iy,iz), '-', ix, recon(ix,iy,iz), '--')
axis([1 nx -1 16]), legend('true', 'recon', 2)
title 'MIDDLE SLICE'
iz=12;
subplot(233)
plot(ix, x(ix,iy,iz), '-', ix, recon(ix,iy,iz), '--')
axis([1 nx -1 16]), legend('true', 'recon')
title(sprintf('SLICE %d', iz))
%figure
%subplot(2,1,1)
%ix1 = nx/2;
%iy1 = ny/2;
%iy = 1:ny;
%ix = 1:nx;
%iz = nx/2;
%plot(iy, x(ix1, iy, iz), '-', iy, recon(ix1,iy,iz), '--');
%title('plot of x vs. y on middle slice')


text:
function w = dbnc(dsTrain, varargin)

	prtrace(mfilename);

	% Settings for the different training algorithms.
    if(~exist('nnet','dir'))
		error('Neural network toolbox not found')
    end

    cInput_types = {'gaussian', 'logsigm'};

    %argument definition
    p = inputParser;   % Create instance of inputParser class.
    p.addRequired('dsTrain', @(x)( isdataset(x)));

    p.addParamValue('input_type', 'logsigm', @(x)( any(strcmpi(x,cInput_types))) );
    p.addParamValue('num_hidden', 3, @(x)( isnumeric(x) && x > 0 ));
    p.addParamValue('size_hidden', [], @(x)( (all(isnumeric(x)) && all(x > 0)) || isempty(x) ) );
    p.addParamValue('Initialization_dataset', [], @(x)( isdataset(x) )  );
    p.addParamValue('Validation_dataset', [], @(x)( isdataset(x) )  );

	% Check arguments
    try
        p.parse( dsTrain, varargin{:} );
    catch MyError
        rethrow(MyError);
    end

    dsTrain = p.Results.dsTrain;
    dsInitialization = p.Results.Initialization_dataset;
    dsValidation = p.Results.Validation_dataset;
    num_hidden = p.Results.num_hidden;
    size_hidden = p.Results.size_hidden;
    input_type = p.Results.input_type;

    % Dont know why this variable uses a lot of bytes to store at disk.
    clear p

    %Check the datasets
    islabtype(dsTrain,'crisp');
    isvaldfile(dsTrain,1,2); 							% At least 1 object per class, 2 classes
    dsTrain = testdatasize(dsTrain);
    %force domain to be in [0-1]
    w_scale_domain = scalem(dsTrain, 'domain');
    dsTrain = dsTrain * w_scale_domain;

    if( isempty(dsValidation) )
        dsValidation = dsTrain;
    else
        dsValidation = testdatasize(dsValidation);
        iscomdset(dsTrain, dsValidation);   							% Check whether training and tuning set match
        dsValidation = dsValidation * w_scale_domain;
    end

    if( isempty(dsInitialization) )
        dsInitialization = dsTrain;
    else
        dsInitialization = dsInitialization * w_scale_domain;
    end

    [train_m, train_k, train_c] = getsize(dsTrain);
    [ train_labels, train_lablist] = getnlab(dsTrain);

    if( isempty(size_hidden) )
        % num_hidden layers each of train_k size by default
        size_hidden = repmat(train_k, num_hidden, 1);
    end

    %train the first RBM on data
    top = +dsInitialization;
    %train all other RBM's on top of each other
    for ii=1:num_hidden
        [model{ii} , top ] = rbmBB(top, size_hidden(ii));
    end
    clear top

    model{num_hidden}.Wc = 0.1*randn(train_c, size(model{num_hidden}.W,2) );
    model{num_hidden}.cc = 0.1*randn(1, train_c);
    model{num_hidden}.labels = 1:train_c;

    %fine-tuning of the DBN
    targets = zeros(train_m, train_c);
    for ii = 1:train_c
        targets(train_labels == ii, ii) = 1;
    end

    for epoch = 1:20

        %%%%%%%%%%%%%%% PERFORM CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        max_iter=3;

        if epoch<6  % First update top-level weights holding other weights fixed.

            top = +dsTrain;
            for ii=1:num_hidden
                top = rbmVtoH(model{ii}, top);
            end

            VV = colvec([ model{num_hidden}.Wc'; model{num_hidden}.cc ]);
            Dim = [ size_hidden(end); train_c];
            XX = minimize(VV,'CG_CLASSIFY_INIT',max_iter,Dim,top,targets);
            aux = reshape(XX, size_hidden(end)+1, train_c);
            model{num_hidden}.Wc = aux(1:end-1,:)';
            model{num_hidden}.cc = aux(end,:);

        else

            VV = [];
            for ii=1:num_hidden
                VV = [VV; colvec([model{ii}.W; model{ii}.b ]) ];
            end
            VV = [VV; colvec([model{ii}.Wc'; model{ii}.cc ])];

            Dim = [ train_k; colvec(size_hidden); train_c];

            XX = minimize(VV,'CG_CLASSIFY',max_iter,Dim, +dsTrain,targets);

            xxx = 0;
            for ii=1:num_hidden
                aux = reshape(XX(xxx+1:xxx+(Dim(ii)+1)*Dim(ii+1)),Dim(ii)+1,Dim(ii+1));
                model{ii}.W = aux(1:end-1,:);
                model{ii}.b = aux(end,:);
                xxx = xxx+(Dim(ii)+1)*Dim(ii+1);
            end
            aux = reshape(XX(xxx+1:xxx+(Dim(ii+1)+1)*Dim(ii+2)),Dim(ii+1)+1,Dim(ii+2));
            model{ii}.Wc = aux(1:end-1,:)';
            model{ii}.cc = aux(end,:);

        end

    end

    % Create mapping.

    w = w_scale_domain*mapping('dbnPredict', 'trained', model, train_lablist, train_k, train_c);
    w = setname(w, 'DBN classifier');
%     w = setcost(w,a);

return

author: marianux

repo: ecg-kit

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


text:
function a = dif2_llt ( n )

%*****************************************************************************80
%
%% DIF2_LLT returns the Cholesky factor of the DIF2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    a(i,i) = sqrt ( i + 1 ) / sqrt ( i );
  end

  for i = 2 : n
    a(i,i-1) = - sqrt ( i - 1 ) / sqrt ( i );
  end

  return
end


text:
function call_price=bs_european_call(S, K, r, sigma, time)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% European put option using Black-Scholes' formula
%
%
% Reference:
%
% John Hull, "Options, Futures and other Derivative Securities",
% Prentice-Hall, second edition, 1993.
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%   S:      spot price
%   K:      strike price
%   r:      interest rate
%   sigma:  volatility
%   time:   time to maturity
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

time_sqrt = sqrt(time);

d1 = (log(S/K)+r*time)/(sigma*time_sqrt)+0.5*sigma*time_sqrt;
d2 = d1-(sigma*time_sqrt);

call_price  = S*normcdf(d1)-K*exp(-r*time)*normcdf(d2);


text:
function [xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points(X,om,T,f,c,k)

%project_points.m
%
%[xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points(X,om,T,f,c,k)
%
%Projects a 3D structure onto the image plane.
%
%INPUT: X: 3D structure in the world coordinate frame (3xN matrix for N points)
%       (om,T): Rigid motion parameters between world coordinate frame and camera reference frame
%               om: rotation vector (3x1 vector); T: translation vector (3x1 vector)
%       f: camera focal length in units of horizontal and vertical pixel units (2x1 vector)
%       c: principal point location in pixel units (2x1 vector)
%       k: Distortion coefficients (radial and tangential) (4x1 vector)
%
%OUTPUT: xp: Projected pixel coordinates (2xN matrix for N points)
%        dxpdom: Derivative of xp with respect to om ((2N)x3 matrix)
%        dxpdT: Derivative of xp with respect to T ((2N)x3 matrix)
%        dxpdf: Derivative of xp with respect to f ((2N)x2 matrix if f is 2x1, or (2N)x1 matrix is f is a scalar)
%        dxpdc: Derivative of xp with respect to c ((2N)x2 matrix)
%        dxpdk: Derivative of xp with respect to k ((2N)x4 matrix)
%
%Definitions:
%Let P be a point in 3D of coordinates X in the world reference frame (stored in the matrix X)
%The coordinate vector of P in the camera reference frame is: Xc = R*X + T
%where R is the rotation matrix corresponding to the rotation vector om: R = rodrigues(om);
%call x, y and z the 3 coordinates of Xc: x = Xc(1); y = Xc(2); z = Xc(3);
%The pinehole projection coordinates of P is [a;b] where a=x/z and b=y/z.
%call r^2 = a^2 + b^2.
%The distorted point coordinates are: xd = [xx;yy] where:
%
%xx = a * (1 + kc(1)*r^2 + kc(2)*r^4)      +      2*kc(3)*a*b + kc(4)*(r^2 + 2*a^2);
%yy = b * (1 + kc(1)*r^2 + kc(2)*r^4)      +      kc(3)*(r^2 + 2*b^2) + 2*kc(4)*a*b;
%
%The left terms correspond to radial distortion, the right terms correspond to tangential distortion
%
%Fianlly, convertion into pixel coordinates: The final pixel coordinates vector xp=[xxp;yyp] where:
%
%xxp = f(1)*xx + c(1)
%yyp = f(2)*yy + c(2)
%
%
%NOTE: About 90 percent of the code takes care fo computing the Jacobian matrices
%
%
%Important function called within that program:
%
%rodrigues.m: Computes the rotation matrix corresponding to a rotation vector
%
%rigid_motion.m: Computes the rigid motion transformation of a given structure

if nargin < 6,
   k = zeros(4,1);
   if nargin < 5,
      c = zeros(2,1);
      if nargin < 4,
         f = ones(2,1);
         if nargin < 3,
            T = zeros(3,1);
            if nargin < 2,
               om = zeros(3,1);
               if nargin < 1,
                  error('Need at least a 3D structure to project (in project_points.m)');
                  return;
               end;
            end;
         end;
      end;
   end;
end;

[m,n] = size(X);

[Y,dYdom,dYdT] = rigid_motion(X,om,T);

inv_Z = 1./Y(3,:);

x = (Y(1:2,:) .* (ones(2,1) * inv_Z)) ;

bb = (-x(1,:) .* inv_Z)'*ones(1,3);
cc = (-x(2,:) .* inv_Z)'*ones(1,3);

dxdom = zeros(2*n,3);
dxdom(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(1:3:end,:) + bb .* dYdom(3:3:end,:);
dxdom(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(2:3:end,:) + cc .* dYdom(3:3:end,:);

dxdT = zeros(2*n,3);
dxdT(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(1:3:end,:) + bb .* dYdT(3:3:end,:);
dxdT(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(2:3:end,:) + cc .* dYdT(3:3:end,:);

% Add distortion:

r2 = x(1,:).^2 + x(2,:).^2;

dr2dom = 2*((x(1,:)')*ones(1,3)) .* dxdom(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdom(2:2:end,:);
dr2dT = 2*((x(1,:)')*ones(1,3)) .* dxdT(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdT(2:2:end,:);

r4 = r2.^2;

dr4dom = 2*((r2')*ones(1,3)) .* dr2dom;
dr4dT = 2*((r2')*ones(1,3)) .* dr2dT;

% Radial distortion:

cdist = 1 + k(1) * r2 + k(2) * r4;

dcdistdom = k(1) * dr2dom + k(2) * dr4dom;
dcdistdT = k(1) * dr2dT+ k(2) * dr4dT;
dcdistdk = [ r2' r4' zeros(n,2)];

xd1 = x .* (ones(2,1)*cdist);

dxd1dom = zeros(2*n,3);
dxd1dom(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdom;
dxd1dom(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdom;
coeff = (reshape([cdist;cdist],2*n,1)*ones(1,3));
dxd1dom = dxd1dom + coeff.* dxdom;

dxd1dT = zeros(2*n,3);
dxd1dT(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdT;
dxd1dT(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdT;
dxd1dT = dxd1dT + coeff.* dxdT;

dxd1dk = zeros(2*n,4);
dxd1dk(1:2:end,:) = (x(1,:)'*ones(1,4)) .* dcdistdk;
dxd1dk(2:2:end,:) = (x(2,:)'*ones(1,4)) .* dcdistdk;

% tangential distortion:

a1 = 2.*x(1,:).*x(2,:);
a2 = r2 + 2*x(1,:).^2;
a3 = r2 + 2*x(2,:).^2;

delta_x = [k(3)*a1 + k(4)*a2 ;
   k(3) * a3 + k(4)*a1];

ddelta_xdx = zeros(2*n,2*n);
aa = (2*k(3)*x(2,:)+6*k(4)*x(1,:))'*ones(1,3);
bb = (2*k(3)*x(1,:)+2*k(4)*x(2,:))'*ones(1,3);
cc = (6*k(3)*x(2,:)+2*k(4)*x(1,:))'*ones(1,3);

ddelta_xdom = zeros(2*n,3);
ddelta_xdom(1:2:end,:) = aa .* dxdom(1:2:end,:) + bb .* dxdom(2:2:end,:);
ddelta_xdom(2:2:end,:) = bb .* dxdom(1:2:end,:) + cc .* dxdom(2:2:end,:);

ddelta_xdT = zeros(2*n,3);
ddelta_xdT(1:2:end,:) = aa .* dxdT(1:2:end,:) + bb .* dxdT(2:2:end,:);
ddelta_xdT(2:2:end,:) = bb .* dxdT(1:2:end,:) + cc .* dxdT(2:2:end,:);

ddelta_xdk = zeros(2*n,4);
ddelta_xdk(1:2:end,3) = a1';
ddelta_xdk(1:2:end,4) = a2';
ddelta_xdk(2:2:end,3) = a3';
ddelta_xdk(2:2:end,4) = a1';

xd2 = xd1 + delta_x;

dxd2dom = dxd1dom + ddelta_xdom ;
dxd2dT = dxd1dT + ddelta_xdT;
dxd2dk = dxd1dk + ddelta_xdk ;

% Pixel coordinates:
if length(f)>1,
    xp = xd2 .* (f * ones(1,n))  +  c*ones(1,n);
    coeff = reshape(f*ones(1,n),2*n,1);
    dxpdom = (coeff*ones(1,3)) .* dxd2dom;
    dxpdT = (coeff*ones(1,3)) .* dxd2dT;
    dxpdk = (coeff*ones(1,5)) .* dxd2dk;
    dxpdf = zeros(2*n,2);
    dxpdf(1:2:end,1) = xd2(1,:)';
    dxpdf(2:2:end,2) = xd2(2,:)';
else
    xp = f * xd2 + c*ones(1,n);
    dxpdom = f  * dxd2dom;
    dxpdT = f * dxd2dT;
    dxpdk = f  * dxd2dk;
    dxpdf = xd2(:);
end;

dxpdc = zeros(2*n,2);
dxpdc(1:2:end,1) = ones(n,1);
dxpdc(2:2:end,2) = ones(n,1);

return;

% Test of the Jacobians:

n = 10;

X = 10*randn(3,n);
om = randn(3,1);
T = [10*randn(2,1);40];
f = 1000*rand(2,1);
c = 1000*randn(2,1);
k = 0.5*randn(4,1);

[x,dxdom,dxdT,dxdf,dxdc,dxdk] = project_points(X,om,T,f,c,k);

% Test on om: OK

dom = 0.000000001 * norm(om)*randn(3,1);
om2 = om + dom;

[x2] = project_points(X,om2,T,f,c,k);

x_pred = x + reshape(dxdom * dom,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on T: OK!!

dT = 0.0001 * norm(T)*randn(3,1);
T2 = T + dT;

[x2] = project_points(X,om,T2,f,c,k);

x_pred = x + reshape(dxdT * dT,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on f: OK!!

df = 0.001 * norm(f)*randn(2,1);
f2 = f + df;

[x2] = project_points(X,om,T,f2,c,k);

x_pred = x + reshape(dxdf * df,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on c: OK!!

dc = 0.01 * norm(c)*randn(2,1);
c2 = c + dc;

[x2] = project_points(X,om,T,f,c2,k);

x_pred = x + reshape(dxdc * dc,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on k: OK!!

dk = 0.001 * norm(4)*randn(4,1);
k2 = k + dk;

[x2] = project_points(X,om,T,f,c,k2);

x_pred = x + reshape(dxdk * dk,2,n);

norm(x2-x)/norm(x2 - x_pred)

repo: ekfmonoslam