text:
function varargout = ProjectSetup(varargin)
% PROJECTSETUP MATLAB code for ProjectSetup.fig
%      PROJECTSETUP, by itself, creates a new PROJECTSETUP or raises the existing
%      singleton*.
%
%      H = PROJECTSETUP returns the handle to a new PROJECTSETUP or the handle to
%      the existing singleton*.
%
%      PROJECTSETUP('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PROJECTSETUP.M with the given input arguments.
%
%      PROJECTSETUP('Property','Value',...) creates a new PROJECTSETUP or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ProjectSetup_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ProjectSetup_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ProjectSetup

% Last Modified by GUIDE v2.5 11-Dec-2020 11:13:31

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ProjectSetup_OpeningFcn, ...
                   'gui_OutputFcn',  @ProjectSetup_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% -------------------------------------------------------------------------
% --- Executes just before ProjectSetup is made visible.
function ProjectSetup_OpeningFcn(hObject, ~, handles, varargin)
% ProjectSetup(p1,v1,...)
%
% ProjectSetup configures a JAABA project. It can optionally accept an
% existing project configuration (basicParamsStruct). It returns a
% configured project in the supplied 'handleobj'. A project configuration
% currently looks like a stripped-down Macguffin and consists of
%   - Names of behaviors/classifiers
%   - Name of experimental movie files
%   - etc.
%
% Call ProjectSetup as follows:
% hobj = HandleObj;
% uiwait(ProjectSetup('handleobj',hobj,...<other args>...));
% processResultsOrDoWhatever(hobj.data); % hobj.data is a Macguffin
%
% P-V arguments:
% * basicParamsStruct: Either [], or Macguffin (or Macguffin-like struct).
% Supplies initial values for project parameters. If basicParamStruct==[],
% then we are creating a new project.
% * handleobj: A HandleObj.

handles.output = hObject;

[figureJLabel, ...
 basicParamsStruct,...
 defaultmoviefilename,...
 defaultmovieindexfilename,...
 defaulttrxfilename,...
 defaultlbl_file,...
 handles.handleobj] = ...
   myparse(varargin,...
           'figureJLabel',[],...
           'basicParamsStruct',[],...
           'defaultmoviefilename',0,...
           'defaultmovieindexfilename',0,...
           'defaulttrxfilename',0,...
           'defaultlbl_file',0,...
           'handleobj',[]);
handles.new = isempty(basicParamsStruct);
handles.figureJLabel = figureJLabel;
handles.basicParamsStruct = basicParamsStruct;

adjustColorsIfMac(hObject);

% Need to derive the basic and advanced sizes (part of the model) from
% the current figure dimensions
handles = setBasicAndAdvancedSizesToMatchFigure(handles);
handles.mode = 'basic';
handles = updateFigurePosition(handles);
centerOnParentFigure(hObject,figureJLabel);
% ADD THIS FOR DEBUGGING WITH AUTOMATION
% set(handles.figureProjectSetup,'Position',[413   588   585   347]);

% Initialize the list of possible feature lexicon names
[handles.featureLexiconNameList,xmlList] = getFeatureLexiconListsFromXML();
handles.featureLexiconNameListIsST = cellfun(@xmlParamsIsST,xmlList);
assert(numel(handles.featureLexiconNameList)==numel(handles.featureLexiconNameListIsST));

% Populate the featureLexiconName popuplist with options
set(handles.featureconfigpopup,'String',handles.featureLexiconNameList);

if isempty(basicParamsStruct)
  % New project
  warnst = warning('off','MATLAB:structOnObject');
  handles.basicParamsStruct = struct(Macguffin('flies')); % the default featureLexiconName
  warning(warnst);
  if ischar(defaultmoviefilename)
    handles.basicParamsStruct.file.moviefilename = defaultmoviefilename;
  end
  if ischar(defaultmovieindexfilename)
    handles.basicParamsStruct.file.movieindexfilename = defaultmovieindexfilename;
  end
  if ischar(defaulttrxfilename)
    handles.basicParamsStruct.file.trxfilename = defaulttrxfilename;
  end
  if ischar(defaultlbl_file),
    handles.basicParamsStruct.lbl_file = defaultlbl_file;
  end
else
  warnst = warning('off','MATLAB:structOnObject');
  handles.basicParamsStruct = struct(basicParamsStruct);
  warning(warnst);

  % Quirk: basicParamsStruct.behaviors
  % For existing projects basicParamsStruct.behaviors.names includes the
  % no-behavior names (eg 'None', 'No_Chase', etc). In ProjectSetup, the
  % user currently only specifies the positive names.
  % So, within ProjectSetup, we truncate .behaviors.names to just the real
  % names; as a quirk, we leave .behaviors.labelcolors unmodified, so
  % .labelcolors DOES include no-behaviors, as the user may want to change
  % those colors.
  behnames = Labels.verifyBehaviorNames(handles.basicParamsStruct.behaviors.names);
  handles.basicParamsStruct.behaviors.names = behnames;
end

if ~isfield(handles.basicParamsStruct,'extra'),
  handles.basicParamsStruct.extra = struct;
end
if ~isfield(handles.basicParamsStruct.extra,'usePastOnly'),
  handles.basicParamsStruct.extra.usePastOnly = false;
end

% Update the configuration table in the GUI
updateConfigTable(handles);

set(hObject,'name',fif(handles.new,'New...','Edit...'));
onoff = fif(handles.new,'on','off');
set(handles.featureconfigpopup,'enable',onoff);
set(handles.pushbutton_copy,'visible',onoff);
set(handles.editName,'enable',onoff);
set(handles.edittrxfilename,'enable',onoff);
set(handles.pushbutton_perframe,'visible',onoff);

% Set the current score-as-feature file
fileNameList = {handles.basicParamsStruct.scoreFeatures(:).classifierfile};
handles.indexOfScoreFeaturesFile = fif(isempty(fileNameList),[],1);

guidata(hObject, handles);

% Update all the text fields in the figure
updateEditsListboxesAndPopupmenus(handles);

% Make the figure visible
set(hObject,'Visible','on');

return

% -------------------------------------------------------------------------
% --- Outputs from this function are returned to the command line.
function varargout = ProjectSetup_OutputFcn(hObject, ~, handles)  %#ok
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
%varargout{2} = handles.outfile;
%delete(handles.figureProjectSetup);
return

% -------------------------------------------------------------------------
function handles=setBasicAndAdvancedSizesToMatchFigure(handles)
% Derive the basic and advanced sizes (part of the model) from
% the current figure dimensions
curPos = get(handles.figureProjectSetup,'Position');
tablePos = get(handles.config_table,'Position');
reducedWidth = tablePos(1)-15;
reducedHeight = curPos(4);
handles.advancedSize = curPos(3:4);
handles.basicSize = [reducedWidth reducedHeight];
return

% -------------------------------------------------------------------------
function handles = updateFigurePosition(handles)
% Update the figure to match the current mode (which is in the model)
if strcmp(handles.mode,'basic')
  curPos = get(handles.figureProjectSetup,'Position');
  set(handles.figureProjectSetup,'Position',[curPos(1:2) handles.basicSize]);
  set(handles.togglebutton_advanced,'Value',0);
else
  curPos = get(handles.figureProjectSetup,'Position');
  set(handles.figureProjectSetup,'Position',[curPos(1:2) handles.advancedSize]);
  set(handles.togglebutton_advanced,'Value',1);
end
return

% -------------------------------------------------------------------------
function updateEditsListboxesAndPopupmenus(handles)

% update the behavior name editbox
if isfield(handles.basicParamsStruct.behaviors,'names');
  names = handles.basicParamsStruct.behaviors.names;
  if numel(names)>0
    namestr = String.cellstr2CommaSepList(names);
  else
    namestr = '';
  end
  set(handles.editName,'String',namestr);
else
  set(handles.editName,'String','');
end

fnames = {'scorefilename','moviefilename','trxfilename'};
boxnames = {'editscorefilename','editmoviefilename','edittrxfilename'};
for ndx = 1:numel(fnames)
  curf = fnames{ndx};
  curbox = boxnames{ndx};
  if isfield(handles.basicParamsStruct.file,curf);
    val = handles.basicParamsStruct.file.(curf);
    if iscellstr(val)
      val = String.cellstr2CommaSepList(val);
    end
    set(handles.(curbox),'String',val);
  else
    set(handles.(curbox),'String','');
  end
end

% Update the select feature dictionary name
indexOfFeatureLexicon = find(strcmp(handles.basicParamsStruct.featureLexiconName,handles.featureLexiconNameList));
if isempty(indexOfFeatureLexicon),
  indexOfFeatureLexicon = 1;
end
set(handles.featureconfigpopup,'Value',indexOfFeatureLexicon);
tfST = handles.featureLexiconNameListIsST(indexOfFeatureLexicon);
if tfST
  visVal = 'on';
else
  visVal = 'off';
end
set(handles.txST,'Visible',visVal);

% % Update the list of scores-an-inputs
% fileNameList = {handles.basicParamsStruct.scoreFeatures(:).classifierfile};
% set(handles.listbox_inputscores,'String',fileNameList);
% set(handles.listbox_inputscores,'Value',handles.indexOfScoreFeaturesFile);

% % Disble the Remove button iff the list of scores-as-inputs is empty
% set(handles.pushbutton_removelist,'enable',offIff(isempty(fileNameList)));

return

% -------------------------------------------------------------------------
function handles = fileNameEditTwiddled(handles,editName)
% Deal with one of the file name edits being modified.
% This is like a controller method.
varName=editName(5:end);  % trim 'edit' off of start
newFileName=strtrim(get(handles.(editName),'String'));
allowEmpty = ismember(varName,{'moviefilename'});
allowMulti = ismember(varName,{'scorefilename','moviefilename'});
[handles,whatHappened]=setFileName(handles,varName,newFileName,allowEmpty,allowMulti);
warnIfInvalid(varName,whatHappened);
updateFileNameEdit(handles,editName);
updateConfigTable(handles);
return

% -------------------------------------------------------------------------
function [handles,whatHappened] = setFileName(handles,varName,newFileName,allowEmpty,allowMulti)
% Set the given file name in the model, if it's a valid name.  whatHappened
% is a string that can be 'notChanged', 'emptyEntry','changed', or
% 'invalidEntry', depending.

fileName=handles.basicParamsStruct.file.(varName);

if allowEmpty && isempty(newFileName),
  handles.basicParamsStruct.file.(varName)=newFileName;
  if strcmp(newFileName,fileName),
    whatHappened='notChanged';
  else
    whatHappened='changed';
  end
  return;
end

if isequal(newFileName,fileName)
  whatHappened='notChanged';
elseif isempty(newFileName)
  whatHappened='emptyEntry';
elseif ~allowMulti && IsNiceFileName(newFileName)
  whatHappened='changed';
  handles.basicParamsStruct.file.(varName)=newFileName;
elseif allowMulti
  newFileNameCellStr = String.commaSepList2CellStr(newFileName);
  if all(cellfun(@IsNiceFileName,newFileNameCellStr))
    whatHappened = 'changed';
    handles.basicParamsStruct.file.(varName)=newFileNameCellStr;
  else
    whatHappened = 'invalidEntry';
  end
else
  whatHappened='invalidEntry';
end

return

% -------------------------------------------------------------------------
function warnIfInvalid(varName,whatHappened)
% Throw up a warning dialog if whatHappened equals 'invalidEntry'
if isequal(whatHappened,'invalidEntry')
  message=...
    sprintf(['The name specified for %s cannot have special characters.'...
             'Please use only alphanumeric characters and underscore.'],varName);
  title='Invalid File Name';
  uiwait(warndlg(message,title,'modal'));
end
return

% -------------------------------------------------------------------------
function updateFileNameEdit(handles,editName)
  % Update the named file name edit to match the model.
  varName=editName(5:end);  % trim 'edit' off of start
  val = handles.basicParamsStruct.file.(varName);
  if iscellstr(val)
    val = String.cellstr2CommaSepList(val);
  end
  set(handles.(editName),'String',val);
return

% -------------------------------------------------------------------------
function updateConfigTable(handles)
% Update the config table (a GUI element) to match the current "model"
% state
basicParamsStruct = handles.basicParamsStruct;
fields2remove = {'featureLexicon','windowFeaturesParams','scoreFeatures', ...
                 'sublexiconPFNames','labels','gtLabels','expDirNames', ...
                 'gtExpDirNames', 'classifierStuff', 'version','classifierfile',...
                 'windowFeaturesParams','aptInfo'};
if ~handles.new
  editModeRemove = {'featureLexiconName' 'behaviors.names' 'behaviors.type' ...
    'behaviors.nbeh' 'file.trxfilename' 'extra.perframe'};
  fields2remove = [fields2remove editModeRemove];
end
basicParamsStruct = structrmfield(basicParamsStruct,fields2remove);
data = GetParamsAsTable(basicParamsStruct);
set(handles.config_table,'Data',data);
return

% -------------------------------------------------------------------------
function data = GetParamsAsTable(basicParamsStruct)
data = addToList(basicParamsStruct,{},'');
idx = cellfun(@iscell,data(:,2));
if any(idx),
  for i = find(idx(:)'),
    if all(cellfun(@ischar,data{i,2})),
      data{i,2} = sprintf('%s,',data{i,2}{:});
      if numel(data{i,2})>0 && data{i,2}(end) == ',',
        data{i,2} = data{i,2}(1:end-1);
      end
    end
  end
end

if any(cellfun(@iscell,data(:,2))),
  data = {};
  return;
end

% -------------------------------------------------------------------------
function list = addToList(curStruct,list,pathTillNow)
if isempty(fieldnames(curStruct)), return; end
fnames = fieldnames(curStruct);
for ndx = 1:numel(fnames)
  if isstruct(curStruct.(fnames{ndx})),
    list = addToList(curStruct.(fnames{ndx}),list,[pathTillNow fnames{ndx} '.']);
  else
    list{end+1,1} = [pathTillNow fnames{ndx}]; %#ok<AGROW>
    param = curStruct.(fnames{ndx});
    if isnumeric(param)
      if isempty(param)
        q='';
      else
        q = num2str(param(1));
        for jj = 2:numel(param)
          q = [q ',' num2str(param(jj))]; %#ok<AGROW>
        end
      end
      list{end,2} = q;
    else
      list{end,2} = param;
    end
  end
end

% -------------------------------------------------------------------------
function editName_Callback(hObject, eventdata, handles)  %#ok
% hObject    handle to editName (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editName as text
%        str2double(get(hObject,'String')) returns contents of editName as a double
name = String.commaSepList2CellStr(get(hObject,'String'));
if any(cellfun(@(x)isempty(regexp(x,'^[a-zA-Z][\w_,]*$','once','start')),name));
   uiwait(warndlg(['The behavior name cannot have special characters.'...
                   'Please use only alphanumeric characters and _']));
   updateEditsListboxesAndPopupmenus(handles);
   return
end

handles.basicParamsStruct.behaviors.names = name;
if numel(name)>1
  % multiple classifiers; set labelcolors
  nlabels = 2*numel(name);
  clrs = Labels.cropOrAugmentLabelColors(zeros(1,0),nlabels,'lines');
  handles.basicParamsStruct.behaviors.labelcolors = clrs;
end
handles.basicParamsStruct.file.scorefilename = cellfun(@(x)sprintf('scores_%s.mat',x),name,'uni',0);
guidata(hObject,handles);
updateEditsListboxesAndPopupmenus(handles);
updateConfigTable(handles);
return

function behaviorsStruct = enforceBehaviorParamConstraints(behaviorsStruct,dowarn)

if ~exist('dowarn','var')
  dowarn = false;
end

% Add no-behavior names
behaviorsStruct.names = Labels.behnames2labelnames(behaviorsStruct.names);

nlabels = numel(behaviorsStruct.names);
if numel(behaviorsStruct.labelcolors)~=3*nlabels
  if dowarn
    warning('ProjectSetup:behaviorParams',...
      'Specified label colors not consistent with number of behaviors. Updating.');
  end
  behaviorsStruct.labelcolors = Labels.cropOrAugmentLabelColors(behaviorsStruct.labelcolors,nlabels,'lines');
end

% -------------------------------------------------------------------------
function editName_CreateFcn(hObject, eventdata, handles)  %#ok
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
return

% -------------------------------------------------------------------------
% --- Executes on selection change in featureconfigpopup.
function featureconfigpopup_Callback(hObject, eventdata, handles)  %#ok
% hObject    handle to featureconfigpopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
iFeatureLexicon=get(hObject,'Value');
featureLexiconNameNew = handles.featureLexiconNameList{iFeatureLexicon};
% Save the current values of the main four fields
% (Note that any changes in the advanced part of the window will be lost.)
behaviorName=handles.basicParamsStruct.behaviors.names;
movieFileName=handles.basicParamsStruct.file.moviefilename;
movieIndexFileName=handles.basicParamsStruct.file.movieindexfilename;
trackFileName=handles.basicParamsStruct.file.trxfilename;
scoreFileName=handles.basicParamsStruct.file.scorefilename;
% Replace basicParamsStruct with one appropriate to the new feature lexicon
% name
old=warning('query','MATLAB:structOnObject');
warning('off','MATLAB:structOnObject');  % turn off annoying warning
handles.basicParamsStruct = struct(Macguffin(featureLexiconNameNew));
warning(old);  % restore annoying warning
% Now restore fields from the old basicParamsStruct
handles.basicParamsStruct.behaviors.names=behaviorName;
handles.basicParamsStruct.file.moviefilename=movieFileName;
handles.basicParamsStruct.file.movieindexfilename=movieIndexFileName;
handles.basicParamsStruct.file.trxfilename=trackFileName;
handles.basicParamsStruct.file.scorefilename=scoreFileName;
guidata(hObject,handles);
updateConfigTable(handles);
updateEditsListboxesAndPopupmenus(handles);
return

% -------------------------------------------------------------------------
% --- Executes during object creation, after setting all properties.
function featureconfigpopup_CreateFcn(hObject, ~, ~)  %#ok
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -------------------------------------------------------------------------
function editscorefilename_Callback(hObject, eventdata, handles)  %#ok
editName=get(hObject,'tag');
handles = fileNameEditTwiddled(handles,editName);
guidata(hObject,handles);

% -------------------------------------------------------------------------
function editscorefilename_CreateFcn(hObject, eventdata, handles)  %#ok
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -------------------------------------------------------------------------
function editmoviefilename_Callback(hObject, eventdata, handles)  %#ok
editName=get(hObject,'tag');
handles = fileNameEditTwiddled(handles,editName);
guidata(hObject,handles);

% -------------------------------------------------------------------------
function editmoviefilename_CreateFcn(hObject, eventdata, handles)  %#ok
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -------------------------------------------------------------------------
function edittrxfilename_Callback(hObject, eventdata, handles)  %#ok
editName=get(hObject,'tag');
handles = fileNameEditTwiddled(handles,editName);
guidata(hObject,handles);

% -------------------------------------------------------------------------
% --- Executes during object creation, after setting all properties.
function edittrxfilename_CreateFcn(hObject, eventdata, handles)  %#ok
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% -------------------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)  %#ok
handles.outfile = 0;
guidata(hObject,handles);
delete(findAncestorFigure(hObject));
%uiresume(handles.figureProjectSetup);
return

% -------------------------------------------------------------------------
function pushbutton_done_Callback(hObject, ~, handles)  %#ok
behaviorName = handles.basicParamsStruct.behaviors.names;
if isempty(behaviorName)
  uiwait(errordlg('You must enter a behavior name.', ...
                  'No behavior name', ...
                  'modal'));
  return
end

% Check for an empty track file name
trackFileName = handles.basicParamsStruct.file.trxfilename;
if isempty(trackFileName)
  uiwait(errordlg('You must enter a track file name.', ...
                  'No track file name', ...
                  'modal'));
  return
end

% Check for an empty score file name
scoreFileName = handles.basicParamsStruct.file.scorefilename;
if isempty(scoreFileName)
  uiwait(errordlg('You must enter a score file name.', ...
                  'No score file name', ...
                  'modal'));
  return
end

nbehavior = numel(handles.basicParamsStruct.behaviors.names);
assert(iscellstr(handles.basicParamsStruct.file.scorefilename));
nscorefile = numel(handles.basicParamsStruct.file.scorefilename);
if nbehavior~=nscorefile
  uiwait(errordlg('The number of score files must match the number of behaviors.', ...
                  'Invalid score file names', ...
                  'modal'));
  return
end

% if movie is an mjpg, ask for an index file
if ischar(handles.basicParamsStruct.file.moviefilename)
  [~,~,ext] = fileparts(handles.basicParamsStruct.file.moviefilename);
else
  [~,~,ext] = fileparts(handles.basicParamsStruct.file.moviefilename{1});
end
if ismember(lower(ext),{'.mjpg','.mjpeg'}),

  res = questdlg('Is this an indexed MJPG file? If so, you will be prompted to select the corresponding index file location.');
  if strcmpi(res,'Yes'),
    if isfield(handles.basicParamsStruct.file,'movieindexfilename') && ...
        ischar(handles.basicParamsStruct.file.movieindexfilename),
      defaultindexfile = handles.basicParamsStruct.file.movieindexfilename;
    else
      defaultindexfile = 'index.txt';
    end
    res = inputdlg({'MJPG index file name: '},'MJPG index file name',1,{defaultindexfile});
    % cancel
    if isempty(res),
      return;
    end
    handles.basicParamsStruct.file.movieindexfilename = res{1};
  else
    handles.basicParamsStruct.file.movieindexfilename = 0;
  end

else
  handles.basicParamsStruct.file.movieindexfilename = 0;
end

% Finish initialization of basicParamsStruct and return as Macguffin in handleobj
handles.basicParamsStruct.behaviors = ...
  enforceBehaviorParamConstraints(handles.basicParamsStruct.behaviors,true);

if handles.new
  for i = nbehavior:-1:1
    cs(i,1) = ClassifierStuff();
    wfp{i,1} = handles.basicParamsStruct.windowFeaturesParams{1};
  end
  handles.basicParamsStruct.classifierStuff = cs;
  handles.basicParamsStruct.windowFeaturesParams = wfp;
end

basicParamsStruct = handles.basicParamsStruct;
% AL: Not sure why we return a Macguffin vs just the struct, prob doesn't
% make a difference
if basicParamsStruct.fromAPT
  behaviorName=basicParamsStruct.behaviors.names;
  labelColors=basicParamsStruct.behaviors.labelcolors;
  movieFileName=basicParamsStruct.file.moviefilename;
  movieIndexFileName=basicParamsStruct.file.movieindexfilename;
  trackFileName=basicParamsStruct.file.trxfilename;
  scoreFileName=basicParamsStruct.file.scorefilename;
  wfparams=basicParamsStruct.windowFeaturesParams;
  cstuff=basicParamsStruct.classifierStuff;
  everythingParams = Macguffin(basicParamsStruct.featureLexiconName,basicParamsStruct.aptInfo);
  everythingParams.behaviors.names=behaviorName;
  everythingParams.behaviors.labelcolors=labelColors;
  everythingParams.file.moviefilename=movieFileName;
  everythingParams.file.movieindexfilename=movieIndexFileName;
  everythingParams.file.trxfilename=trackFileName;
  everythingParams.file.scorefilename=scoreFileName;
  everythingParams.windowFeaturesParams=wfparams;
  everythingParams.classifierStuff=cstuff;
else
  everythingParams = Macguffin(basicParamsStruct);
end
assert(~isempty(handles.handleobj));
handles.handleobj.data = everythingParams;

% Delete the ProjectSetup window
delete(get(hObject,'parent'));

% -------------------------------------------------------------------------
% --- Executes on button press in pushbutton_copy.
function pushbutton_copy_Callback(hObject, eventdata, handles) %#ok
% hObject    handle to pushbutton_copy (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[fnames,pname] = uigetfile('*.jab','Select a project to copy settings from');
if fnames == 0; return; end;

list = {'Target Type','Behavior Name and File Names','Window Features',...
  'List of perframe features','Classifier files used as input',...
  'Classifier Settings',...
  'Advanced Parameters'};

dstr = {'Select the parameters to copy.','Use Control/Command click to select multiple entries'};
[sel,ok] = listdlg('PromptString',dstr,'ListSize',[350 120],'Name','Parameters to copy',...
  'ListString',list);
sellist = list(sel);

if ok == 0, return; end

try

  origparams = load(fullfile(pname,fnames),'-mat');

if ismember('Target Type',sellist)
    handles.basicParamsStruct.featureLexiconName = origparams.x.featureLexiconName;
    handles.basicParamsStruct.featureLexicon = origparams.x.featureLexicon;
end

if ismember('Behavior Name and File Names',sellist)
  behnames = origparams.x.behaviors.names;
  isNone=strcmpi('none',behnames);
  behnames=behnames(~isNone);
  handles.basicParamsStruct.behaviors.names = behnames;
  handles.basicParamsStruct.file = origparams.x.file;
end

if ismember('List of perframe features',sellist)
  if ~isequal(handles.basicParamsStruct.featureLexicon,origparams.x.featureLexicon),
    res = questdlg(['Target type (or one of the target settings) are not the same for the current project '...
        'and the original project. Are you sure you want to import the list of perframe features?'],...
        'Import list of perframe features','Yes','No','No');
      if strcmp(res,'Yes')
        handles.basicParamsStruct.sublexiconPFNames= origparams.x.sublexiconPFNames;
      end
  else
    handles.basicParamsStruct.sublexiconPFNames= origparams.x.sublexiconPFNames;

  end
end

if ismember('Window Features',sellist),
  if ~isequal(handles.basicParamsStruct.featureLexicon,origparams.x.featureLexicon) || ...
    ~isequal(handles.basicParamsStruct.sublexiconPFNames,origparams.x.sublexiconPFNames)
    res = questdlg(['Target type or list of perframe features (or one of the traget settings) are not the same for the current project '...
      'and the original project. Are you sure you want to import the window features?'],...
      'Import Window features','Yes','No','No');
    if strcmp(res,'Yes')
      handles.basicParamsStruct.windowFeaturesParams = origparams.x.windowFeaturesParams;
    end
  else
      handles.basicParamsStruct.windowFeaturesParams = origparams.x.windowFeaturesParams;
  end
end

if ismember('Classifier files used as input', sellist);
  handles.basicParamsStruct.scoreFeatures = origparams.x.scoreFeatures;
end

if ismember('Classifier Settings',sellist),
  clf_params_specific = {'classifierStuff.type'
    'classifierStuff.postProcessParams'
    'classifierStuff.trainingParams'
    'classifierStuff.savewindowdata'
    'classifierStuff.selFeatures'
    'classifierStuff.predictOnlyCurrentFly'
    };

  for str = clf_params_specific(:)',
    try %#ok<TRYNC>
      eval(sprintf('handles.basicParamsStruct.%s = origparams.x.%s;',str{1},str{1}));
    end
  end

end

if ismember('Advanced Parameters',sellist),
  adv_params_specific = {'behaviors.labelcolors'
    'behaviors.unknowncolor'
    'trxGraphicParams'
    'extra'};

  for str = adv_params_specific(:)',
    try %#ok<TRYNC>
      eval(sprintf('handles.basicParamsStruct.%s = origparams.x.%s;',str{1},str{1}));
    end
  end

end

guidata(hObject,handles);
updateEditsListboxesAndPopupmenus(handles);
updateConfigTable(handles);

catch ME,
  uiwait(warndlg(sprintf('Could not copy settings from previous project:%s',ME.message)));
end
return

% -------------------------------------------------------------------------
% --- Executes on button press in togglebutton_advanced.
function togglebutton_advanced_Callback(hObject, eventdata, handles)  %#ok
handles.mode=fif(get(hObject,'Value'),'advanced','basic');
handles = updateFigurePosition(handles);
guidata(hObject,handles);
return

% -------------------------------------------------------------------------
% --- Executes on button press in pushbutton_perframe.
function pushbutton_perframe_Callback(hObject, eventdata, handles)  %#ok
featureLexiconName=handles.basicParamsStruct.featureLexiconName;
sublexiconPFNames = handles.basicParamsStruct.sublexiconPFNames;

if ~get(handles.radiobutton_apt,'Value')
  featureLexicon=featureLexiconFromFeatureLexiconName(featureLexiconName);
  [featureLexiconPFNames,isSelected,missingPFNames] = ...
    GetAllPerframeList(featureLexicon,sublexiconPFNames);

  % If some to-be-calculated PFs are not in the lexicon, warn the user.
  if ~isempty(missingPFNames),
    list = naturalLanguageListFromStringList(missingPFNames);
    if length(missingPFNames)==1
      wstr = ...
        sprintf(['Perframe feature %s is defined in the project file ' ...
                 'but is not defined for the current target type.\n' ...
                 'Ignoring it.'],...
                list);
    else
      wstr = ...
        sprintf(['Perframe features %s are defined in the project file ' ...
                 'but are not defined for the current target type.\n' ...
                 'Ignoring them.'],...
                list);
    end
    uiwait(warndlg(wstr));
  end
else
  featureLexicon = handles.basicParamsStruct.featureLexicon;
  [featureLexiconPFNames,isSelected,missingPFNames] = ...
    GetAllPerframeList(featureLexicon,sublexiconPFNames);
  assert(isempty(missingPFNames),'No PF feature should be missing');
end
% Put up the dialog that allows users to pick which PFs will be computed

[isSelected,ok] = listdlg('ListString',featureLexiconPFNames, ...
                          'InitialValue',find(isSelected), ...
                          'Name','Selecte perframe features',...
                          'PromptString',{sprintf('%s-click to',fif(ismac(),'Command','Control')), ...
                                          'select/deselect perframe features'}, ...
                                          'ListSize',[250,700]);

% If the user clicked OK, update the list of to-be-calculated PFs
if ok,
  sublexiconPFNames=featureLexiconPFNames(isSelected);
  handles.basicParamsStruct.sublexiconPFNames = sublexiconPFNames;
  guidata(hObject,handles);
end

return

% -------------------------------------------------------------------------
function [featureLexiconPFNames,isSelected,missingPFNames]= ...
  GetAllPerframeList(featureLexicon,sublexiconPFNames)
% Returns the list of all per-frame feature names in the lexicon indicated
% by featureLexiconName.  Also returns isSelected, a boolean array of the
% same size as featureLexiconPFNames, which is true iff that PF is in the
% sublexiconPFNames list.  Finally, missingPFNames is a list of all PF
% names in sublexiconPFNames that are not present in
% featureLexiconPFNames.

%featureLexiconName=handles.basicParamsStruct.featureLexiconName;
%featureconfigfile = basicParamsStruct.file.featureconfigfile;
%basicParamsStruct = ReadXMLParams(featureconfigfile);
%allPfList = fieldnames(basicParamsStruct.perframe);
featureLexiconPFNames = fieldnames(featureLexicon.perframe);
nFeatureLexiconPFNames=length(featureLexiconPFNames);
isSelected = false(nFeatureLexiconPFNames,1);
missingPFNames = {};
%toBeCalculatedPFNames = basicParamsStruct.toBeCalculatedPFNames;
if isempty(sublexiconPFNames),
  isSelected = true(nFeatureLexiconPFNames,1);
  missingPFNames = {};
else
  for ndx = 1:length(sublexiconPFNames),
    sublexiconPFName=sublexiconPFNames{ndx};
    allndx = find(strcmp(sublexiconPFName,featureLexiconPFNames));
    if isempty(allndx)
      missingPFNames{end+1} = sublexiconPFName;  %#ok
    else
      isSelected(allndx) = true;
    end
  end
end
return

% -------------------------------------------------------------------------
% --- Executes on button press in pushbutton_addadvanced.
function pushbutton_addadvanced_Callback(hObject, eventdata, handles)  %#ok
in = inputdlg({'Configuration Parameter Name','Configuration Parameter Value'});
if isempty(in) || length(in) < 2
  return;
end

handles = AddConfig(handles,in{1},in{2});
updateConfigTable(handles);
updateEditsListboxesAndPopupmenus(handles);
guidata(hObject,handles);

% -------------------------------------------------------------------------
% --- Executes on button press in pushbutton_removeadvanced.
function pushbutton_removeadvanced_Callback(hObject, eventdata, handles)  %#ok
data = get(handles.config_table,'Data');

jUIScrollPane = findjobj(handles.config_table);
jUITable = jUIScrollPane.getViewport.getView;
allndx = jUITable.getSelectedRows + 1;
if numel(allndx)==1 && allndx <1, return, end

for ndx = allndx(:)'
  handles = RemoveConfig(handles,data{ndx,1});
end
updateConfigTable(handles);
updateEditsListboxesAndPopupmenus(handles);
guidata(hObject,handles);
return

% -------------------------------------------------------------------------
% --- Executes when entered data in editable cell(s) in config_table.
function config_table_CellEditCallback(hObject, eventdata, handles)  %#ok
% hObject    handle to config_table (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)
if isempty(eventdata.Indices); return; end
data = get(handles.config_table,'Data');
ndx = eventdata.Indices(1);
if eventdata.Indices(2) == 2
  handles = EditConfigValue(handles,data{ndx,1},data{ndx,2});
else
  handles = EditConfigName(handles,eventdata.PreviousData,eventdata.NewData);
end
updateEditsListboxesAndPopupmenus(handles);
guidata(hObject,handles);

% -------------------------------------------------------------------------
function handles = RemoveConfig(handles,name)

[fieldNamePath,dotLastField] = splitext(name);
if isempty(dotLastField)
  handles.basicParamsStruct = rmfield(handles.basicParamsStruct,fieldNamePath);
else
  lastField=dotLastField(2:end);  %#ok
  evalStr = sprintf(...
    'handles.basicParamsStruct.%s = rmfield(handles.basicParamsStruct.%s,lastField);',...
    fieldNamePath,fieldNamePath);
  eval(evalStr);
end

% -------------------------------------------------------------------------
function handles = EditConfigName(handles,oldName,newName)
eval_str = sprintf(...
  'value = handles.basicParamsStruct.%s;',...
  oldName);
eval(eval_str);
handles = AddConfig(handles,newName,value);
handles = RemoveConfig(handles,oldName);
return

% -------------------------------------------------------------------------
function handles = EditConfigValue(handles,name,value)

switch name
  case 'behaviors.names'
    value = regexp(value,',','split'); %#ok<NASGU>
  otherwise
end

% try to convert it to a number if possible
valuen = regexp(value,',','split');
valuen = str2double(valuen);
if all(~isnan(valuen)),
  value = valuen; %#ok<NASGU>
end
eval_str = sprintf('handles.basicParamsStruct.%s = value;',name);
eval(eval_str);
return

% -------------------------------------------------------------------------
function handles = AddConfig(handles,name,value)

% If value is a string, convert to a double matrix, if that works
if ischar(value) && ~isempty(str2num(value)),  %#ok<ST2NM>
  value = str2num(value);  %#ok
end
evalStr=sprintf('handles.basicParamsStruct.%s = value;',name);
eval(evalStr);

return

% -------------------------------------------------------------------------
function figureProjectSetup_CloseRequestFcn(hObject, eventdata, handles)  %#ok
% Hint: delete(hObject) closes the figure
%uiresume(hObject);
pushbutton_cancel_Callback(hObject, eventdata, handles)
return

% -------------------------------------------------------------------------
% --- Executes on key press with focus on figureProjectSetup and none of its controls.
function figureProjectSetup_WindowKeyPressFcn(hObject, eventdata, handles)  %#ok
% hObject    handle to figureProjectSetup (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)

if isequal(eventdata.Key,'f') && isControlLike(eventdata.Modifier)
  handles.basicParamsStruct.file.moviefilename = 'movie.ufmf';
  handles.basicParamsStruct.file.trxfilename = 'registered_trx.mat';
  guidata(hObject,handles);  % write the changes to the figure
  updateEditsListboxesAndPopupmenus(handles);
end

return

% --- Executes when selected object is changed in uibuttongroup1.
function uibuttongroup1_SelectionChangedFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in uibuttongroup1
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~get(handles.radiobutton_apt,'Value')
  set(handles.featureconfigpopup,'Enable','on');
  set(handles.text_apt,'Visible','Off');
  set(handles.pushbutton_edit_APT,'Visible','off');
  featureconfigpopup_Callback(handles.featureconfigpopup, eventdata, handles)
  return;
end
% Do stuff if init from APT project is selected
if isfield(handles.basicParamsStruct,'lbl_file') && ...
    ischar(handles.basicParamsStruct.lbl_file),
  defaultlbl_file = handles.basicParamsStruct.lbl_file;
else
  defaultlbl_file = '';
end
debug = false;
if debug
  Q = load('~/temp/aptprojdata');
  aptStruct = Q.qq;
else
  [fname,fpath] = uigetfile('*.lbl','Select APT project File..',defaultlbl_file);
  if fname == 0
    set(handles.radiobutton_apt,'Value',0);
    set(handles.radiobutton_list,'Value',1);
    return;
  end
  % Get the detail of APT
  lbl_file = fullfile(fpath,fname);
  aptStruct = APTProject('lbl_file',lbl_file);
end

if isempty(aptStruct.featureLexicon)
  set(handles.radiobutton_apt,'Value',0);
  set(handles.radiobutton_list,'Value',1);
  return;
end
% Update the GUI.
set(handles.featureconfigpopup,'Enable','off');
set(handles.text_apt,'Visible','On');
set(handles.text_apt,'String',aptStruct.featureLexiconName);
set(handles.pushbutton_edit_APT,'Visible','on');

% Save the current values of the main four fields
% (Note that any changes in the advanced part of the window will be lost.)
behaviorName=handles.basicParamsStruct.behaviors.names;
movieFileName=handles.basicParamsStruct.file.moviefilename;
movieIndexFileName=handles.basicParamsStruct.file.movieindexfilename;
trackFileName=handles.basicParamsStruct.file.trxfilename;
scoreFileName=handles.basicParamsStruct.file.scorefilename;
% Replace basicParamsStruct with one appropriate to the new feature lexicon
% name
old=warning('query','MATLAB:structOnObject');
warning('off','MATLAB:structOnObject');  % turn off annoying warning
handles.basicParamsStruct = struct(Macguffin(handles.basicParamsStruct,aptStruct));
warning(old);  % restore annoying warning
% Now restore fields from the old basicParamsStruct
handles.basicParamsStruct.behaviors.names=behaviorName;
handles.basicParamsStruct.file.moviefilename=movieFileName;
handles.basicParamsStruct.file.movieindexfilename=movieIndexFileName;
handles.basicParamsStruct.file.trxfilename=trackFileName;
handles.basicParamsStruct.file.scorefilename=scoreFileName;
guidata(hObject,handles);
updateConfigTable(handles);
updateEditsListboxesAndPopupmenus(handles);

% --- Executes on button press in radiobutton_list.
function radiobutton_list_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_list (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_list

% --- Executes on button press in radiobutton_apt.
function radiobutton_apt_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_apt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_apt

% --- Executes on button press in checkbox_st.
function checkbox_st_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_st (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_st

if ~get(hObject,'Value')
   handles.basicParamsStruct.stFeatures = false;
   guidata(hObject,handles);
   return;
end

fname = uigetdir('.','Select an example JAABA experiment directory');
if fname == 0
  mov_file = '';
  trx_file = '';
  set(hObject,'Value',0);
  return
else
  mov_file = fullfile(fname,handles.basicParamsStruct.file.moviefilename);
  trx_file = fullfile(fname,handles.basicParamsStruct.file.trxfilename);
  if ~exist(mov_file,'file')
      errmsg = sprintf('Movie file %s does not exist',mov_file);
      errordlg(errmsg);
      set(hObject,'Value',0);
      return;
  end
  if ~exist(trx_file,'file'),trx_file=''; end
end

set(handles.figureProjectSetup, 'pointer', 'watch')
[useST,stInfo] = spaceTimeOptions('mov_file',mov_file,...
                    'trx_file',trx_file,...
                    'prev_st',handles.basicParamsStruct.stInfo);

set(handles.figureProjectSetup, 'pointer', 'arrow')
if useST
    handles.basicParamsStruct.stFeatures = true;
    handles.basicParamsStruct.stInfo = stInfo;
else
    handles.basicParamsStruct.stFeatures = false;
    set(hObject,'Value',0);
end
guidata(hObject,handles);

% --- Executes on button press in pushbutton_edit_APT.
function pushbutton_edit_APT_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_edit_APT (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Do stuff if init from APT project is selected
[fname,fpath] = uigetfile('*.lbl','Select APT project File..',handles.basicParamsStruct.aptInfo.lbl_file);
if fname == 0,
  return;
end
% Get the detail of APT
lbl_file = fullfile(fpath,fname);
aptStruct = APTProject('lbl_file',lbl_file,'aptStruct',handles.basicParamsStruct.aptInfo);
if isempty(aptStruct.featureLexicon)
  set(handles.radiobutton_apt,'Value',0);
  set(handles.radiobutton_list,'Value',1);
  return;
end

author: kristinbranson

repo: JAABA


text:
%% DEMO 7:  Algorithms 02. SART
%
%
% In this demo the usage of the algorithms on the SART family among with
% their options are presented.
%
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
%% Initialize
clear;
close all;

%% Define Geometry
geo=defaultGeometry('nVoxel',[128;128;128]);

%% Load data and generate projections
% see previous demo for explanation
angles=linspace(0,2*pi,100);
head=headPhantom(geo.nVoxel);
projections=Ax(head,geo,angles,'interpolated');
noise_projections=addCTnoise(projections);
%% SART family of algorithms
%
% There are 3 algorithms in this damily included in TIGRE: SART,SIRT and
% OS-SART.
%
% The main difference between them is the update process.
%   SART: Updates the image projection by projection
%   SIRT: Updates the image using the whole set of projections at once
%   OS-SART: middle ground. Updates the image using a subset of the
%            projections
%
%  Of these algorithms, SART is generally the one reaching a better image
%  (less L2 error) for the same amount of iterations, and SIRT is the
%  worst (still relatively similar). However, SART needs increased
%  computational time per iteration, as it needs to update the image very often,
%  while SIRT only updates the emage ones for the whole sets of projections.
%  OS-SART lies in the middle, reaching similar convergence (L2 error per
%  iteration) than SART but with less computational time than SART.
%
%% Usage, with optional parameters.
% In the three algorithms, there are 4 mandatory input arguments:
% Projections, geometry, angles and number of iterations.
%
%
% Optional arguments for all of them
%==========================================================================
% 'lambda': hyperparameter. The update will be multiplied by this number
% every iteration, to make the steps bigger or smaller. Default: 1
%
lambda=1;

% 'lambdared': reduction multiplier for the hyperparameter.
% lambda=lambda*lambdared every iterations, so the steps can be smaller
% the further the update. Default=0.99
lambdared=0.999;

% 'Init' : Initialization method. Possible options are
%          'none' (default). There will be no initialization method, just
%                 the algorithm
%
%          'FDK'  Initialize the image with the result of FDK algorithm
%
%          'multigrid' Initialize using the multigrid method. The image
%                      will be solved in a small scale, and the size of it
%                      will increase until the desired size is reached.
%
%          'image'     Initialzies with a user given image. Not recoomended
%                      unless you really know what you are doing.

initmode='none';

% 'InitImg' : related to init. The image to use for initializing the
% algorithm.

% 'verbose': boolean to make the algorithm display (or not) running state.
%            default true.

verbose=true;
% 'QualMeas'     Asks the algorithm for a set of quality measurement
%                parameters. Input should contain a cell array of desired
%                quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                These will be computed in each iteration.
qualmeas={'RMSE'};

% SIRT and SART both have no extra input parameters.
% =========================================================================
[imgSIRT,errL2SIRT,qualitySIRT]=SIRT(projections,geo,angles,20,...
                            'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas);
[imgSART,errL2SART,qualitySART]=SART(projections,geo,angles,20,...
                            'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas);
% OS-SART
% ========================================================================
% Additionally OS-SART includes a couple of other parameters, related to
% the subsets.
%
%   'BlockSize':   Sets the projection block size used simultaneously. If
%                  BlockSize = 1 OS-SART becomes SART and if  BlockSize = size(angles,2)
%                  then OS-SART becomes SIRT. Default is 20.
blcks=8;
% 'OrderStrategy':  Chooses the subset ordering strategy. Options are
%                  'ordered' :uses them in the input order, but divided
%                  'random'  : orders them randomply
%                  'angularDistance': chooses the next subset with the
%                                     biggest angular distance with the
%                                     ones used.  (default)
order='angularDistance';
[imgOSSART,errL2OSSART,qualityOSSART]=OS_SART(projections,geo,angles,20,...
                            'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas,...
                             'BlockSize',blcks,'OrderStrategy',order);
%% Lets have a brief show of the results
% set(0,'DefaultTextInterpreter', 'latex')

subplot(211)
plot(log10([errL2SIRT;[errL2OSSART nan(1,length(errL2SIRT)-length(errL2OSSART))];[errL2SART nan(1,length(errL2SIRT)-length(errL2SART))]]'));
title('Convergence')
xlabel('Iteration')
ylabel('$ log_{10}(|Ax-b|) $','interpreter','latex')
legend('SIRT','OS-SART','SART')
subplot(212)
plot(log10([qualitySIRT;[qualityOSSART nan(1,length(qualitySIRT)-length(qualityOSSART))];[qualitySART nan(1,length(qualitySIRT)-length(qualitySART))]]'));
title('Evolution of RMSE')
legend('SIRT','OS-SART','SART')
xlabel('Iteration')
ylabel('$ log_{10}(RMSE) $','interpreter','latex')

%% plot the results

% It is clear that SART will get to better results for the same amoutn of
% iterations, however, it takes x7 more time to run.

% SART
% OS-SART
% SIRT

plotImg([imgSIRT;  imgOSSART; imgSART;],'Dim','Z','Savegif','sarts.gif');

% plot error
plotImg(abs([head-imgSIRT; head-imgOSSART; head-imgSART; ]),'Dim','Z');

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function y = log_normcdf( x )

%LOG_NORMCDF   Internal CVX version.

error(nargchk(1,1,nargin));
if ~isreal( x ),
    error( 'Argument must be real.' );
end

persistent a b nb ob
if isempty( a ),
    a =sqrt( [ 0.018102332171520
               0.011338501342044
               0.072727608432177
               0.184816581789135
               0.189354610912339
               0.023660365352785 ] );
    a = sparse(diag(a));
    b = [3 2.5 2 1 -1 -2]';
    nb = length(b);
    ob = ones(nb,1);
end

cx = cvx_isconstant( x );
sx = size(x);
nx = prod(sx);
y  = a * ( b * ones(1,nx) - ob * reshape( x, 1, nx ) );
if cx,
    y = cvx( sum( cvx_constant( max( y, 0 ) ) .^ 2 ) );
else
    y = sum_square_pos( y );
end
y = - reshape( y, sx );

repo: REST


text:
function point_num = sphere_llq_grid_point_count ( lat_num, long_num )

%*****************************************************************************80
%
%% SPHERE_LLQ_GRID_POINT_COUNT counts points for a SPHERE LLQ grid.
%
%  Discussion:
%
%    A SPHERE LLQ grid imposes a grid of quadrilaterals on a sphere,
%    using latitude and longitude lines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LAT_NUM, LONG_NUM, the number of latitude
%    and longitude lines to draw.  The latitudes do not include the North and
%    South poles, which will be included automatically, so LAT_NUM = 5, for
%    instance, will result in points along 7 lines of latitude.
%
%    Output, integer POINT_NUM, the number of grid points.
%
  point_num = 2 + lat_num * long_num;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mu = meanmv(X, dim);
% handles missing values (NaN).

error(nargchk(1,2,nargin));

if nargin < 2
  dim = min(find(size(X)~=1));
  if isempty(dim), dim = 1; end
end

% NaNs to zero
Imv = isnan(X);
X(Imv) = 0;

% Sum and divide by the number of non-missing values
div = sum(~Imv, dim);
div(div==0) = 1;
mu = sum(X, dim) ./ div;
%mu(isinf(xmean)) = NaN;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
%
% [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
%
% convert a binary (or multi-valued) volume to tetrahedral mesh
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 img: a volumetric binary image
%	 ix,iy,iz: subvolume selection indices in x,y,z directions
%	 opt: as defined in vol2surf.m
%	 maxvol: target maximum tetrahedral elem volume
%                when method='cgalmesh', maxvol can specify the target
%                for each label (subregion index) by the following syntax
%                'label1=size1:label2=size2:...'
%	 dofix: 1: perform mesh validation&repair, 0: skip repairing
%	 method: 'cgalsurf' or omit: use CGAL surface mesher
%		 'simplify': use binsurface and then simplify
%		 'cgalmesh': use CGAL 3.5 3D mesher for direct mesh generation [new]
%
%		 generally speaking, 'cgalmesh' is the most robust path
%		 if you want to product meshes from binary or multi-region
%		 volumes, however, its limitations include 1) only accept
%		 uint8 volume, and 2) can not extract meshes from gray-scale
%		 volumes. If ones goal is to process a gray-scale volume,
%		 he/she should use the 'cgalsurf' option. 'simplify' approach
%		 is not recommended unless other options has failed.
%	 isovalues: a list of isovalues where the levelset is defined
%
% output:
%	 node: output, node coordinates of the tetrahedral mesh
%	 elem: output, element list of the tetrahedral mesh, the last
%	       column is the region ID
%	 face: output, mesh surface element list of the tetrahedral mesh
%	       the last column denotes the boundary ID
%        region: optional output. if opt.autoregion is set to 1, region
%              saves the interior points for each closed surface component
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin>=8)
	if(strcmp(method,'cgalmesh'))
		vol=img(ix,iy,iz);
		if(length(unique(vol(:)))>64 & dofix==1)
			error([ 'it appears that you are processing a ' ...
                                'grayscale image. Currently cgalmesher ' ...
                                'does not support grayscale images. ' ...
                                'Please use "cgalsurf" method to mesh a grayscale ' ...
                                'volume. If you are certain to run cgalmesher ' ...
                                'on your data, please set dofix=0 and run this again.' ]);
		end
		[node elem,face]=cgalv2m(vol,opt,maxvol);
		return;
	end
end

%first, convert the binary volume into isosurfaces
if(nargin==8)
	[no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method);
elseif(nargin==9)
	[no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues);
else
        [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,'cgalsurf');
end
%then, create volumetric mesh from the surface mesh
if(nargin>=8)
   if(strcmp(method,'cgalpoly'))
	[node,elem,face]=cgals2m(no(:,1:3),el(:,1:3),opt,maxvol);
        return;
   end
end

[node,elem,face]=surf2mesh(no,el,[],[],1,maxvol,regions,holes);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Filtering A Delta Function Input Signal Example Part 2
%
% This example illustrates how numerical aliasing can be avoided by
% spatially smoothing the source mask.
%
% author: Bradley Treeby
% date: 19th January 2010
% last update: 24th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 256;           % number of grid points in the x (row) direction
dx = 10e-3/Nx;      % grid point spacing in the x direction [m]
kgrid = makeGrid(Nx, dx);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]

% create a time array
num_time_steps = 1024;
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
kgrid.t_array = 0:dt:dt*(num_time_steps - 1);

% define a single element source
source_offset = 50;
source.p_mask = zeros(Nx, 1);
source.p_mask(1 + source_offset, 1) = 1;

% spatially smooth the source mask maintaining the maximum magnitude
source.p_mask = smooth(kgrid, source.p_mask, true);

% threshold out small values
source.p_mask(source.p_mask < 0.05) = 0;

% define a delta function input pulse
temporal_offset = 100;      % [time steps]
source_magnitude = 2;       % [au]
source_func = zeros(size(kgrid.t_array));
source_func(temporal_offset) = source_magnitude;

% assign and scale the input pulse based on the source mask
source.p(1:sum(source.p_mask ~= 0),:) = source.p_mask(source.p_mask ~= 0)*source_func;

% force the source mask to be binary
source.p_mask(source.p_mask ~= 0) = 1;

% define a single element sensor
sensor.mask = zeros(Nx, 1);
sensor.mask(end - source_offset, 1) = 1;

% run the simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PMLSize', 30);

% compute the amplitude spectra of the recorded time series
[f, output_as] = spect(sensor_data, 1/dt);

% extract the maximum frequency supported by the grid (two points per
% wavelength)
f_max = kgrid.k_max * min(medium.sound_speed(:)) / (2*pi);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the input and recorded time series
figure;
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
plot(kgrid.t_array*scale, source.p(floor(sum(source.p_mask(:))/2) + 1, :), 'k-', kgrid.t_array*scale, sensor_data, 'b-');
xlabel(['Time [' prefix 's]']);
ylabel('Pressure [au]');
legend('input pulse', 'recorded pulse');

% plot the amplitude spectrum
[f_sc, scale, prefix] = scaleSI(max(f));
figure;
plot(f*scale, output_as, 'b-');
xlabel(['Frequency [' prefix 'Hz]']);
ylabel('Amplitude [au]');

% plot the maximum frequency supported by the grid
ylim = get(gca, 'YLim');
hold on;
line([f_max*scale, f_max*scale], [0, ylim(2)], 'LineStyle','--', 'Color', 'k');
legend('amplitude spectrum of recorded pulse', 'maximum frequency supported by grid');

repo: Aero-matlab


text:
function A = myspconvert(A, nrows, ncols, tol)

% function A = myspconvert(A, nrows, ncols, tol)
% Purpose: convert a triplet sparse matrix into a Matlab sparse
%          matrix and crop entries smaller than tol

ids = find(abs(A(:,3))>tol);
A = A(ids, :);

% convert to sparse matrix
A = spconvert(A);

% pad to required size
if(size(A,1)<nrows | size(A,2)<ncols)
A(nrows,ncols) = 0;
end
return


text:
% Test file for chebfun constructor (basic).

function pass = test_constructor_basic(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% Some basic test functions:
FF = {@sin, @(x) [sin(x), cos(x)], @(x) [sin(x), cos(x), exp(-x)]};

for j = 1:numel(FF);
    % Initialise k:
    k = 0;

    % Pick the test function:
    F = FF{j};

    % Test on [-1 1]:
    f = chebfun(F, [-1, 1], pref);
    xx = linspace(-1, 1);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 10*eps*vscale(f);
    pass(j, k+2) = err < 50*pref.chebfuneps;
    k = k + 2;

    % Test on [-1 1] (no domain passed):
    f = chebfun(F, pref);
    xx = linspace(-1, 1);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 10*eps*vscale(f);
    pass(j, k+2) = err < 500*pref.chebfuneps;
    k = k + 2;

    % Test on [0 10000]:
    f = chebfun(F, [0, 10000], pref);
    xx = linspace(0, 10000);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 1e4*eps*vscale(f);
    pass(j, k+2) = err < 1e2*hscale(f)*pref.chebfuneps;
    k = k + 2;

    % Test on piecewise domain:
    f = chebfun(F, [-1, 0, .5, sqrt(pi/4), 1], pref);
    xx = linspace(-1, 1);
    err = norm(feval(f, xx) - F(xx), inf);
    pass(j, k+1) = err < 10*eps*vscale(f);
    pass(j, k+2) = err < 100*pref.chebfuneps;
    k = k + 2;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ex_ssc_omp(md, num_samples_per_digit)
    if nargin < 1
        error('database must be loaded.');
    end
    if nargin < 2
        num_samples_per_digit = 50;
    end
    close all;
    clc;
    check_digits(md, num_samples_per_digit, @ssc_omp, 'ssc_omp');
end

author: indigits

repo: sparse-plex


text:
function qwgw_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests QWGW for the Chebyshev Type 1 weight.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2014
%
%  Author:
%
%    John Burkardt
%

%
%  Set the quadrature interval and number of points.
%
  a = -1.0;
  b = +1.0;
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  Compute points and weights for Gauss quadrature\n' );
  fprintf ( 1, '  with the Chebyshev Type 1 weight w(x) = 1/sqrt(1-x^2).\n' );
  fprintf ( 1, '  Order N = %d\n', n );
  fprintf ( 1, '  Interval = [%g,%g]\n', a, b );
%
%  Set the recursion coefficients.
%
  aj = zeros ( n, 1 );
  bj = zeros ( n, 1 );

  aj(1:n) = 0.0;

  bj(1) = 1.0 / 2.0;
  for j = 2 : n - 1
    bj(j) = 1.0 / 4.0;
  end
  bj(n) = 0.0;

  bj(1:n) = sqrt ( bj(1:n) );

  mu0 = pi;
%
%  Compute the points and weights.
%
  [ x, w ] = sgqf ( n, aj, bj, mu0 );

  r8vec_print ( n, x, '  Abscissas:' );
  r8vec_print ( n, w, '  Weights:' );

  return
end


text:
%%
% \documentclass[12pt]{article}
%
% \title{ODEbox: A Toolbox for Ordinary Differential Equations\\
% Sturm Liouville Example 2\\
% The Mathieu Differential Equation}
%
% \author{Matthew Harker and Paul O'Leary\\
% Institute for Automation\\
% University of Leoben\\
% A-8700 Leoben,
% Austria\\
% URL: automation.unileoben.ac.at\\
% \\
% Original: January 9, 2013\\
% $\copyright$ 2013\\
% \\
% Last Modified: \today}
%%
% \section{The Mathieu Differential Equation}
%
% This example addresses the solution of the Mathieu differential equation.
% This equation occurs in conjunction with the modelling of the vibration
% of a ellptical membrane. The problem has no known analytical solution.
% %
% \begin{equation}
%    -\ddot{y} + 2 \,r \, \cos( 2 \, x )= \lambda \, y
%    \hspace{5mm}
%    \text{with}
%    \hspace{5mm}
%    y(0) = 0}
%   \hspace{5mm}
%   \text{and}
%   \hspace{5mm}
%   \dot{y}(\pi) = 0
% \end{equation}
% %
% in the closed interval $0 \leq x \leq \pi$.
%
%%
% \section{tidy up the Workspace}
%
clear all;
close all;
setUpGraphics;
%%
% \section{Define the Nodes and Compute the Basis Functions}
%
% Define the number of nodes and the number of basis functions used.
%
noPts = 1000;
noBfs = round( noPts/2);
%%
% Compute the Chebyshev points, but scaled to form a closed interval
% $0 \leq x \leq \pi$.
%
x = dopNodes( noPts, 'chebyends' );
x = pi * (x + 1)/2;
%%
% Synthesize a complete set of basis functions
%
B = dop(x);
%%
% \subsection{Generate the Differentiating Matrix}
%
% Synthesize the differentiating matrix with support length $l_s$
%
ls = 13;
D = dopDiffLocal( x, ls, ls );
%%
% \subsection{Define the Constraints and Compute the Admissible Functions}
%
% Define the constraints
%
C = zeros( noPts, 2 );
C(1,1) = 1;
C(end,end) = 1;
%%
% Compute the set of constrained basis functions, i.e., admissible
% functions.
%
Bc = dopConstrain( C, B );
%%
% Trucate the basis functions
%
Bc = Bc(:,1:noBfs);
%
%%
% \section{Setup the Sturm-Liouville Matrix Linear Differential Operator}
%
% Setup the Linear differential operator for the Mathieu differential
% equation.
%
c1 = 0;
c2 = -50;
gx = diag(c1 + c2 * cos(2*x)) ;
L = Bc' * (D * D - gx )* Bc;
%
%%
% \section{Solve the Eigenvector Problem}
%
% Solve the eigenvalue problem
%
[Vec, Val] = eig( L );
vals = -diag( Val );
%%
% and sort the solutions.
%
[vals, inds] = sort(vals);
Vec = Vec(:,inds);
%%
% and compute the final eigenfunctions.
%
sols = Bc * Vec;
%%
% \section{Close Pair of Eigenvalues}
%
% The Mathieu differential equation as paramatized here is known to produce
% a close pair of eigenvalues. We now show this pair.
%
noEVals = 2;
for k=1:noEVals
    numStr = num2str(vals(k),'%2.10E') ;
    str = ['\lambda_{', int2str(k-1), '} = ', numStr]
end;
%noV = length(vals);
%n = [1:noV]';
%valsT = n.^2;
%
%%
% \section{Plot a Few Eigenfunctions}
%
setUpGraphics(10)
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.18 0.17 0.8 0.8];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig1 = figure;
plot( x, sols(:,1), 'r');
hold on
plot( x, sols(:,2), 'g');
plot( x, sols(:,3), 'b');
plot( x, sols(:,4), 'k');
grid on;
xlabel('$$x$$');
ylabel('$$y8x)$$');
axis([0,pi,-0.1,0.1]);
%
% \caption{The first $4$ eigen functions of the Mathieu differential equation.}
%%
%
setUpGraphics(10)
FigureSize=[1 1 10 10];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.16 0.17 0.8 0.8];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig2 = figure;
P1 = [0.16 0.4 0.8 0.55] ;
A = axes('position',P1);
imagesc( log10(abs(Vec) ));
colorbar;
ylabel(['Basis function number $$m$$, $$n_b = ',int2str(noBfs),'$$']);
P2 = [0.16 0.1 0.615 0.25] ;
A = axes('position',P2);
plot( vals(1:noBfs), 'k');
range = axis;
axis([0,noBfs,-1000,900000]);
grid on;
xlabel('Eigenvector number $$n$$');
ylabel('$$\lambda$$');
%
% \caption{The Rayleigh-Ritz Spectrum of the Eigenfunvtions with respect to
% the admissible functions and the eigenvalues for the Mathieu differential
% equation.}
%%
% \section{Save the Figure to Disk.}
fileType = 'eps';
printFigure( fig1, ['SLEx2EigFnsNb',int2str(noBfs)], fileType);
printFigure( fig2, ['SLEx2SpectrumNb',int2str(noBfs)], fileType);

author: Sable

repo: mcbench-benchmarks


text:
function [elc, lab] = elec1020_locate(pnt, dhk, nas, ini, lpa, rpa, feedback)

% ELEC1020_LOCATE determines 10-20 (20%, 10% and 5%) electrode positions
% on a scalp surface that is described by its surface triangulation

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

if nargin<7
  feedback = false;
end

% determine the approximate location of the vertex
ori = (lpa+rpa+nas+ini)/4;      % center of head
ver =  cross(rpa-lpa, nas-ini); % orientation
ver = ver /sqrt(norm(ver));     % make correct length
ver = ori + 0.7*ver;            % location from center of head

if feedback
  figure
  ft_plot_mesh(struct('pos', pnt, 'tri', dhk), 'edgecolor', 'none', 'facecolor', 'skin')
  lighting gouraud
  material dull
  lightangle(0, 90);
  alpha 0.9
  ft_plot_mesh(nas, 'vertexsize', 30)
  ft_plot_mesh(lpa, 'vertexsize', 30)
  ft_plot_mesh(ini, 'vertexsize', 30)
  ft_plot_mesh(rpa, 'vertexsize', 30)
  ft_plot_mesh(ver, 'vertexsize', 30)
  grid on
  hold on
  view([1 1 0.5])
end

% point near LPA that is at 50% of left lower contour
[cnt1, cnt2] = elec1020_follow(pnt, dhk, nas, lpa, ini, feedback);
mle = elec1020_fraction(cnt1, cnt2, 0.5);

% point near RPA that is at 50% of right lower contour
[cnt1, cnt2] = elec1020_follow(pnt, dhk, nas, rpa, ini, feedback);
mre = elec1020_fraction(cnt1, cnt2, 0.5);

% determine two points that approximate the vertex
[cnt1, cnt2] = elec1020_follow(pnt, dhk, nas, ver, ini, feedback);
ver1 = elec1020_fraction(cnt1, cnt2, 0.5);
[cnt1, cnt2] = elec1020_follow(pnt, dhk, mle, ver, mre, feedback);
ver2 = elec1020_fraction(cnt1, cnt2, 0.5);

% refined estimate is the average of these two
ver = (ver1+ver2)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start contouring
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ant-post contour through vertex
fprintf('constructing vertical ant-post contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, nas, ver, ini, feedback);
Nz   = elec1020_fraction(cnt1, cnt2,  0/20);
NFpz = elec1020_fraction(cnt1, cnt2,  1/20);
Fpz  = elec1020_fraction(cnt1, cnt2,  2/20);
AFpz = elec1020_fraction(cnt1, cnt2,  3/20);
AFz  = elec1020_fraction(cnt1, cnt2,  4/20);
AFFz = elec1020_fraction(cnt1, cnt2,  5/20);
Fz   = elec1020_fraction(cnt1, cnt2,  6/20);
FFCz = elec1020_fraction(cnt1, cnt2,  7/20);
FCz  = elec1020_fraction(cnt1, cnt2,  8/20);
FCCz = elec1020_fraction(cnt1, cnt2,  9/20);
Cz   = elec1020_fraction(cnt1, cnt2, 10/20);
CCPz = elec1020_fraction(cnt1, cnt2, 11/20);
CPz  = elec1020_fraction(cnt1, cnt2, 12/20);
CPPz = elec1020_fraction(cnt1, cnt2, 13/20);
Pz   = elec1020_fraction(cnt1, cnt2, 14/20);
PPOz = elec1020_fraction(cnt1, cnt2, 15/20);
POz  = elec1020_fraction(cnt1, cnt2, 16/20);
POOz = elec1020_fraction(cnt1, cnt2, 17/20);
Oz   = elec1020_fraction(cnt1, cnt2, 18/20);
OIz  = elec1020_fraction(cnt1, cnt2, 19/20);
Iz   = elec1020_fraction(cnt1, cnt2, 20/20);

% left-right through vertex
fprintf('constructing C contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, mle, ver, mre, feedback);
T9   = elec1020_fraction(cnt1, cnt2,  0/20);
T9h  = elec1020_fraction(cnt1, cnt2,  1/20);
T7   = elec1020_fraction(cnt1, cnt2,  2/20);
T7h  = elec1020_fraction(cnt1, cnt2,  3/20);
C5   = elec1020_fraction(cnt1, cnt2,  4/20);
C5h  = elec1020_fraction(cnt1, cnt2,  5/20);
C3   = elec1020_fraction(cnt1, cnt2,  6/20);
C3h  = elec1020_fraction(cnt1, cnt2,  7/20);
C1   = elec1020_fraction(cnt1, cnt2,  8/20);
C1h  = elec1020_fraction(cnt1, cnt2,  9/20);
Cz   = elec1020_fraction(cnt1, cnt2, 10/20);
C2h  = elec1020_fraction(cnt1, cnt2, 11/20);
C2   = elec1020_fraction(cnt1, cnt2, 12/20);
C4h  = elec1020_fraction(cnt1, cnt2, 13/20);
C4   = elec1020_fraction(cnt1, cnt2, 14/20);
C6h  = elec1020_fraction(cnt1, cnt2, 15/20);
C6   = elec1020_fraction(cnt1, cnt2, 16/20);
T8h  = elec1020_fraction(cnt1, cnt2, 17/20);
T8   = elec1020_fraction(cnt1, cnt2, 18/20);
T10h = elec1020_fraction(cnt1, cnt2, 19/20);
T10  = elec1020_fraction(cnt1, cnt2, 20/20);

% horizontal ant-post through T7
fprintf('constructing horizontal left contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, Fpz, T7, Oz, feedback);
Fp1h = elec1020_fraction(cnt1, cnt2,  1/20);
Fp1  = elec1020_fraction(cnt1, cnt2,  2/20);
AFp7 = elec1020_fraction(cnt1, cnt2,  3/20);
AF7  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF7 = elec1020_fraction(cnt1, cnt2,  5/20);
F7   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT7 = elec1020_fraction(cnt1, cnt2,  7/20);
FT7  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT7 = elec1020_fraction(cnt1, cnt2,  9/20);
T7   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP7 = elec1020_fraction(cnt1, cnt2, 11/20);
TP7  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP7 = elec1020_fraction(cnt1, cnt2, 13/20);
P7   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO7 = elec1020_fraction(cnt1, cnt2, 15/20);
PO7  = elec1020_fraction(cnt1, cnt2, 16/20);
POO7 = elec1020_fraction(cnt1, cnt2, 17/20);
O1   = elec1020_fraction(cnt1, cnt2, 18/20);
O1h  = elec1020_fraction(cnt1, cnt2, 19/20);

% horizontal ant-post through T8
fprintf('constructing horizontal right contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, Fpz, T8, Oz, feedback);
Fp2h = elec1020_fraction(cnt1, cnt2,  1/20);
Fp2  = elec1020_fraction(cnt1, cnt2,  2/20);
AFp8 = elec1020_fraction(cnt1, cnt2,  3/20);
AF8  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF8 = elec1020_fraction(cnt1, cnt2,  5/20);
F8   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT8 = elec1020_fraction(cnt1, cnt2,  7/20);
FT8  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT8 = elec1020_fraction(cnt1, cnt2,  9/20);
T8   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP8 = elec1020_fraction(cnt1, cnt2, 11/20);
TP8  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP8 = elec1020_fraction(cnt1, cnt2, 13/20);
P8   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO8 = elec1020_fraction(cnt1, cnt2, 15/20);
PO8  = elec1020_fraction(cnt1, cnt2, 16/20);
POO8 = elec1020_fraction(cnt1, cnt2, 17/20);
O2   = elec1020_fraction(cnt1, cnt2, 18/20);
O2h  = elec1020_fraction(cnt1, cnt2, 19/20);

fprintf('constructing AFp contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, AFp7, AFpz, AFp8, feedback);
AFp7h  = elec1020_fraction(cnt1, cnt2,  1/16);
AFp5   = elec1020_fraction(cnt1, cnt2,  2/16);
AFp5h  = elec1020_fraction(cnt1, cnt2,  3/16);
AFp3   = elec1020_fraction(cnt1, cnt2,  4/16);
AFp3h  = elec1020_fraction(cnt1, cnt2,  5/16);
AFp1   = elec1020_fraction(cnt1, cnt2,  6/16);
AFp1h  = elec1020_fraction(cnt1, cnt2,  7/16);
AFp2h  = elec1020_fraction(cnt1, cnt2,  9/16);
AFp2   = elec1020_fraction(cnt1, cnt2, 10/16);
AFp4h  = elec1020_fraction(cnt1, cnt2, 11/16);
AFp4   = elec1020_fraction(cnt1, cnt2, 12/16);
AFp6h  = elec1020_fraction(cnt1, cnt2, 13/16);
AFp6   = elec1020_fraction(cnt1, cnt2, 14/16);
AFp8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing AF contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, AF7, AFz, AF8, feedback);
AF7h  = elec1020_fraction(cnt1, cnt2,  1/16);
AF5   = elec1020_fraction(cnt1, cnt2,  2/16);
AF5h  = elec1020_fraction(cnt1, cnt2,  3/16);
AF3   = elec1020_fraction(cnt1, cnt2,  4/16);
AF3h  = elec1020_fraction(cnt1, cnt2,  5/16);
AF1   = elec1020_fraction(cnt1, cnt2,  6/16);
AF1h  = elec1020_fraction(cnt1, cnt2,  7/16);
AF2h  = elec1020_fraction(cnt1, cnt2,  9/16);
AF2   = elec1020_fraction(cnt1, cnt2, 10/16);
AF4h  = elec1020_fraction(cnt1, cnt2, 11/16);
AF4   = elec1020_fraction(cnt1, cnt2, 12/16);
AF6h  = elec1020_fraction(cnt1, cnt2, 13/16);
AF6   = elec1020_fraction(cnt1, cnt2, 14/16);
AF8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing AFF contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, AFF7, AFFz, AFF8, feedback);
AFF7h  = elec1020_fraction(cnt1, cnt2,  1/16);
AFF5   = elec1020_fraction(cnt1, cnt2,  2/16);
AFF5h  = elec1020_fraction(cnt1, cnt2,  3/16);
AFF3   = elec1020_fraction(cnt1, cnt2,  4/16);
AFF3h  = elec1020_fraction(cnt1, cnt2,  5/16);
AFF1   = elec1020_fraction(cnt1, cnt2,  6/16);
AFF1h  = elec1020_fraction(cnt1, cnt2,  7/16);
AFF2h  = elec1020_fraction(cnt1, cnt2,  9/16);
AFF2   = elec1020_fraction(cnt1, cnt2, 10/16);
AFF4h  = elec1020_fraction(cnt1, cnt2, 11/16);
AFF4   = elec1020_fraction(cnt1, cnt2, 12/16);
AFF6h  = elec1020_fraction(cnt1, cnt2, 13/16);
AFF6   = elec1020_fraction(cnt1, cnt2, 14/16);
AFF8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing F contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, F7, Fz, F8, feedback);
F7h  = elec1020_fraction(cnt1, cnt2,  1/16);
F5   = elec1020_fraction(cnt1, cnt2,  2/16);
F5h  = elec1020_fraction(cnt1, cnt2,  3/16);
F3   = elec1020_fraction(cnt1, cnt2,  4/16);
F3h  = elec1020_fraction(cnt1, cnt2,  5/16);
F1   = elec1020_fraction(cnt1, cnt2,  6/16);
F1h  = elec1020_fraction(cnt1, cnt2,  7/16);
F2h  = elec1020_fraction(cnt1, cnt2,  9/16);
F2   = elec1020_fraction(cnt1, cnt2, 10/16);
F4h  = elec1020_fraction(cnt1, cnt2, 11/16);
F4   = elec1020_fraction(cnt1, cnt2, 12/16);
F6h  = elec1020_fraction(cnt1, cnt2, 13/16);
F6   = elec1020_fraction(cnt1, cnt2, 14/16);
F8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing FFC contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, FFT7, FFCz, FFT8, feedback);
FFT7h  = elec1020_fraction(cnt1, cnt2,  1/16);
FFC5   = elec1020_fraction(cnt1, cnt2,  2/16);
FFC5h  = elec1020_fraction(cnt1, cnt2,  3/16);
FFC3   = elec1020_fraction(cnt1, cnt2,  4/16);
FFC3h  = elec1020_fraction(cnt1, cnt2,  5/16);
FFC1   = elec1020_fraction(cnt1, cnt2,  6/16);
FFC1h  = elec1020_fraction(cnt1, cnt2,  7/16);
FFC2h  = elec1020_fraction(cnt1, cnt2,  9/16);
FFC2   = elec1020_fraction(cnt1, cnt2, 10/16);
FFC4h  = elec1020_fraction(cnt1, cnt2, 11/16);
FFC4   = elec1020_fraction(cnt1, cnt2, 12/16);
FFC6h  = elec1020_fraction(cnt1, cnt2, 13/16);
FFC6   = elec1020_fraction(cnt1, cnt2, 14/16);
FFT8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing FC contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, FT7, FCz, FT8, feedback);
FT7h  = elec1020_fraction(cnt1, cnt2,  1/16);
FC5   = elec1020_fraction(cnt1, cnt2,  2/16);
FC5h  = elec1020_fraction(cnt1, cnt2,  3/16);
FC3   = elec1020_fraction(cnt1, cnt2,  4/16);
FC3h  = elec1020_fraction(cnt1, cnt2,  5/16);
FC1   = elec1020_fraction(cnt1, cnt2,  6/16);
FC1h  = elec1020_fraction(cnt1, cnt2,  7/16);
FC2h  = elec1020_fraction(cnt1, cnt2,  9/16);
FC2   = elec1020_fraction(cnt1, cnt2, 10/16);
FC4h  = elec1020_fraction(cnt1, cnt2, 11/16);
FC4   = elec1020_fraction(cnt1, cnt2, 12/16);
FC6h  = elec1020_fraction(cnt1, cnt2, 13/16);
FC6   = elec1020_fraction(cnt1, cnt2, 14/16);
FT8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing FCC contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, FTT7, FCCz, FTT8, feedback);
FTT7h  = elec1020_fraction(cnt1, cnt2,  1/16);
FCC5   = elec1020_fraction(cnt1, cnt2,  2/16);
FCC5h  = elec1020_fraction(cnt1, cnt2,  3/16);
FCC3   = elec1020_fraction(cnt1, cnt2,  4/16);
FCC3h  = elec1020_fraction(cnt1, cnt2,  5/16);
FCC1   = elec1020_fraction(cnt1, cnt2,  6/16);
FCC1h  = elec1020_fraction(cnt1, cnt2,  7/16);
FCC2h  = elec1020_fraction(cnt1, cnt2,  9/16);
FCC2   = elec1020_fraction(cnt1, cnt2, 10/16);
FCC4h  = elec1020_fraction(cnt1, cnt2, 11/16);
FCC4   = elec1020_fraction(cnt1, cnt2, 12/16);
FCC6h  = elec1020_fraction(cnt1, cnt2, 13/16);
FCC6   = elec1020_fraction(cnt1, cnt2, 14/16);
FTT8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing CCP contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, TTP7, CCPz, TTP8, feedback);
TTP7h  = elec1020_fraction(cnt1, cnt2,  1/16);
CCP5   = elec1020_fraction(cnt1, cnt2,  2/16);
CCP5h  = elec1020_fraction(cnt1, cnt2,  3/16);
CCP3   = elec1020_fraction(cnt1, cnt2,  4/16);
CCP3h  = elec1020_fraction(cnt1, cnt2,  5/16);
CCP1   = elec1020_fraction(cnt1, cnt2,  6/16);
CCP1h  = elec1020_fraction(cnt1, cnt2,  7/16);
CCP2h  = elec1020_fraction(cnt1, cnt2,  9/16);
CCP2   = elec1020_fraction(cnt1, cnt2, 10/16);
CCP4h  = elec1020_fraction(cnt1, cnt2, 11/16);
CCP4   = elec1020_fraction(cnt1, cnt2, 12/16);
CCP6h  = elec1020_fraction(cnt1, cnt2, 13/16);
CCP6   = elec1020_fraction(cnt1, cnt2, 14/16);
TTP8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing CP contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, TP7, CPz, TP8, feedback);
TP7h  = elec1020_fraction(cnt1, cnt2,  1/16);
CP5   = elec1020_fraction(cnt1, cnt2,  2/16);
CP5h  = elec1020_fraction(cnt1, cnt2,  3/16);
CP3   = elec1020_fraction(cnt1, cnt2,  4/16);
CP3h  = elec1020_fraction(cnt1, cnt2,  5/16);
CP1   = elec1020_fraction(cnt1, cnt2,  6/16);
CP1h  = elec1020_fraction(cnt1, cnt2,  7/16);
CP2h  = elec1020_fraction(cnt1, cnt2,  9/16);
CP2   = elec1020_fraction(cnt1, cnt2, 10/16);
CP4h  = elec1020_fraction(cnt1, cnt2, 11/16);
CP4   = elec1020_fraction(cnt1, cnt2, 12/16);
CP6h  = elec1020_fraction(cnt1, cnt2, 13/16);
CP6   = elec1020_fraction(cnt1, cnt2, 14/16);
TP8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing CPP contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, TPP7, CPPz, TPP8, feedback);
TPP7h  = elec1020_fraction(cnt1, cnt2,  1/16);
CPP5   = elec1020_fraction(cnt1, cnt2,  2/16);
CPP5h  = elec1020_fraction(cnt1, cnt2,  3/16);
CPP3   = elec1020_fraction(cnt1, cnt2,  4/16);
CPP3h  = elec1020_fraction(cnt1, cnt2,  5/16);
CPP1   = elec1020_fraction(cnt1, cnt2,  6/16);
CPP1h  = elec1020_fraction(cnt1, cnt2,  7/16);
CPP2h  = elec1020_fraction(cnt1, cnt2,  9/16);
CPP2   = elec1020_fraction(cnt1, cnt2, 10/16);
CPP4h  = elec1020_fraction(cnt1, cnt2, 11/16);
CPP4   = elec1020_fraction(cnt1, cnt2, 12/16);
CPP6h  = elec1020_fraction(cnt1, cnt2, 13/16);
CPP6   = elec1020_fraction(cnt1, cnt2, 14/16);
TPP8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing P contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, P7, Pz, P8, feedback);
P7h  = elec1020_fraction(cnt1, cnt2,  1/16);
P5   = elec1020_fraction(cnt1, cnt2,  2/16);
P5h  = elec1020_fraction(cnt1, cnt2,  3/16);
P3   = elec1020_fraction(cnt1, cnt2,  4/16);
P3h  = elec1020_fraction(cnt1, cnt2,  5/16);
P1   = elec1020_fraction(cnt1, cnt2,  6/16);
P1h  = elec1020_fraction(cnt1, cnt2,  7/16);
P2h  = elec1020_fraction(cnt1, cnt2,  9/16);
P2   = elec1020_fraction(cnt1, cnt2, 10/16);
P4h  = elec1020_fraction(cnt1, cnt2, 11/16);
P4   = elec1020_fraction(cnt1, cnt2, 12/16);
P6h  = elec1020_fraction(cnt1, cnt2, 13/16);
P6   = elec1020_fraction(cnt1, cnt2, 14/16);
P8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing PPO contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, PPO7, PPOz, PPO8, feedback);
PPO7h  = elec1020_fraction(cnt1, cnt2,  1/16);
PPO5   = elec1020_fraction(cnt1, cnt2,  2/16);
PPO5h  = elec1020_fraction(cnt1, cnt2,  3/16);
PPO3   = elec1020_fraction(cnt1, cnt2,  4/16);
PPO3h  = elec1020_fraction(cnt1, cnt2,  5/16);
PPO1   = elec1020_fraction(cnt1, cnt2,  6/16);
PPO1h  = elec1020_fraction(cnt1, cnt2,  7/16);
PPO2h  = elec1020_fraction(cnt1, cnt2,  9/16);
PPO2   = elec1020_fraction(cnt1, cnt2, 10/16);
PPO4h  = elec1020_fraction(cnt1, cnt2, 11/16);
PPO4   = elec1020_fraction(cnt1, cnt2, 12/16);
PPO6h  = elec1020_fraction(cnt1, cnt2, 13/16);
PPO6   = elec1020_fraction(cnt1, cnt2, 14/16);
PPO8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing PO contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, PO7, POz, PO8, feedback);
PO7h  = elec1020_fraction(cnt1, cnt2,  1/16);
PO5   = elec1020_fraction(cnt1, cnt2,  2/16);
PO5h  = elec1020_fraction(cnt1, cnt2,  3/16);
PO3   = elec1020_fraction(cnt1, cnt2,  4/16);
PO3h  = elec1020_fraction(cnt1, cnt2,  5/16);
PO1   = elec1020_fraction(cnt1, cnt2,  6/16);
PO1h  = elec1020_fraction(cnt1, cnt2,  7/16);
PO2h  = elec1020_fraction(cnt1, cnt2,  9/16);
PO2   = elec1020_fraction(cnt1, cnt2, 10/16);
PO4h  = elec1020_fraction(cnt1, cnt2, 11/16);
PO4   = elec1020_fraction(cnt1, cnt2, 12/16);
PO6h  = elec1020_fraction(cnt1, cnt2, 13/16);
PO6   = elec1020_fraction(cnt1, cnt2, 14/16);
PO8h  = elec1020_fraction(cnt1, cnt2, 15/16);

fprintf('constructing POO contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, POO7, POOz, POO8, feedback);
POO7h  = elec1020_fraction(cnt1, cnt2,  1/16);
POO5   = elec1020_fraction(cnt1, cnt2,  2/16);
POO5h  = elec1020_fraction(cnt1, cnt2,  3/16);
POO3   = elec1020_fraction(cnt1, cnt2,  4/16);
POO3h  = elec1020_fraction(cnt1, cnt2,  5/16);
POO1   = elec1020_fraction(cnt1, cnt2,  6/16);
POO1h  = elec1020_fraction(cnt1, cnt2,  7/16);
POO2h  = elec1020_fraction(cnt1, cnt2,  9/16);
POO2   = elec1020_fraction(cnt1, cnt2, 10/16);
POO4h  = elec1020_fraction(cnt1, cnt2, 11/16);
POO4   = elec1020_fraction(cnt1, cnt2, 12/16);
POO6h  = elec1020_fraction(cnt1, cnt2, 13/16);
POO6   = elec1020_fraction(cnt1, cnt2, 14/16);
POO8h  = elec1020_fraction(cnt1, cnt2, 15/16);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start contouring the low electrode locations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% low horizontal ant-post through T9
fprintf('constructing low horizontal left contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, Nz, T9, Iz, feedback);
AFp9 = elec1020_fraction(cnt1, cnt2,  3/20);
AF9  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF9 = elec1020_fraction(cnt1, cnt2,  5/20);
F9   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT9 = elec1020_fraction(cnt1, cnt2,  7/20);
FT9  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT9 = elec1020_fraction(cnt1, cnt2,  9/20);
T9   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP9 = elec1020_fraction(cnt1, cnt2, 11/20);
TP9  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP9 = elec1020_fraction(cnt1, cnt2, 13/20);
P9   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO9 = elec1020_fraction(cnt1, cnt2, 15/20);
PO9  = elec1020_fraction(cnt1, cnt2, 16/20);
POO9 = elec1020_fraction(cnt1, cnt2, 17/20);
I1   = elec1020_fraction(cnt1, cnt2, 18/20);
I1h  = elec1020_fraction(cnt1, cnt2, 19/20);

[cnt1, cnt2] = elec1020_follow(pnt, dhk, NFpz, T9h, OIz, feedback);
AFp9h = elec1020_fraction(cnt1, cnt2,  3/20);
AF9h  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF9h = elec1020_fraction(cnt1, cnt2,  5/20);
F9h   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT9h = elec1020_fraction(cnt1, cnt2,  7/20);
FT9h  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT9h = elec1020_fraction(cnt1, cnt2,  9/20);
T9h   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP9h = elec1020_fraction(cnt1, cnt2, 11/20);
TP9h  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP9h = elec1020_fraction(cnt1, cnt2, 13/20);
P9h   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO9h = elec1020_fraction(cnt1, cnt2, 15/20);
PO9h  = elec1020_fraction(cnt1, cnt2, 16/20);
POO9h = elec1020_fraction(cnt1, cnt2, 17/20);
OI1   = elec1020_fraction(cnt1, cnt2, 18/20);
OI1h  = elec1020_fraction(cnt1, cnt2, 19/20);

% low horizontal ant-post through T10
fprintf('constructing low horizontal right contour\n');
[cnt1, cnt2] = elec1020_follow(pnt, dhk, Nz, T10, Iz, feedback);
AFp10 = elec1020_fraction(cnt1, cnt2,  3/20);
AF10  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF10 = elec1020_fraction(cnt1, cnt2,  5/20);
F10   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT10 = elec1020_fraction(cnt1, cnt2,  7/20);
FT10  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT10 = elec1020_fraction(cnt1, cnt2,  9/20);
T10   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP10 = elec1020_fraction(cnt1, cnt2, 11/20);
TP10  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP10 = elec1020_fraction(cnt1, cnt2, 13/20);
P10   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO10 = elec1020_fraction(cnt1, cnt2, 15/20);
PO10  = elec1020_fraction(cnt1, cnt2, 16/20);
POO10 = elec1020_fraction(cnt1, cnt2, 17/20);
I2    = elec1020_fraction(cnt1, cnt2, 18/20);
I2h   = elec1020_fraction(cnt1, cnt2, 19/20);

[cnt1, cnt2] = elec1020_follow(pnt, dhk, NFpz, T10h, OIz, feedback);
AFp10h = elec1020_fraction(cnt1, cnt2,  3/20);
AF10h  = elec1020_fraction(cnt1, cnt2,  4/20);
AFF10h = elec1020_fraction(cnt1, cnt2,  5/20);
F10h   = elec1020_fraction(cnt1, cnt2,  6/20);
FFT10h = elec1020_fraction(cnt1, cnt2,  7/20);
FT10h  = elec1020_fraction(cnt1, cnt2,  8/20);
FTT10h = elec1020_fraction(cnt1, cnt2,  9/20);
T10h   = elec1020_fraction(cnt1, cnt2, 10/20);
TTP10h = elec1020_fraction(cnt1, cnt2, 11/20);
TP10h  = elec1020_fraction(cnt1, cnt2, 12/20);
TPP10h = elec1020_fraction(cnt1, cnt2, 13/20);
P10h   = elec1020_fraction(cnt1, cnt2, 14/20);
PPO10h = elec1020_fraction(cnt1, cnt2, 15/20);
PO10h  = elec1020_fraction(cnt1, cnt2, 16/20);
POO10h = elec1020_fraction(cnt1, cnt2, 17/20);
OI2    = elec1020_fraction(cnt1, cnt2, 18/20);
OI2h   = elec1020_fraction(cnt1, cnt2, 19/20);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect all the computed electrode locations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lab = {
  'LPA'
  'RPA'
  'NAS'
  'INI'
  'Nz'
  'Fp1'
  'Fpz'
  'Fp2'
  'AF9'
  'AF7'
  'AF5'
  'AF3'
  'AF1'
  'AFz'
  'AF2'
  'AF4'
  'AF6'
  'AF8'
  'AF10'
  'F9'
  'F7'
  'F5'
  'F3'
  'F1'
  'Fz'
  'F2'
  'F4'
  'F6'
  'F8'
  'F10'
  'FT9'
  'FT7'
  'FC5'
  'FC3'
  'FC1'
  'FCz'
  'FC2'
  'FC4'
  'FC6'
  'FT8'
  'FT10'
  'T9'
  'T7'
  'C5'
  'C3'
  'C1'
  'Cz'
  'C2'
  'C4'
  'C6'
  'T8'
  'T10'
  'TP9'
  'TP7'
  'CP5'
  'CP3'
  'CP1'
  'CPz'
  'CP2'
  'CP4'
  'CP6'
  'TP8'
  'TP10'
  'P9'
  'P7'
  'P5'
  'P3'
  'P1'
  'Pz'
  'P2'
  'P4'
  'P6'
  'P8'
  'P10'
  'PO9'
  'PO7'
  'PO5'
  'PO3'
  'PO1'
  'POz'
  'PO2'
  'PO4'
  'PO6'
  'PO8'
  'PO10'
  'O1'
  'Oz'
  'O2'
  'I1'
  'Iz'
  'I2'
  'AFp9h'
  'AFp7h'
  'AFp5h'
  'AFp3h'
  'AFp1h'
  'AFp2h'
  'AFp4h'
  'AFp6h'
  'AFp8h'
  'AFp10h'
  'AFF9h'
  'AFF7h'
  'AFF5h'
  'AFF3h'
  'AFF1h'
  'AFF2h'
  'AFF4h'
  'AFF6h'
  'AFF8h'
  'AFF10h'
  'FFT9h'
  'FFT7h'
  'FFC5h'
  'FFC3h'
  'FFC1h'
  'FFC2h'
  'FFC4h'
  'FFC6h'
  'FFT8h'
  'FFT10h'
  'FTT9h'
  'FTT7h'
  'FCC5h'
  'FCC3h'
  'FCC1h'
  'FCC2h'
  'FCC4h'
  'FCC6h'
  'FTT8h'
  'FTT10h'
  'TTP9h'
  'TTP7h'
  'CCP5h'
  'CCP3h'
  'CCP1h'
  'CCP2h'
  'CCP4h'
  'CCP6h'
  'TTP8h'
  'TTP10h'
  'TPP9h'
  'TPP7h'
  'CPP5h'
  'CPP3h'
  'CPP1h'
  'CPP2h'
  'CPP4h'
  'CPP6h'
  'TPP8h'
  'TPP10h'
  'PPO9h'
  'PPO7h'
  'PPO5h'
  'PPO3h'
  'PPO1h'
  'PPO2h'
  'PPO4h'
  'PPO6h'
  'PPO8h'
  'PPO10h'
  'POO9h'
  'POO7h'
  'POO5h'
  'POO3h'
  'POO1h'
  'POO2h'
  'POO4h'
  'POO6h'
  'POO8h'
  'POO10h'
  'OI1h'
  'OI2h'
  'Fp1h'
  'Fp2h'
  'AF9h'
  'AF7h'
  'AF5h'
  'AF3h'
  'AF1h'
  'AF2h'
  'AF4h'
  'AF6h'
  'AF8h'
  'AF10h'
  'F9h'
  'F7h'
  'F5h'
  'F3h'
  'F1h'
  'F2h'
  'F4h'
  'F6h'
  'F8h'
  'F10h'
  'FT9h'
  'FT7h'
  'FC5h'
  'FC3h'
  'FC1h'
  'FC2h'
  'FC4h'
  'FC6h'
  'FT8h'
  'FT10h'
  'T9h'
  'T7h'
  'C5h'
  'C3h'
  'C1h'
  'C2h'
  'C4h'
  'C6h'
  'T8h'
  'T10h'
  'TP9h'
  'TP7h'
  'CP5h'
  'CP3h'
  'CP1h'
  'CP2h'
  'CP4h'
  'CP6h'
  'TP8h'
  'TP10h'
  'P9h'
  'P7h'
  'P5h'
  'P3h'
  'P1h'
  'P2h'
  'P4h'
  'P6h'
  'P8h'
  'P10h'
  'PO9h'
  'PO7h'
  'PO5h'
  'PO3h'
  'PO1h'
  'PO2h'
  'PO4h'
  'PO6h'
  'PO8h'
  'PO10h'
  'O1h'
  'O2h'
  'I1h'
  'I2h'
  'AFp9'
  'AFp7'
  'AFp5'
  'AFp3'
  'AFp1'
  'AFpz'
  'AFp2'
  'AFp4'
  'AFp6'
  'AFp8'
  'AFp10'
  'AFF9'
  'AFF7'
  'AFF5'
  'AFF3'
  'AFF1'
  'AFFz'
  'AFF2'
  'AFF4'
  'AFF6'
  'AFF8'
  'AFF10'
  'FFT9'
  'FFT7'
  'FFC5'
  'FFC3'
  'FFC1'
  'FFCz'
  'FFC2'
  'FFC4'
  'FFC6'
  'FFT8'
  'FFT10'
  'FTT9'
  'FTT7'
  'FCC5'
  'FCC3'
  'FCC1'
  'FCCz'
  'FCC2'
  'FCC4'
  'FCC6'
  'FTT8'
  'FTT10'
  'TTP9'
  'TTP7'
  'CCP5'
  'CCP3'
  'CCP1'
  'CCPz'
  'CCP2'
  'CCP4'
  'CCP6'
  'TTP8'
  'TTP10'
  'TPP9'
  'TPP7'
  'CPP5'
  'CPP3'
  'CPP1'
  'CPPz'
  'CPP2'
  'CPP4'
  'CPP6'
  'TPP8'
  'TPP10'
  'PPO9'
  'PPO7'
  'PPO5'
  'PPO3'
  'PPO1'
  'PPOz'
  'PPO2'
  'PPO4'
  'PPO6'
  'PPO8'
  'PPO10'
  'POO9'
  'POO7'
  'POO5'
  'POO3'
  'POO1'
  'POOz'
  'POO2'
  'POO4'
  'POO6'
  'POO8'
  'POO10'
  'OI1'
  'OIz'
  'OI2'
  'T3' % this is now called T7
  'T4' % this is now called T8
  'T5' % this is now called P7
  'T6' % this is now called P8
  'M1' % left mastoid
  'M2' % right mastoid
  'A1' % left ear lobe
  'A2' % right ear lobe
  };

% allow for upper case electrode position labels
NAS = nas;
INI = ini;
LPA = lpa;
RPA = rpa;

% it would be possible to assign locations to these old locations
% but there are no locations determined for M1/2 and A1/2
if false
  T3 = T7;
  T4 = T8;
  T5 = P7;
  T6 = P8;
end

% assign the known local electrode positions to the output matrix
nlab = numel(lab);
elc = ones(nlab,3) * nan;
for i=1:nlab
  if exist(lab{i}, 'var')
    eval(sprintf('elc(%d,:) = %s;', i, lab{i}));
  else
    fprintf('not placing electrode %s\n', lab{i});
  end
end

% remove unknown electrode positions
sel = ~isnan(elc(:,1));
elc = elc(sel, :);
lab = lab(sel);

if feedback
  elec = [];
  elec.elecpos = elc;
  elec.label = lab;
  ft_plot_sens(elec)
end


text:
function mtrExportFiberGroupToMetrotrac(outFile, fgFile, faNiftiFile)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% mtrExportFiberGroupToMetrotrac(outFile, [fgFile], [faNiftiFile])
%%
%% outFile - filename for metrotrac pathways
%%
%%
%% Author: Anthony Sherbondy
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ieNotDefined('fgFile')
    [f,p] = uigetfile({'*.mat';'*.*'},'Select a fiber group file for input...');
    if(isnumeric(f)), disp('Conversion canceled.'); return; end
    fgFile = fullfile(p,f);
end
fg = dtiNewFiberGroup();
fg = dtiReadFibers(fgFile);

if ieNotDefined('faNiftiFile')
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a FA file for reference to DTI space input...');
    if(isnumeric(f)), disp('Conversion canceled.'); return; end
    faNiftiFile = fullfile(p,f);
end
nifa = niftiRead(faNiftiFile);
mmPerVox = nifa.pixdim;
xformToAcpc = nifa.qto_xyz;
scene_dim = size(nifa.data);

% Matrix for converting fiber group coords into the metrotrac pathway space
xformFromAcpc = diag([(mmPerVox(:)') 1])*inv(xformToAcpc);

% HACK: Unchecked and can be bogus because DTIQuery doesn't use this
% Get ACPC from transform
ACPC = [xformToAcpc(1,4)/xformToAcpc(1,1) xformToAcpc(2,4)/xformToAcpc(2,2) xformToAcpc(3,4)/xformToAcpc(3,3)];

% Create blank database for finding some of the sizes to write out
pdb = mtrPathwayDatabase();

numstats = length( pdb.pathway_statistic_headers );

% Redo offset
offset = 5*4 + 6*8 + numstats * getStatisticsHeaderSizeWordAligned(pdb); % 5 uint,6 doubles
%5* sizeof (int) + 6 * sizeof(double) + pathways->getNumPathStatistics()*sizeof (DTIPathwayStatisticHeader);
% bool == char, but not guaranteed on all platforms

fid = fopen(outFile, 'wb');

% Writing header
fwrite(fid,offset,'uint');

fwrite(fid,scene_dim,'uint');
fwrite(fid,mmPerVox,'double');
fwrite(fid,ACPC,'double');
fwrite(fid,numstats,'uint');

% Write Stat Header
saveStatisticsHeader(pdb,fid);

%Writing path info
numpaths = length(fg.fibers);
fwrite(fid,numpaths,'uint');

for ff = 1:numpaths
    path_offset = 3 * 4  + numstats*8; % 3 int, numstats double
    fwrite(fid,path_offset,'int');

    %nn = size( fg.fibers{ff}.coords, 2 );
    nn = size( fg.fibers{ff}, 2 );

    fwrite(fid,nn,'int');
    fwrite(fid,1,'int'); % algo type HACK
    fwrite(fid,1,'int'); % seedpointindex HACK ???

    % Stats
    if( numstats ~= 0 )
        error('Can not handle nonzero statistics with path.');
    end

    % Writing path nodes
    %pos = fg.fibers{ff}.coords;
    pos = fg.fibers{ff};
    pos = mrAnatXformCoords(xformFromAcpc, pos')';
    % Change from 1 based to 0 based positions in mm space
    pos = pos - repmat(mmPerVox(:),1,size(pos,2));
    fwrite(fid,pos,'double');

    % HACK: Don't worry numstats has to be zero right now
% %     % Writing stats values per position
% %     for as = 1:numstats
% %         if( this.pathway_statistic_headers(as).is_computed_per_point )
% %             fwrite(fid,this.pathways(p).point_stat_array(as,:),'double');
% %         end
% %     end

end

disp('Saved Database.');
fclose(fid);

function pos = asMatrixStruct(this)
pos = [this.xpos; this.ypos; this.zpos];
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function timers = cvx_toc

%CVX_TOC Returns the current timing information.
%
%   CVX collects the following timing information:
%   --- Elapsed time since CVX_TIC was last called. If CVX_TIC was never
%       called, this measures the time since CVX was first initialized.
%   --- Elapsed time spent between top-level CVX_BEGIN and CVX_END comamnds
%   --- Elapsed time spent within the top-level CVX_END itself. This
%       includes final problem extraction, presolving, and solving.
%   --- Elapsed time spent calling the numerical solver.
%   To retrieve the current totals for these numbers, type CVX_TOC.
%   CVX_TIC resets these numbers to zero.

global cvx___
if isempty( cvx___ ),
    error( 'CVX has not yet been used (or the global workspace has been cleared).' );
end
tstart = tic;
timers = double(cvx___.timers);
dbltim = isa( timers, 'double' );
if dbltim,
    tstart = double(tstart);
end
timers(1) = tstart - timers(1);
if isempty(cvx___.increment),
    switch computer,
        case { 'MACI', 'MACI64' },
            cvx___.increment = 1e9;
        otherwise,
            t1 = tic; pause(0.25); t2 = tic; t3 = toc(t1); t4=toc(t2);
            if dbltim,
                cvx___.increment = (double(t2)-double(t1))/(double(t3)-0.5*double(t4));
            else
                two = uint64(2);
                cvx___.increment = double( two * ( t2 - t1 ) ) / ( two * t3 - t4 );
            end
    end
end
timers = [ timers(1), timers(1:3) - timers(2:4), timers(4) ];
timers = timers / cvx___.increment;
if nargout == 0,
    fprintf( 'Total time:      %g sec\n', timers(1) );
    fprintf( 'Outside of CVX:  %g sec\n', timers(2) );
    fprintf( 'Model building:  %g sec\n', timers(3) );
    fprintf( 'Presolving:      %g sec\n', timers(4) );
    fprintf( 'Solving:         %g sec\n', timers(5) );
    clear timers
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX