text:
function res = viewHelpParameter(paramIn)
% Maps the paramIn to the help file, displaying an explanation of what each
% viewGet/Set parameter should do.
% This function should never
% be called directly, but is instead wrapped by viewMapParameterField.
%
%    res = viewHelpParameter(paramIn);
%
% Add help functionality for viewGet/Set.
%
% By using this function, we can get help from the program itself when
% trying to call a certain field. This embeds knowledge of what each field
% does into the program, rather than into people's heads.
%
% res returns a multi-line string that, when printed, will display across
% multiple lines, similar to the multi-line comments previously stored in
% the Matlab code.
%
% Examples:
%   viewHelpParameter('name')
%   viewHelpParameter('curdt')

global DictViewHelp;

if isempty(DictViewHelp)
    DictViewHelp = containers.Map;

    DictViewHelp('homedir') = ['Return full path to directory.' char(10) ...
        'homedir = viewGet(vw, ''Home Directory'');'];
    DictViewHelp('sessionname') = ['Retrun name of session, such as ''BW090616-8Bars-14deg''.' char(10) ...
        'sessionName = viewGet(vw, ''session name'');'];
    DictViewHelp('subject') = ['Return name of subject, such as ''Wandell''' char(10) ...
        'subject = viewGet(vw, ''subject'')'];
    DictViewHelp('name') = ['Return name of view, such as ''INPLANE{1}''' char(10) ...
        'name = viewGet(vw, ''view name'');'];
    DictViewHelp('annotation') = ['Return description of currently selected scan (string, such as' char(10) ...
        '''14 Deg 8 Bars with blanks'')' char(10) ...
        'annotation = viewGet(vw, ''annotation'');'];
    DictViewHelp('annotations') = ['Return a cell array with descriptions of each scan in current' char(10) ...
        'dataType' char(10) ...
        'annotations = viewGet(vw, ''annotations'');'];
    DictViewHelp('viewtype') = ['Return the view type (''Gray'', ''Volume'', ''Inplane'', ''Flat'', etc) ' char(10) ...
        'viewType = viewGet(vw, ''View Type'');'];
    DictViewHelp('subdir') = ['Return the sub directory name (not the full path) with data for' char(10) ...
        'current view ' char(10) ...
        'subdir = viewGet(vw, ''sub directory'');'];
    DictViewHelp('curscan') = ['Return the currently selected scan number' char(10) ...
        'curscan = viewget(vw, ''Current Scan'');'];
    DictViewHelp('curslice') = ['Return the current slice number. This is the actual slice number' char(10) ...
        'if we are in the Inplane view. It is the plane number (sag, cor,' char(10) ...
        'or axi) of the currently selected plane if we are in the Volume' char(10) ...
        'view. And it is 1 or 2 in the Flat view (for left or right).' char(10) ...
        'curslice = viewGet(vw, ''Current Slice'');'];
    DictViewHelp('nscans') = ['Return the number of scans in the currently selected dataTYPE' char(10) ...
        'nscans = viewGet(vw, ''Number of Scans'');'];
    DictViewHelp('nslices') = ['Return the number of slices in the current view struct' char(10) ...
        'nslices = viewGet(vw, ''Number of Slices'');'];
    DictViewHelp('montageslices') = ['Return the current subset of slices that are visible in the GUI.' char(10) ...
        'montageSlices = viewGet(vw, ''Montage Slices'');' char(10) ...
        '' char(10) ...
        'only for some view types: inplane montage, flat level' char(10) ...
        'Well, this seems not to be true since it returns a value in the' char(10) ...
        'Volume view.'];
    DictViewHelp('dtname') = ['Return the name of the currently selected dataTYPE' char(10) ...
        'dtName = viewGet(vw, ''Data TYPE Name'');' char(10) ...
        '' char(10) ...
        'ras 06/07: note that I ''stole'' the ''datatype'' alias from' char(10) ...
        'the dtnumber property: I think the name is more relevant'];
    DictViewHelp('curdt') = ['Return the number of the currently selected dataTYPE' char(10) ...
        'dtNum = viewGet(vw, ''Current Data TYPE'');'];
    DictViewHelp('dtstruct') = ['Return the currently selected dataTYPE struct' char(10) ...
        'dtStruct = viewGet(vw, ''DT struct'');       '];
    DictViewHelp('coherence') = ['Coherence for all voxels, all scans in current dataTYPE        ' char(10) ...
        'co = viewGet(vw, ''Coherence'');        '];
    DictViewHelp('scanco') = ['Coherence for single scan' char(10) ...
        'scanco = viewGet(vw, ''scan coherence'', 1);        '];
    DictViewHelp('phase') = ['Phase for all voxels, all scans in current dataTYPE    ' char(10) ...
        'ph = viewGet(vw, ''Phase'');'];
    DictViewHelp('scanph') = ['Phase values for single scan' char(10) ...
        'viewGet(vw,''Scan Phase'',1);'];
    DictViewHelp('amplitude') = ['Amplitude for all voxels, all scans in current dataTYPE    ' char(10) ...
        'amp = viewGet(vw, ''Amplitude'');'];
    DictViewHelp('scanamp') = ['Amplitude values for single scan (selected scan or specified' char(10) ...
        'scan).' char(10) ...
        'scan = 1; scanAmp = viewGet(vw,''scan Amp'', scan);' char(10) ...
        'scanAmp = viewGet(vw,''scan Amp'');'];
    DictViewHelp('refph') = ['Return the reference phase used for computing phase-referred' char(10) ...
        'coherence. Should be [0 2*pi]?' char(10) ...
        'refph = viewGet(vw,''reference phase'');'];
    DictViewHelp('ampmap') = ['Return the colormap currently used to display amplutitude data.' char(10) ...
        'Should be 3 x numColors.' char(10) ...
        'ampMap = viewGet(vw, ''amplitude color map'');'];
    DictViewHelp('coherencemap') = ['Return the colormap currently used to display coherence data.' char(10) ...
        'Should be 3 x numColors.' char(10) ...
        'cohMap = viewGet(vw, ''coherence color map'');'];
    DictViewHelp('correlationmap') = ['Return the colormap currently used to display correlation data.' char(10) ...
        'Should be 3 x numColors. ' char(10) ...
        'corMap = viewGet(vw, ''correlation color map'');' char(10) ...
        '' char(10) ...
        '[Q: what is a correlation map and how does it differ from' char(10) ...
        'coherence map?]'];
    DictViewHelp('cothresh') = ['Return the coherence threshold. Should be in [0 1].' char(10) ...
        'cothresh = viewGet(vw, ''Coherence Threshold'');'];
    DictViewHelp('phwin') = ['Return  phWindow values from phWindow sliders (non-hidden views)' char(10) ...
        'or from the view.settings.phWin field (hidden views). If can''t' char(10) ...
        'find either, defaults to [0 2*pi].' char(10) ...
        'phwin = viewGet(vw, ''phase window'');'];
    DictViewHelp('twparams') = ['Return travelling wave parameters.' char(10) ...
        'twparams = viewGet(vw, ''Travelling Wave Parameters'');        '];
    DictViewHelp('cmap') = ['Return the colormap for whichever data view (co, ph, amp, map) is' char(10) ...
        'currently selected.' char(10) ...
        'cmap = viewGet(vw, ''color map'');'];
    DictViewHelp('cmapcolor') = ['Return color portion of current color overlay map.' char(10) ...
        'cmapColor = viewGet(vw, ''cmap color'');'];
    DictViewHelp('cmapgrayscale') = ['Return grayscale portion of current color overlay map' char(10) ...
        'cmapGray = viewGet(vw, ''cmap grayscale'');'];
    DictViewHelp('map') = ['Return the parameter map for the current data type. Map is cell' char(10) ...
        'array 1 x nscans.' char(10) ...
        'map = viewGet(vw, ''map'');'];
    DictViewHelp('mapwin') = ['Return mapWindow values from mapWindow sliders (non-hidden views)' char(10) ...
        'or from the view.settings.mapWin field (hidden views).        ' char(10) ...
        'mapWin = viewGet(vw, ''Map Window'');'];
    DictViewHelp('mapname') = ['Return the name of the current paramter map (string), e.g.,' char(10) ...
        '''eccentricty''. ' char(10) ...
        'mapName = viewGet(vw, ''Map Name'');'];
    DictViewHelp('mapunits') = ['Return the map units for the current paramter map (string), e.g.,' char(10) ...
        '''degrees''. The map units are for display only; they are not used' char(10) ...
        'for calculations.' char(10) ...
        'mapUnits = viewGet(vw, ''Map Units'');'];
    DictViewHelp('mapclip') = ['Return map clip values. These are the clip values for the' char(10) ...
        'colorbar. They are not the clip values in the slider (which are' char(10) ...
        'called ''mapwin''). Values outside of mapclip are colored the same' char(10) ...
        'as the minimum or maximum value according to the color lookup' char(10) ...
        'table. Values outside of mapwin are not shown at all.' char(10) ...
        'mapClip = viewGet(vw, ''Map Window'');'];
    DictViewHelp('scanmap') = ['Return the parameter map for the currently selected or the' char(10) ...
        'specified scan.' char(10) ...
        'scanMap = viewGet(vw, ''scan map'')' char(10) ...
        'scan = 1; scanMap = viewGet(vw, ''scan map'', scan);'];
    DictViewHelp('anatomy') = ['Return the anatomical underlay image.'];
    DictViewHelp('anatomymap') = ['Return the colormap for the anatomical underlay image.' char(10) ...
        'anataomyMap = viewGet(vw, ''Anatomy Map'');'];
    DictViewHelp('anatomynifti') = ['Return the actual nifti struct stored in anat'];
    DictViewHelp('anatclip') = ['Return anatomy clipping values from anatMin and anatMax sliders.' char(10) ...
        'anataomyClip = viewGet(vw, ''Anatomy Clip'');'];
    DictViewHelp('anatslicedim') = ['Returns the dimension of the matrix that is associated with' char(10) ...
        'slices'];
    DictViewHelp('anatslicedims') = ['Returns the dimensions of each 2D array making up each slice'];
    DictViewHelp('anatsize') = ['Load an anatomy if one does not already exist'];
    DictViewHelp('anatomycurrentslice') = ['Return the anatomical underlay image for only one slice' char(10) ...
        'anat = viewGet(vw, ''Anatomy Current Slice'', curSlice);'];
    DictViewHelp('anatsizexyz') = ['The class data store the planes in a different order from the' char(10) ...
        'vw.anat.  If you want sizes that work for the class data, call' char(10) ...
        'this size.' char(10) ...
        'anatSizeXYZ = viewGet(vw, ''Anatomy Size For Class'');'];
    DictViewHelp('brightness') = ['Return the value of the brightness slider bar for the anatomical' char(10) ...
        'underlay image.' char(10) ...
        'brightness = viewGet(vw, ''brightness'');'];
    DictViewHelp('contrast') = ['Return the value of the contrast slider bar for the anatomical' char(10) ...
        'underlay image.' char(10) ...
        'contrast = viewGet(vw, ''contrast'');'];
    DictViewHelp('mmpervox') = ['Return the size of a voxel in mm' char(10) ...
        'mmPerVox = viewGet(vw, ''mm per voxel'');'];
    DictViewHelp('mmpervolvox') = ['Return the size of a gray/volume voxel in mm. ' char(10) ...
        'If we are in the gray or volume view, this is no different from' char(10) ...
        'viewGet(vw, ''mm per vox''). If we are in a different viewType' char(10) ...
        '(inplane, flat),  we need to initialize a hidden gray view' char(10) ...
        'and then read voxel size.' char(10) ...
        'mmpervolvox = viewGet(vw, ''mm per vol vox'');        '];
    DictViewHelp('gray') = '';
        DictViewHelp('gray') = '';
        DictViewHelp('ngraylayers') = ['Return the number of gray layers in the classification file.' char(10) ...
        'Assumes a Gray or Volume view.' char(10) ...
        'nGrayLayers = viewGet(vw, ''Number Gray Layers'');'];
        DictViewHelp('gray') = '';
        DictViewHelp('gray') = '';
        DictViewHelp('scannerxform') = ['Return the transform matrix from INPLANE coordinates (x,y,z, indices) to' char(10) ...
        'scanner coordinates (in mm). Inplane dicoms are needed to get' char(10) ...
        'this information.' char(10) ...
        'The transform will yield: scannerCoords = scannerXform*ipCoords;' char(10) ...
        'scannerxform = viewGet(vw, ''Scanner transform'');' char(10) ...
        'scannerxform = viewGet(vw, ''scannerXform'', [rawFile]);'];
        DictViewHelp('b0dir') = ['Return the direction of the B0 field from a scan as a unit vector' char(10) ...
        'b0vector = viewGet(vw, ''b0 direction'');'];
        DictViewHelp('b0angle') = ['Return the direction of b0 field from a scan as an angle in degrees' char(10) ...
        'b0angle = viewGet(vw, ''b0 angle'');'];
        DictViewHelp('rois') = ['Return ROIs as struct. Includes all ROIs and all fields' char(10) ...
        'ROIs = viewGet(vw, ''ROIs'');'];
        DictViewHelp('roistruct') = ['Return selected or specified ROI as struct. ' char(10) ...
        'ROI = viewGet(vw, ''ROI struct'');' char(10) ...
        'ROI = viewGet(vw, ''ROI struct'', 1);'];
        DictViewHelp('roicoords') = ['Return the coordinates of the currently selected or the specified' char(10) ...
        'ROI' char(10) ...
        'roiCoords = viewGet(vw, ''ROI coords'');' char(10) ...
        'ROI = 1; roiCoords = viewGet(vw, ''ROI coords'', ROI);'];
        DictViewHelp('roiindices') = ['Return the indices of all voxels in the currently selected or the' char(10) ...
        'specified ROI. Only implemented for gray / vol view. Could be' char(10) ...
        'implemented for other views. ' char(10) ...
        'roiIndices = viewGet(vw, ''ROI indices'');' char(10) ...
        'ROI = 1; roiIndices = viewGet(vw, ''ROI indices'', ROI);    '];
        DictViewHelp('gray') = '';
        DictViewHelp('gray') = '';
        DictViewHelp('roivertinds') = ['Return mesh indices for each voxel in a gray view ROI' char(10) ...
        'roiVertexIndices = ...' char(10) ...
        'viewGet(vw, ''ROI Vertex Indices'', [ROI],[mrmPrefs]);        '];
        DictViewHelp('roivertinds') = '';
        DictViewHelp('roivertinds') = '';
        DictViewHelp('roivertinds') = '';
        DictViewHelp('roivertinds') = '';
        DictViewHelp('roiname') = ['Return the name of the currently selected or the requested ROI.' char(10) ...
        'roiName = viewGet(vw, ''ROI name'');' char(10) ...
        'roi = 1; roiName = viewGet(vw, ''ROI name'', roi);'];
        DictViewHelp('roimodified') = ['Return the modification date of the currently selected or the' char(10) ...
        'requested ROI.' char(10) ...
        'roiName = viewGet(vw, ''ROI modified'');' char(10) ...
        'roi = 1; roiName = viewGet(vw, ''ROI modified'', roi);'];
        DictViewHelp('allroinames') = ['Return a cell array of the names of all currently loaded ROIs.' char(10) ...
        'roiNames = viewGet(vw, ''All ROI names'');'];
        DictViewHelp('nrois') = ['Return the number of currently loaded ROIs.' char(10) ...
        'nrois = viewGet(vw, ''n ROIs'');'];
        DictViewHelp('selectedroi') = ['Return the number of the currently selected ROI.' char(10) ...
        'selectedROI = viewGet(vw, ''Selected ROI'');'];
        DictViewHelp('filledperimeter') = ['Return the value of the field "filledPerimeter". Not sure exactly' char(10) ...
        'how this works, but we think this field is 1 if ROIs are drawn as' char(10) ...
        'filled perimeter and 0 if not, but perhaps it can take other' char(10) ...
        'values?' char(10) ...
        'filledPerimeter = viewGet(vw, ''filled perimeter'');'];
        DictViewHelp('selroicolor') = ['Return the color of the currently selected or the requested ROI.' char(10) ...
        'This can be a Matlab character for a color (''c'', ''w'', ''b'', etc)' char(10) ...
        'or an RGB triplet.' char(10) ...
        'selcol = viewGet(vw, ''Selected ROI color'');' char(10) ...
        'roi = 1; selcol = viewGet(vw, ''Selected ROI color'', roi);'];
        DictViewHelp('prevcoords') = ['Return the coords of the previous ROI state. ' char(10) ...
        'But what does this mean? How does it work?' char(10) ...
        'prevCoords = viewGet(vw, ''previous coords'');'];
        DictViewHelp('roistodisplay') = ['Return the number of each ROI that has been selected to be' char(10) ...
        'displayed. This can be one ROI only (the currently selected ROI),' char(10) ...
        'all ROIs, or a subset. A vector of index values is returned.' char(10) ...
        'See roiSetOptions for an explanation of this format.' char(10) ...
        'roisToDisplay = viewGet(vw, ''ROIs To Display'');'];
        DictViewHelp('roidrawmethod') = ['Return the method for drawing ROIs, such as ''perimeter'',' char(10) ...
        '''filledPerimeter''.' char(10) ...
        'ROIdrawMethod = viewGet(vw, ''ROI Draw Method'');'];
        DictViewHelp('showrois') = ['Return the boolean to indicate whether ROIs are drawn or hidden.' char(10) ...
        'Applies to views such as Volume, Gray, Inplane, Flat, and also' char(10) ...
        'meshes. ' char(10) ...
        '[OK, not quite a Boolean because I see that values can be -2, for' char(10) ...
        'example. What are all the possible values and what do they mean?] ' char(10) ...
        'showROIs = viewGet(vw, ''Show ROIs'');        '];
        DictViewHelp('hidevolumerois') = ['Return the Boolean ''hideVolumeROIs''. If true, then ROIs in' char(10) ...
        'Gray/Volume view are not drawn (but ROIs on meshes may still be' char(10) ...
        'drawn.) The reason we might want to do this is because redrawing' char(10) ...
        'ROIs in the Volume view can be very slow. This field works' char(10) ...
        'indepently from the field ''showROIs'', which controls whether ROIs' char(10) ...
        'are displayed in any view and on meshes. So if showROIs is true' char(10) ...
        'and hideVolumeROIs is false, ROIs will show on a mesh, an inplane' char(10) ...
        'view, or a flat view, but not a Gray/Volume view. If showROIs is' char(10) ...
        'false, no ROIs will be shown anywhere.' char(10) ...
        'hidevolumerois = viewGet(vw, ''Hide Volume ROIs'');'];
        DictViewHelp('maskrois') = ['Return the Boolean ''maskROIs''. If maskROIs is true, then maps' char(10) ...
        'that are displayed on meshes are masked out everywhere outside of' char(10) ...
        'the currently displayed ROIs. The currently displayed ROIs might' char(10) ...
        'be just the currently selected ROI, all the ROIs, or a subset of' char(10) ...
        'ROIs. (See the case ''roisToDisplay''.)' char(10) ...
        'maskROIs = viewGet(vw, ''mask ROIs'');'];
        DictViewHelp('tseriesdir') = ['Return tSeries directory for a view; make it if it does not' char(10) ...
        'exist.' char(10) ...
        'makeIt = 0; tDir = viewGet(vw,''tSeriesDir'',makeIt);' char(10) ...
        'makeIt = 1; tDir = viewGet(vw,''tSeriesDir'',makeIt);' char(10) ...
        'tDir = viewGet(vw,''tSeriesDir'')'];
        DictViewHelp('datasize') = ['Return the size of the data arrays, i.e., size of co for' char(10) ...
        'a single scan.' char(10) ...
        'dataSize = viewGet(vw, ''Data Size'');' char(10) ...
        'scan = 1; dataSize = viewGet(vw, ''Data Size'', scan);'];
        DictViewHelp('dim') = ['Return the dimension of data in current slice or specificed slice' char(10) ...
        'dim = viewGet(vw, ''Slice Dimension'')' char(10) ...
        'scan = 1; dim = viewGet(vw, ''Slice Dimension'', scan)        '];
        DictViewHelp('functionalslicedim') = ['Return the dimension of functional data in current slice or' char(10) ...
        'specificed slice' char(10) ...
        'dim = viewGet(vw, ''Slice Dimension'')' char(10) ...
        'scan = 1; dim = viewGet(vw, ''Slice Dimension'', scan)'];
        DictViewHelp('tseries') = ['Return the time series of all data currently loaded into the view' char(10) ...
        'struct.' char(10) ...
        'tseries = viewGet(vw, ''time series'');'];
        DictViewHelp('tseriesslice') = ['Return the time series for the currently selected slice if it is' char(10) ...
        'loaded into the view struct (return blank otherwise).' char(10) ...
        'tseries = viewGet(vw, ''Time Series Slice'');'];
        DictViewHelp('tseriesscan') = ['Return the time series for the current scan (if it is loaded into' char(10) ...
        'the view struct; return blank if it is not loaded).' char(10) ...
        'tseriesScan = viewGet(vw, ''time series scan'');'];
        DictViewHelp('tr') = ['Return the scan TR in seconds' char(10) ...
        'tr = viewGet(vw,''tr'')' char(10) ...
        'scan = 1; tr = viewGet(vw,''tr'',scan)'];
        DictViewHelp('nframes') = ['Return the number of time frames in the current or specified' char(10) ...
        'scan.' char(10) ...
        'nframes = viewGet(vw,''nFrames'');' char(10) ...
        'scan = 1; nframes = viewGet(vw,''nFrames'',scan);'];
        DictViewHelp('ncycles') = ['Return the number of cycles in the current or specified scan' char(10) ...
        '(assuming scan is set up for coranal). ' char(10) ...
        'nycles = viewGet(vw,''ncycles'')' char(10) ...
        'scan = 1; nycles = viewGet(vw,''ncycles'', scan)'];
        DictViewHelp('framestouse') = ['Return a vector of time frames in the current or specified' char(10) ...
        'scan to be used for coranal (block) analyses' char(10) ...
        'frames = viewGet(vw,''frames to use'');' char(10) ...
        'scan = 1; frames = viewGet(vw,''frames to use'',scan);'];
        DictViewHelp('rmfile') = ['Return the path to the currently loaded retinotopy model.' char(10) ...
        'rmFile = viewGet(vw, ''retinotopy model file'');'];
        DictViewHelp('rmmodel') = ['Return the currently loaded retinotopy model struct.' char(10) ...
        'rm = viewGet(vw, ''retinotopy model'');'];
        DictViewHelp('rmcurrent') = ['Return the currently selected retinotopy model struct. Note that' char(10) ...
        'there may be multiple models loaded.' char(10) ...
        'rm = viewGet(vw, ''rm current model'');'];
        DictViewHelp('rmmodelnames') = ['Return the description of currently loaded retinotopy models.' char(10) ...
        'models = viewGet(vw, ''rm model names'');' char(10) ...
        'models = viewGet(vw, ''retinotopy model names'');'];
        DictViewHelp('rmparams') = ['Return the retinotopy model parameters.' char(10) ...
        'params = viewGet(vw, ''Retinotopy Parameters'');'];
        DictViewHelp('rmstimparams') = ['Return the retinotopy model stimulus parameters. This is a subset' char(10) ...
        'of the retinopy model parameters.' char(10) ...
        'stimParams = viewGet(vw, ''RM Stimulus Parameters'');'];
        DictViewHelp('rmmodelnum') = ['Return the retinotopy model number that is currently selected.' char(10) ...
        '(There may be more than one model loaded.)' char(10) ...
        'modelNum = viewGet(vw, ''Retinotopy Model Number'');'];
        DictViewHelp('rmhrf') = ['Return the hrf struct for the current retinopy model. This struct' char(10) ...
        'contains a descriptor (such as ''two gammas (SPM style)'') and the' char(10) ...
        'parameters associated with this function.' char(10) ...
        'rmhrf = viewGet(vw, ''Retinotopy model HRF'');'];
        DictViewHelp('one gamma (boynton style)') = '';
        DictViewHelp('two gammas (spm style)') = '';
        DictViewHelp('impulse') = '';
        DictViewHelp('impulse') = '';
        DictViewHelp('allmeshes') = ['Return the structs for all currently loaded meshes.' char(10) ...
        'allmeshes = viewGet(vw, ''all meshes'');'];
        DictViewHelp('allmeshids') = ['Return the ID list for all meshes. IDs are numbers generated by' char(10) ...
        'mrMesh that are associated with each new mesh session. They are' char(10) ...
        'typically 4 digit numbers starting at 1001. (why??)' char(10) ...
        'idList = viewGet(vw,''All Window IDs'');'];
        DictViewHelp('mesh') = ['Return the mesh structure for the selected or the requested mesh.' char(10) ...
        'If the mesh number is specificied, it indexes the cell array of' char(10) ...
        'meshes currently attached to the view structure. The mesh number' char(10) ...
        'bears no relation to the mesh ID number, which is generated by' char(10) ...
        'the mesh server.' char(10) ...
        'msh = viewGet(vw, ''mesh'');' char(10) ...
        'meshnum = 1; msh = viewGet(vw, ''mesh'', meshnum);'];
        DictViewHelp('currentmesh') = ['Return the mesh structure for the selected mesh. This is' char(10) ...
        'redundant with the case ''mesh''.' char(10) ...
        'msh = viewGet(vw, ''Current Mesh'');'];
        DictViewHelp('meshn') = ['Return the number of the currently selected mesh (index into the' char(10) ...
        'cell array of meshes)' char(10) ...
        'msh = viewGet(vw, ''current mesh number'');'];
        DictViewHelp('meshdata') = ['I think this is supposed to return the data displayed on the' char(10) ...
        'current mesh, but I have just tried it and it doesn''t seem to' char(10) ...
        'work. So what does it do?' char(10) ...
        'meshData = viewGet(vw, ''current mesh data'');'];
        DictViewHelp('nmesh') = ['Return the number of meshes currently attached to the view' char(10) ...
        'struct.' char(10) ...
        'nmesh = viewGet(vw, ''Number of Meshes'');'];
        DictViewHelp('meshnames') = ['Return the name of all meshes currently attached to the view' char(10) ...
        'struct.' char(10) ...
        'meshNames = viewGet(vw, ''mesh names'');'];
        DictViewHelp('meshdir') = ['Return the directory in which the currently selected mesh' char(10) ...
        'resides. Default to anat dir if not found.' char(10) ...
        'meshDir = viewGet(vw, ''mesh directory'');        '];
        DictViewHelp('nodes') = ['Return the array of nodes. Only gray views have nodes. See help' char(10) ...
        'for mrManDist.m for a description of the node structure. In' char(10) ...
        'brief, nodes are 8 x nvoxels. The first 3 rows correspond to the' char(10) ...
        'voxel location and the next 5 correspond to gray graph-related' char(10) ...
        'data.' char(10) ...
        'nodes = viewGet(vw, ''nodes'');'];
        DictViewHelp('xyznodes') = ['Return the xyz coordinates of the gray voxels as found in nodes' char(10) ...
        'array. Assumes a Gray view. See case ''nodes'' and help for' char(10) ...
        'mrManDist for more information.' char(10) ...
        '' char(10) ...
        'Must call this sagittal, axial coronal or whatever the mapping is' char(10) ...
        'ras, 06/07 -- I believe it''s [cor axi sag]. coords is [axi cor sag].' char(10) ...
        '' char(10) ...
        'xyzNodes = viewGet(vw, ''xyz nodes'');'];
        DictViewHelp('nodegraylevel') = ['Return the gray level of each voxel as determined by the nodes' char(10) ...
        'array. Assumes a Gray view. See case ''nodes'' and help for' char(10) ...
        'mrManDist for more information.' char(10) ...
        'nodeGrayLevel = viewGet(vw, ''gray level'');'];
        DictViewHelp('nnodes') = ['Return the number of nodes. Assumes a Gray view. See case ''nodes''' char(10) ...
        'and help for mrManDist for more information. ' char(10) ...
        'nNodes = viewGet(vw, ''number of nodes'');'];
        DictViewHelp('edges') = ['Return the edge structure of the gray graph. Assumes a Gray view.' char(10) ...
        'See help for mrManDist for more information. ' char(10) ...
        'edges = viewGet(vw, ''edges'');'];
        DictViewHelp('nedges') = ['Return the number of edges in the gray graph. Assumes a Gray' char(10) ...
        'view. See case ''edges'' and help for mrManDist for more' char(10) ...
        'information.' char(10) ...
        'nEdges = viewGet(vw, ''number of edges'');      '];
        DictViewHelp('allleftnodes') = ['Return the subset of nodes in the Gray graph that are in the left' char(10) ...
        'hemisphere. See mrgGrowGray and mrManDist.' char(10) ...
        'allLeftNodes = viewGet(vw, ''all left nodes'');'];
        DictViewHelp('allleftedges') = ['Return the subset of edges in the Gray graph that are in the left' char(10) ...
        'hemisphere. See mrgGrowGray and mrManDist.' char(10) ...
        'allLeftEdges = viewGet(vw, ''all left edges'');'];
        DictViewHelp('allrightnodes') = ['Return the subset of nodes in the Gray graph that are in the' char(10) ...
        'right hemisphere. See mrgGrowGray and mrManDist.' char(10) ...
        'allRightNodes = viewGet(vw, ''all right nodes'');'];
        DictViewHelp('allrightedges') = ['Return the subset of edges in the Gray graph that are in the' char(10) ...
        'right hemisphere. See mrgGrowGray and mrManDist.' char(10) ...
        'allRightEdges = viewGet(vw, ''all right edges'');'];
        DictViewHelp('allnodes') = ['Return all nodes from Gray graph by taking union of allLeftNodes' char(10) ...
        'and allRightNodes. ' char(10) ...
        '' char(10) ...
        'This is NOT necessarily the same as simply returning ''vw.nodes''.' char(10) ...
        'When we install a segmentation, we can either keep all the nodes' char(10) ...
        'in the gray graph, or only those that fall within the functional' char(10) ...
        'field of view (to save space). When we do the latter, the fields' char(10) ...
        'vw.coords, vw.nodes, and vw.edges contain only the coords, nodes,' char(10) ...
        'and eges within the functional field of view. However the fields' char(10) ...
        'vw.allLeftNodes, vw.allLeftEdges, vw.allRightNodes, and' char(10) ...
        'vw.allRightEdges contain the edges and nodes for the entire' char(10) ...
        'hemisphere' char(10) ...
        '' char(10) ...
        'Example: nodes = viewGet(vw, ''all nodes'');'];
        DictViewHelp('alledges') = ['Return all edges from Gray graph by taking union of allLeftEdges' char(10) ...
        'and allRightEdges. See ''allnodes'' for explanation.' char(10) ...
        '' char(10) ...
        'Example: edges = viewGet(vw, ''all edges'');'];
        DictViewHelp('coords') = ['Return all the coordinates in the current view. If in Flat view,' char(10) ...
        'return the coordinates for a particular slice (slice specified in' char(10) ...
        'varargin{1}). If in Inplane view, slice specification is' char(10) ...
        'optional. If in Gray or Volume view, slice specification is' char(10) ...
        'ignored.' char(10) ...
        '<gray, volume or inplane> ' char(10) ...
        'coords = viewGet(vw, ''coords'');' char(10) ...
        '<flat or inplane>' char(10) ...
        'slice  = viewGet(vw, ''current slice''); ' char(10) ...
        'coords = viewGet(vw,''coords'', slice);'];
        DictViewHelp('allcoords') = ['Return all coords from Gray graph, including those that are not' char(10) ...
        'included in the functional field of view. See ''allnodes'' for' char(10) ...
        'explanation. If session was initialized with the option ' char(10) ...
        '''keepAllNodes'' == true, then this call will be identical to' char(10) ...
        'viewGet(vw.coords).' char(10) ...
        '' char(10) ...
        'Example: coords = viewGet(vw, ''all coords'');'];
        DictViewHelp('coordsfilename') = ['Return the path to the file in which coordinates are stored.' char(10) ...
        'Assumes that a gray view has been created (though current view' char(10) ...
        'can be any type).' char(10) ...
        'coordsFileName = viewGet(vw, ''coords file name'');'];
        DictViewHelp('ncoords') = ['Return the number of coordinates in the current view. See case' char(10) ...
        '''coords''.' char(10) ...
        'nCoords = viewGet(vw, ''number of coords'');'];
        DictViewHelp('classfilename') = ['Return the path to either the left or the right gray/white' char(10) ...
        'classification file. ' char(10) ...
        'fname = viewGet(vw, ''class file name'', ''left'');' char(10) ...
        'fname = viewGet(vw, ''class file name'', ''right'');'];
        DictViewHelp('classdata') = ['classFileRight = viewGet(vw,''class data'',''right'');'];
        DictViewHelp('graymatterfilename') = ['grayFile = viewGet(vw,''Gray matter filename'',''right'');'];
        DictViewHelp('datavalindex') = '';
        DictViewHelp('analysisdomain') = '';
        DictViewHelp('graycoords') = ['''graycoords'' is also an alias for coords in volume/gray views:'];
        DictViewHelp('leftpath') = '';
        DictViewHelp('rightpath') = '';
        DictViewHelp('fliplr') = '';
        DictViewHelp('imagerotation') = '';
        DictViewHelp('hemifromcoords') = '';
        DictViewHelp('roihemi') = '';
        DictViewHelp('ishidden') = '';
        DictViewHelp('ui') = '';
        DictViewHelp('fignum') = '';
        DictViewHelp('windowhandle') = '';
        DictViewHelp('displaymode') = '';
        DictViewHelp('anatomymode') = '';
        DictViewHelp('coherencemode') = '';
        DictViewHelp('correlationmode') = '';
        DictViewHelp('phasemode') = '';
        DictViewHelp('amplitudemode') = '';
        DictViewHelp('projectedamplitudemode') = '';
        DictViewHelp('mapmode') = '';
        DictViewHelp('zoom') = '';
        DictViewHelp('crosshairs') = '';
        DictViewHelp('locs') = '';
        DictViewHelp('phasecma') = ['This returns only the color part of the map'];
        DictViewHelp('cmapcurrent') = '';
        DictViewHelp('cmapcurmodeclip') = '';
        DictViewHelp('cmapcurnumgrays') = '';
        DictViewHelp('cmapcurnumcolors') = '';
        DictViewHelp('flipud') = ['Return the boolean indicating whether to invert the image u/d in' char(10) ...
        'the graphical user interface.  It is sometimes convenient to do' char(10) ...
        'this in the Inplane view if the top of the slice corresponds to' char(10) ...
        'the bottom of the brain.' char(10) ...
        'Example: ' char(10) ...
        'flipud = viewGet(vw, ''flip updown'');'];
        DictViewHelp('ampclip') = '';
        DictViewHelp('addmesh') = '';
        DictViewHelp('anatinitialize') = '';
        DictViewHelp('cbarrange') = '';
        DictViewHelp('colorbarhandle') = '';
        DictViewHelp('deletemesh') = '';
        DictViewHelp('initdisplaymodes') = '';
        DictViewHelp('leftclassfile') = '';
        DictViewHelp('leftgrayfile') = '';
        DictViewHelp('mainaxishandle') = '';
        DictViewHelp('refreshfn') = '';
        DictViewHelp('rightclassfile') = '';
        DictViewHelp('rightgrayfile') = '';
        DictViewHelp('roi') = '';
        DictViewHelp('roioptions') = '';
        DictViewHelp('spatialgrid') = '';
        DictViewHelp('spatialgrad') = '';
        DictViewHelp('uiimage') = '';
        DictViewHelp('size') = '';

end %if

if DictViewHelp.isKey(paramIn)
    res = DictViewHelp(paramIn);
else
    error('Dict:ViewHelpError', 'The input %s does not appear to be in the dictionary', paramIn);
end %if

return


text:
function e = end(X,k,n)
%END Last index of indexing expression for sptenmat.
%
%   The expression X(end,:) will call END(X,1,2) to determine
%   the value of the first index.
%
%   See also SPTENMAT, SPTENMAT/SUBSREF, END.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if n > 2
  error('Subscript out of range.');
end
e = size(X,k);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%IKINE_SYM  Symbolic inverse kinematics
%
% Q = R.IKINE_SYM(K, OPTIONS) is a cell array (Cx1) of inverse kinematic
% solutions of the SerialLink object ROBOT.  The cells of Q represent the
% different possible configurations.  Each cell of Q is a vector (Nx1), and
% element J is the symbolic expressions for the J'th joint angle.  The
% solution is in terms of the desired end-point pose of the robot which is
% represented by the symbolic matrix (3x4) with elements
%      nx ox ax tx
%      ny oy ay ty
%      nz oz az tz
% where the first three columns specify orientation and the last column
% specifies translation.
%
% K <= N can have only specific values:
%  - 2 solve for translation tx and ty
%  - 3 solve for translation tx, ty and tz
%  - 6 solve for translation and orientation
%
% Options::
%
% 'file',F    Write the solution to an m-file named F
%
% Example::
%
%         mdl_planar2
%         sol = p2.ikine_sym(2);
%         length(sol)
%         ans =
%               2       % there are 2 solutions
%         s1 = sol{1}  % is one solution
%         q1 = s1(1);      % the expression for q1
%         q2 = s1(2);      % the expression for q2
%
% Notes::
% - Requires the Symbolic Toolbox for MATLAB.
% - This code is experimental and has a lot of diagnostic prints.
% - Based on the classical approach using Pieper's method.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function out = ikine_sym(srobot, N, varargin)

    %
    % Given a robot model the following steps are performed:
    % 1. Convert model to symbolic form
    % 2. Find relevant trig equations and solve them for joint angles
    % 3. Write an M-file to implement the solution
    %      xikine(T)
    %      xikine(T, S) where S is a 3 vector with elements 1 or 2 to select
    %       the first or second solution for the corresponding joint.
    %
    % TODO:
    %  - handle the wrist joints, only first 3 joints so far
    %  - handle base and tool transforms

    opt.file = [];
    opt = tb_optparse(opt, varargin);

    % make a symbolic representation of the passed robot
    srobot = sym(srobot);
    q = srobot.gencoords();

    % test N DOF has an allowable value
    switch N
        case 2
        case 3
        case 6
        otherwise
            error('RTB:ikine_sym:badarg', 'Can only solve for 2,3,6 DOF');
    end

    % define symbolic elements of the homogeneous transform
    syms nx ox ax tx
    syms ny oy ay ty
    syms nz oz az tz
    syms d3

    % inits
    Q = {};
    trigsubOld = [];
    trigsubNew = [];

    % loop over each joint variable
    for j=1:N
        fprintf('----- solving for joint %d\n', j);

        % create some equations to sift through
        [left,right] = pieper(srobot, j, 'left');

        % decide which equations to look at
        if j <= 3
            % for first three joints only focus on translational part
            left = left(1:3, 4); left = left(:);
            right = right(1:3, 4); right = right(:);
        else
            % for last three joints only focus on rotational part
            left = left(1:3, 1:3); left = left(:);
            right = right(1:3, 1:3); right = right(:);
        end

        % substitute sin/cos for preceding joint as S/C, essentially removes
        % the joint variables from the equations and treats them as constants.
        if ~isempty(trigsubOld)
            left = subs(left, trigsubOld, trigsubNew);
            right = subs(right, trigsubOld, trigsubNew);
        end

        % then simplify the LHS
        %   do it after the substitution to prevent sum of angle terms being introduced
        left = simplify(left);

        % search for a solveable equation:
        %    function of current joint variable on the LHS
        %    constant element on the RHS
        k = NaN;
        for i=1:length(left)
            if hasonly(left(i), j) && isconstant(right(i))
                k = i;
                break;
            end
        end

        eq = [];

        if ~isnan(k)
            % create the equation to solve: LHS-RHS == 0
            eq = left(k) - right(k);
        else
            % ok, we weren't lucky, try another strategy

            % find all equations:
            %    function of current joint variable on the LHS

            k = [];
            for i=1:length(left)
                % has qj on the left and constant on the right
                if hasonly(left(i), j)
                    k = [k i];
                end
            end

            % hopefully we found two of them
            if length(k) < 2
                continue;
            end

            % we did, lets see if the sum square RHS is constant
            rhs = simplify(right(k(1))^2 + right(k(2))^2); % was simple
            if isconstant( rhs )
                % it is, let's sum and square the LHS
                fprintf('lets square and add %d %d\n', k);

                eq = simplify( expand( left(k(1))^2 + left(k(2))^2 ) ) - rhs; % was simple
            end
        end

        % expand the list of joint variable subsitutions
        fprintf('subs sin/cos q%d for S/C\n', j);
        trigsubOld = [trigsubOld mvar('sin(q%d)', j) mvar('cos(q%d)', j)];
        trigsubNew = [trigsubNew mvar('S%d', j) mvar('C%d', j)];

        if isempty(eq)
            fprintf('cant solve this equation');
            k
            left(k)==right(k)
            error('cant solve');
        end
        % now solve the equation
        if srobot.links(j).isrevolute()
            % for revolute joint it will be a trig equation, do we know how to solve it?
            Q{j} = solve_joint(eq, j );
            if isempty(Q)
                warning('cant solve this kind of equation');
            end
        else
            fprintf('prismatic case\n')
            q = sym( sprintf('q%d', j) );
            Q{j} = solve( eq == 0, q);
        end
    end

    % final simplification
    %  get rid of C^2+S^2 and C^4, S^4 terms
    fprintf('**final simplification pass\n')

    % create a list of simplifications
    %  substitute S^2 = 1-C^2, S^4=(1-C^2)^2
    tsubOld = [];
    tsubNew = [];
    for j=1:N
        tsubOld = [tsubOld mvar('S%d', j)^2 mvar('S%d', j)^4];
        tsubNew = [tsubNew 1-mvar('C%d', j)^2 (1-mvar('C%d', j)^2)^2];
    end

    for j=1:N
        for k=1:5
            % seem to need to iterate this, not quite sure why
            Q{j} = simplify( expand( subs(Q{j}, tsubOld, tsubNew) ) );
        end
    end

    % Q is a cell array of equations for joint variables
    if nargout > 0
        out = Q;
    end

    if ~isempty(opt.file)
        fprintf('**generate MATLAB code\n')
        gencode(Q);
    end
end

%PIEPER Return a set of equations using Pieper's method
%
% [L,R] = pieper(robot, n, which)
%
% If robot has link matrix A1 A2 A3 A4 then returns 12 equations from equating the coefficients of
%
%  A1' T = A2 A3 A4     n=1, which='left'
%  A2' A1' T = A3 A4    n=2, which='left'
%  A3' A2' A1' T = A4   n=3, which='left'
%
%  T A4' = A1 A2 A3     n=1, which='right'
%  T A4' A3' = A1 A2    n=2, which='right'
%  T A4' A3' A2' = A1   n=3, which='right'
%
% Judicious choice of the equations can lead to joint solutions

function [L,R] = pieper(robot, n, which)

    if nargin < 3
        which = 'left';
    end

    syms nx ox ax tx real
    syms ny oy ay ty real
    syms nz oz az tz real

    T = [nx ox ax tx
        ny oy ay ty
        nx oz az tz
        0  0  0  1 ];

    T = inv(robot.base) * T * inv(robot.tool);

    q = robot.gencoords();

    % Create the symbolic A matrices
    for j=1:robot.n
        A{j} = robot.links(j).A(q(j));
    end

    switch which
        case 'left'
            left = T;
            for j=1:n
                left = inv(A{j}) * left ;
            end

            right = eye(4,4);
            for j=n+1:robot.n
                right = right * A{j};
            end

        case 'right'
            left = T;
            for j=1:n
                left = left * inv(A{robot.n-j+1});
            end

            right = eye(4,4);
            for j=1:(robot.n-n)
                right = right * A{j};
            end
    end

    %     left = simple(left);
    %     right = simple(right);

    if nargout == 0
        left == right
    elseif nargout == 1
        L = left;
    elseif nargout == 2
        L = left;
        R = right;
    end
end

%SOLVE_JOINT Solve a trigonometric equation
%
% S = SOLVE_JOINT(EQ, J) solves the equation EQ=0 for the joint variable qJ.
% The result is a cell array of solutions.
%
% The equations must be of the form:
%  A cos(qJ) + B sin(qJ) = 0
%  A cos(qJ) + B sin(qJ) = C
%
% where A, B, C are arbitrarily complex expressions.  qJ can be the only
% joint variable in the expression.

function s = solve_joint(eq, j)

    sinj = mvar('sin(q%d)', j);
    cosj = mvar('cos(q%d)', j);

    A = getcoef(eq, cosj);
    B = getcoef(eq, sinj);

    if isempty(A) || isempty(B)
        warning('don''t know how to solve this kind of equation');
    end

    C = -simplify(eq - A*cosj - B*sinj);  % was simple

    if C == 0
        % A cos(q) + B sin(q) = 0
        s(2) = atan2(A, -B);
        s(1) = atan2(-A, B);
    else
        % A cos(q) + B sin(q) = C
        r = sqrt(A^2 + B^2 - C^2);
        phi = atan2(A, B);

        s(2) = atan2(C, r) - phi;
        s(1) = atan2(C, -r) - phi;
    end
    if nargout == 0
        try
            eval(s)
        catch
            s
        end
    end
end

%MVAR Create a symbolic variable
%
% V = MVAR(FMT, ARGS) is a symbolic variable created using SPRINTF
%
% eg. mvar('q%d', j)
%
% The symbolic is explicitly declared to be real.

function v = mvar(fmt, varargin)

    if isempty(strfind(fmt, '('))
        % not a function
        v = sym( sprintf(fmt, varargin{:}), 'real' );
    else
        v = sym( sprintf(fmt, varargin{:}) );

    end
end

%HASONLY Determine if an expression contains only certain joint variables
%
% S = HASONLY(E L) is true if the joint variables (q1, q2 etc.) in the expression E
% are listed in the vector L.
%
% Eg. hasonly('sin(q1)*cos(q2)*cos(q4)', [1 2 3]) -> true
% Eg. hasonly('sin(q1)*cos(q2)*cos(q4)', [1]) -> false

function s = hasonly(eq, j)

    q = findq(eq);
    if isempty(q)
        s = false;
    else
        s = all(ismember(j, findq(eq)));
    end
end

%ISCONSTANT Determine if an expression is free of joint variables
%
% S = ISCONSTANT(E) is true if the expression E contains no joint variables such
% q1, q2 etc.

function s = isconstant(eq)
    s = isempty(findq(eq));
end

%FINDQ Find the joint variables in expression
%
% Q = FINDQ(E) returns a list of integers indicating the joint variables found
% in the expression E.  For instance an instance of 'q1' would cause a 1 to be
% returned and so on.
%
% Eg. findq('sin(q1)*cos(q2)+S3') -> [1 2]

function q = findq(s)

    q = [];

    for var=symvar(s)
        if isempty(var)
            break
        end
        varname = char(var);
        if varname(1) == 'q'
            q = [q str2num(varname(2:end))];
        end
    end
end

function coef = getcoef(eq, trig)
    z = children( collect(eq, trig) );
    z = children( z(1) );
    coef = z(1);
end

% Output a joint expression to a file

function s = gencode(Q, filename)

    function s = G(s, fmt, varargin)
        s = strvcat(s, sprintf(fmt, varargin{:}));
    end

    s = 'function q = xikine(T, sol)';
    s = G(s, '  if nargin < 2; sol = ones(1, %d); end', length(Q));
    s = G(s, '  px = T(1,4); py = T(2,4); pz = T(3,4);');

    for j=1:3
        Qj = Q{j};   % cast it to subclass
        if length(Qj) == 1
            s = G(s, '  q(%d) = %s', j, matgen2(Qj));
        elseif length(Qj) == 2
            s = G(s, '  if sol(%d) == 1', j);
            s = G(s, '    q(%d) = %s', j, matgen2(Qj(1)));
            s = G(s, '  else');
            s = G(s, '    q(%d) = %s', j, matgen2(Qj(2)));
            s = G(s, '  end');

        end

        s = G(s, '  S%d = sin(q(%d));', j, j);
        s = G(s, '  C%d = cos(q(%d));', j, j);
        s = G(s, ' ');

    end
    s = G(s, 'end');

    fp = fopen(filename, 'w');
    for i=1:numrows(s)
        fprintf(fp, '%s\n', deblank(s(i,:)));
    end
    fclose(fp);

end

% Generate MATLAB code from an expression
%
% Requires a bit of a hack, a subclass of sym (sym2) to do this

function s = matgen2(e)

    s = matgen(sym2(e));

    k = strfind(s, '=');
    s = deblank( s(k+2:end) );
end

author: Allopart

repo: rbpf-gmapping


text:
function coinWrite(prob,filename,type)
%COINWRITE  Write a Mathematical File From Matlab Values
%
% coinWrite(prob,filename,type) writes the problem prob to the file specified
% by filename and problem type specified by type. If you do not specify a
% file extension it will default to the type specified.
%
%   Available File Types:
%       - MPS  [.mps]
%       - QPS  [.qps]
%       - LP   [.lp]  (Appears to be CPLEX LP version - Simplified)
%
% You may specify a full path to the file, or if you specify a filename
% only, it will be written to the current MATLAB directory.
%
% The routines underneath use COIN-OR utilities for File IO. See
% attached EPL License.

%   Copyright (C) 2011 Jonathan Currie (IPL)

%Quick Checks
if(~exist('type','var') || isempty(type))
    %See if contained within file name
    dots = strfind(filename,'.');
    if(isempty(dots))
        error('You must supply the file type to this function');
    else %hopefully got a good ext, will check below
        type = filename(dots(end)+1:end);
    end
end
if(~ischar(filename))
    error('Filename must be a char array!');
end

%Problem Checks
if(~isempty(prob.fun) || ~isempty(prob.nlcon))
    error('You cannot write Nonlinear problems using this interface');
end
if(isa(prob.f,'function_handle') || isa(prob.H,'function_handle'))
    error('f and H must be vectors / matrices!');
end
if(~isempty(prob.Q))
    error('You cannot write quadratic constraints using this interface');
end
if(prob.sense == -1)
    error('You cannot explicity write maximization problems using this interface'); %fix up later
end

%Determine file type
switch(lower(type))
    case 'mps'
        filetype = 'mps';
    case 'qps'
        filetype = 'qps';
    case 'lp'
        filetype = 'lp';
    otherwise
        error('Unknown file type %s',type);
end

%Check problem type vs file type
if(isfield(prob,'type') && ~isempty(prob.type))
    switch(lower(prob.type))
        case {'lp','bilp','milp'}
            if(strcmpi(filetype,'qps'))
                error('LP/BILP/MILPs should be written to MPS or LP files');
            end
        case {'qp','miqp'}
            if(~strcmpi(filetype,'qps'))
                error('QP/MIQPs should be written to QPS files');
            end
        otherwise
            error('Problem type %s is not supported by this interface',upper(prob.type));
    end
end

% Make sure we have file extension
if(isempty(strfind(filename,['.' filetype])))
    filename = [filename '.' filetype];
end

% If filename includes absolute path, we can skip which
if(~isempty(strfind(filename,':')))
    p = filename;
else %Locate the full path to the file
    p = [cd filesep filename];
end

%Convert to coin problem
if(~isempty(prob.b) || ~isempty(prob.beq))
    [prob.A,prob.rl,prob.ru] = gen2row(prob.A,prob.b,prob.Aeq,prob.beq);
end
if(isempty(prob.A))
    prob.A = spalloc(0,length(prob.f),0); %error thrown internally if not ok
end

%Ensure sparsity
prob.A = sparse(prob.A);
prob.H = sparse(prob.H);

%Ensure Bounds
if(isempty(prob.lb))
    prob.lb = -Inf(size(prob.f));
end
if(isempty(prob.ub))
    prob.ub = Inf(size(prob.f));
end

%Setup Integer Vars
if(isstruct(prob.int))
    prob.int = prob.int.str;
end
if(isempty(prob.int))
    prob.int = zeros(size(prob.f));
elseif(ischar(prob.int)) %char string
    prob.int = lower(prob.int);
    i = zeros(size(prob.f));
    int = prob.int == 'i';
    bin = prob.int == 'b';
    i(int) = 1;
    i(bin) = 1;
    %Some problems have fixed binary vars (i.e. 0 <= bvar_n <= 0), can't
    %just blanket modify the bounds
    lb_idx = prob.lb(bin) ~= 0 & prob.lb(bin) ~= 1;
    ub_idx = prob.ub(bin) ~= 0 & prob.ub(bin) ~= 1;
    if(any(lb_idx))
        prob.lb(lb_idx) = 0;
    end
    if(any(ub_idx))
        prob.ub(ub_idx) = 1;
    end
    prob.int = i;
elseif(length(prob.int) ~= length(prob.f) || any(prob.int > 1)) %find indicies (might be wrong though!)
    i = zeros(size(prob.f));
    i(prob.int) = 1;
    prob.int = i;
end
prob.int = int8(prob.int);

%Setup SOS
if(~isempty(prob.sos))
    [r,c] = size(prob.sos.type);
    if(r > c)
        prob.sos_type = str2num(prob.sos.type); %#ok<ST2NM>
    else
        prob.sos_type = str2num(prob.sos.type'); %#ok<ST2NM>
    end
    if(length(prob.sos.type) == 1)
        prob.sos_index = {prob.sos.index};
        prob.sos_weight = {prob.sos.weight};
    else
        prob.sos_index = prob.sos.index;
        prob.sos_weight = prob.sos.weight;
    end
else
    prob.sos_type = [];
    prob.sos_index = [];
    prob.sos_weight = [];
end
%Setup Objc
if(~isempty(prob.objbias) && prob.objbias == 0)
    prob.objbias = 0;
end

%Call coin-or mex
coinW(prob,p,filetype);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function sparse_grid_mixed_size_tests ( )

%*****************************************************************************80
%
%% SPARSE_GRID_MIXED_SIZE_TESTS calls SPARSE_GRID_MIXED_SIZE_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  tol = sqrt ( eps );

  fprintf ( 1, '\n' );
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_SIZE_TESTS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Call SPARSE_GRID_MIXED_SIZE_TEST with various arguments.\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %e\n', tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 1 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 3 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 4 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 7 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 8 ]';
  alpha = [ 0.0, 1.5 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 2, 9 ]';
  alpha = [ 0.0, 0.5 ]';
  beta = [ 0.0, 1.5 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 6, 4 ]';
  alpha = [ 2.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );

  dim_num = 3;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 1, 2, 5 ]';
  alpha = [ 0.0, 0.0, 0.0 ]';
  beta = [ 0.0, 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Rule 13
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 13, 13 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Rule 16
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 2;
  rule = [ 16, 16 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Dimension 2, Rule 17
%
  dim_num = 2;
  level_max_min = 0;
  level_max_max = 4;
  rule = [ 17, 17 ]';
  alpha = [ 0.0, 0.0 ]';
  beta = [ 0.0, 0.0 ]';
  sparse_grid_mixed_size_test ( dim_num, level_max_min, level_max_max, ...
    rule, alpha, beta, tol );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_SIZE_TESTS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,nCasePerBixel,visBool)
% matRad ompMC monte carlo photon dose calculation wrapper
%
% call
%   dij = matRad_calcPhotonDoseMc(ct,stf,pln,cst,visBool)
%
% input
%   ct:                         matRad ct struct
%   stf:                        matRad steering information struct
%   pln:                        matRad plan meta information struct
%   cst:                        matRad cst struct
%   visBool:                    binary switch to enable visualization
% output
%   dij:                        matRad dij struct
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2018 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg =  MatRad_Config.instance();

tic

% disable visualiazation by default
if nargin < 6
    visBool = false;
end

if nargin < 5
    nCasePerBixel = matRad_cfg.propMC.ompMC_defaultHistories;
    matRad_cfg.dispInfo('Using default number of Histories per Bixel: %d\n',nCasePerBixel);
end

fileFolder = fileparts(mfilename('fullpath'));

if ~matRad_checkMexFileExists('omc_matrad') %exist('matRad_ompInterface','file') ~= 3
    matRad_cfg.dispWarning('Compiled mex interface not found. Trying to compile the ompMC interface on the fly!');
    try
        matRad_compileOmpMCInterface();
    catch MException
        matRad_cfg.dispError('Could not find/generate mex interface for MC dose calculation.\nCause of error:\n%s\n Please compile it yourself (preferably with OpenMP support).',MException.message);
    end
end

matRad_calcDoseInit;

% gaussian filter to model penumbra from (measured) machine output / see diploma thesis siggel 4.1.2
if isfield(machine.data,'penumbraFWHMatIso')
    penumbraFWHM = machine.data.penumbraFWHMatIso;
else
    penumbraFWHM = 5;
    matRad_cfg.dispWarning('photon machine file does not contain measured penumbra width in machine.data.penumbraFWHMatIso. Assuming 5 mm.');
end

sourceFWHM = penumbraFWHM * machine.meta.SCD/(machine.meta.SAD - machine.meta.SCD);
sigmaGauss = sourceFWHM / sqrt(8*log(2)); % [mm]

% set up arrays for book keeping
dij.bixelNum = NaN*ones(dij.totalNumOfBixels,1);
dij.rayNum   = NaN*ones(dij.totalNumOfBixels,1);
dij.beamNum  = NaN*ones(dij.totalNumOfBixels,1);

dij.numHistoriesPerBeamlet = nCasePerBixel;

omcFolder = [matRad_cfg.matRadRoot filesep 'ompMC'];
%omcFolder = [matRad_cfg.matRadRoot filesep 'submodules' filesep 'ompMC'];

%% Setup OmpMC options / parameters

%display options
ompMCoptions.verbose = matRad_cfg.logLevel - 1;

% start MC control
ompMCoptions.nHistories = nCasePerBixel;
ompMCoptions.nSplit = 20;
ompMCoptions.nBatches = 10;
ompMCoptions.randomSeeds = [97 33];

%start source definition
ompMCoptions.spectrumFile       = [omcFolder filesep 'spectra' filesep 'mohan6.spectrum'];
ompMCoptions.monoEnergy         = 0.1;
ompMCoptions.charge             = 0;
ompMCoptions.sourceGeometry     = 'gaussian';
ompMCoptions.sourceGaussianWidth = 0.1*sigmaGauss;

% start MC transport
ompMCoptions.dataFolder   = [omcFolder filesep 'data' filesep];
ompMCoptions.pegsFile     = [omcFolder filesep 'pegs4' filesep '700icru.pegs4dat'];
ompMCoptions.pgs4formFile = [omcFolder filesep 'pegs4' filesep 'pgs4form.dat'];

ompMCoptions.global_ecut = 0.7;
ompMCoptions.global_pcut = 0.010;

% Relative Threshold for dose
ompMCoptions.relDoseThreshold = 1 - matRad_cfg.propDoseCalc.defaultLateralCutOff;

% Output folders
ompMCoptions.outputFolder = [omcFolder filesep 'output' filesep];

% Create Material Density Cube
material = cell(4,5);
material{1,1} = 'AIR700ICRU';
material{1,2} = -1024;
material{1,3} = -974;
material{1,4} = 0.001;
material{1,5} = 0.044;
material{2,1} = 'LUNG700ICRU';
material{2,2} = -974;
material{2,3} = -724;
material{2,4} = 0.044;
material{2,5} = 0.302;
material{3,1} = 'ICRUTISSUE700ICRU';
material{3,2} = -724;
material{3,3} = 101;
material{3,4} = 0.302;
material{3,5} = 1.101;
material{4,1} = 'ICRPBONE700ICRU';
material{4,2} = 101;
material{4,3} = 1976;
material{4,4} = 1.101;
material{4,5} = 2.088;

% conversion from HU to densities & materials
for s = 1:dij.numOfScenarios

    HUcube{s} =  matRad_interp3(dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z,ct.cubeHU{s}, ...
                            dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest');

    % projecting out of bounds HU values where necessary
    if max(HUcube{s}(:)) > material{end,3}
        matRad_cfg.dispWarning('Projecting out of range HU values');
        HUcube{s}(HUcube{s}(:) > material{end,3}) = material{end,3};
    end
    if min(HUcube{s}(:)) < material{1,2}
        matRad_cfg.dispWarning('Projecting out of range HU values');
        HUcube{s}(HUcube{s}(:) < material{1,2}) = material{1,2};
    end

    % find material index
    cubeMatIx{s} = NaN*ones(dij.doseGrid.dimensions,'int32');
    for i = size(material,1):-1:1
        cubeMatIx{s}(HUcube{s} <= material{i,3}) = i;
    end

    % create an artificial HU lookup table
    hlut = [];
    for i = 1:size(material,1)
        hlut = [hlut;material{i,2} material{i,4};material{i,3}-1e-10 material{i,5}]; % add eps for interpolation
    end

    cubeRho{s} = interp1(hlut(:,1),hlut(:,2),HUcube{s});

end

ompMCgeo.material = material;

scale = 10; % to convert to cm

ompMCgeo.xBounds = (dij.doseGrid.resolution.y * (0.5 + [0:dij.doseGrid.dimensions(1)])) ./ scale;
ompMCgeo.yBounds = (dij.doseGrid.resolution.x * (0.5 + [0:dij.doseGrid.dimensions(2)])) ./ scale;
ompMCgeo.zBounds = (dij.doseGrid.resolution.z * (0.5 + [0:dij.doseGrid.dimensions(3)])) ./ scale;

%% debug visualization
if visBool

    figure
    hold on

    axis equal

    % ct box
    ctCorner1 = [ompMCgeo.xBounds(1) ompMCgeo.yBounds(1) ompMCgeo.zBounds(1)];
    ctCorner2 = [ompMCgeo.xBounds(end) ompMCgeo.yBounds(end) ompMCgeo.zBounds(end)];
    plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner1(3)],'k' )
    plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
    plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
    plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' )
    plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner2(3) ctCorner2(3)],'k' )
    plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
    plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
    plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' )
    plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' )
    plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' )
    plot3([ctCorner1(1) ctCorner1(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' )
    plot3([ctCorner2(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' )

    xlabel('x [cm]')
    ylabel('y [cm]')
    zlabel('z [cm]')

    rotate3d on

end

%% Create beamlet source
useCornersSCD = true; %false -> use ISO corners

numOfBixels = [stf(:).numOfRays];
beamSource = zeros(dij.numOfBeams, 3);

bixelCorner = zeros(dij.totalNumOfBixels,3);
bixelSide1 = zeros(dij.totalNumOfBixels,3);
bixelSide2 = zeros(dij.totalNumOfBixels,3);

counter = 0;

for i = 1:dij.numOfBeams % loop over all beams

    % define beam source in physical coordinate system in cm
    beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10;

    for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray!

        counter = counter + 1;

        dij.beamNum(counter)  = i;
        dij.rayNum(counter)   = j;
        dij.bixelNum(counter) = j;

        if useCornersSCD
            beamletCorners = stf(i).ray(j).rayCorners_SCD;
        else
            beamletCorners = stf(i).ray(j).beamletCornersAtIso;
        end

        % get bixel corner and delimiting vectors.
        % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm
        currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ scale;
        bixelCorner(counter,:) = currCorner;
        bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ scale - currCorner;
        bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ scale - currCorner;

        if visBool
            for k = 1:4
                currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10;
                % rays connecting source and ray corner
                plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b')
                % connection between corners
                lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10;
                plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r')
            end
        end

    end

end

ompMCsource.nBeams = dij.numOfBeams;
ompMCsource.iBeam = dij.beamNum(:);

% Switch x and y directions to match ompMC cs.
ompMCsource.xSource = beamSource(:,2);
ompMCsource.ySource = beamSource(:,1);
ompMCsource.zSource = beamSource(:,3);

ompMCsource.nBixels = sum(numOfBixels(:));
ompMCsource.xCorner = bixelCorner(:,2);
ompMCsource.yCorner = bixelCorner(:,1);
ompMCsource.zCorner = bixelCorner(:,3);

ompMCsource.xSide1 = bixelSide1(:,2);
ompMCsource.ySide1 = bixelSide1(:,1);
ompMCsource.zSide1 = bixelSide1(:,3);

ompMCsource.xSide2 = bixelSide2(:,2);
ompMCsource.ySide2 = bixelSide2(:,1);
ompMCsource.zSide2 = bixelSide2(:,3);

if visBool
    plot3(ompMCsource.ySource,ompMCsource.xSource,ompMCsource.zSource,'rx')
end

%% Call the OmpMC interface

%ompMC for matRad returns dose/history * nHistories.
% This factor calibrates to 1 Gy in a %(5x5)cm^2 open field (1 bixel) at
% 5cm depth for SSD = 900 which corresponds to the calibration for the
% analytical base data.
absCalibrationFactor = 3.49056 * 1e12; %Approximate!

%Now we have to calibrate to the the beamlet width.
absCalibrationFactor = absCalibrationFactor * (pln.propStf.bixelWidth/50)^2;

matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation... \n');

outputVariance = matRad_cfg.propMC.ompMC_defaultOutputVariance;

if isfield(pln,'propMC') && isfield(pln.propMC,'outputVariance')
    outputVariance = pln.propMC.outputVariance;
end

%run over all scenarios
for s = 1:dij.numOfScenarios
    ompMCgeo.isoCenter = [stf(:).isoCenter];

    %Run the Monte Carlo simulation and catch  possible mex-interface
    %issues
    try
        %If we ask for variance, a field in the dij will be filled
        if outputVariance
            [dij.physicalDose{s},dij.physicalDose_MCvar{s}] = omc_matrad(cubeRho{s},cubeMatIx{s},ompMCgeo,ompMCsource,ompMCoptions);
        else
            [dij.physicalDose{s}] = omc_matrad(cubeRho{s},cubeMatIx{s},ompMCgeo,ompMCsource,ompMCoptions);
        end
    catch ME
        errorString = [ME.message '\nThis error was thrown by the MEX-interface of ompMC.\nMex interfaces can raise compatability issues which may be resolved by compiling them by hand directly on your particular system.'];
        matRad_cfg.dispError(errorString);
    end

    %Calibrate the dose with above factor
    dij.physicalDose{s} = dij.physicalDose{s} * absCalibrationFactor;
    if isfield(dij,'physicalDose_MCvar')
        dij.physicalDose_MCvar{s} = dij.physicalDose_MCvar{s} * absCalibrationFactor^2;
    end
end

matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n');
matRad_cfg.dispInfo(evalc('toc'));

try
    % wait 0.1s for closing all waitbars
    allWaitBarFigures = findall(0,'type','figure','tag','TMWWaitbar');
    delete(allWaitBarFigures);
    pause(0.1);
catch
end

end

repo: matRad


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtMshClean.m                                 |
%|    #    |   VERSION    : 0.53                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2019                                    |
%| ( === ) |   SYNOPSIS   : Clean tetrahedral mesh                        |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Load mesh
load tetmesh

% Build mesh
mesh1 = msh(X,tet);

% Graphical representation
figure
plot(mesh1)
axis equal
view(45,45)

% Add vertices to tetrahedral mesh
tmp = [zeros(10,3) ; 10*rand(10,3)];
tet = tet + size(tmp,1);
X   = [tmp ; X ; 10*rand(20,3)];

% Build and clean mesh
mesh2 = msh(X,tet);

% Test mesh egality
~isequal(mesh1,mesh2)

% Graphical representation
figure
plot(mesh2)
axis equal
view(45,45)

% Colours
mesh = mshSquare(20,[1 1]);
ctr  = mesh.ctr;
mesh.col(ctr(:,1)<0)  = 1;
mesh.col(ctr(:,1)>=0) = 2;

% Graphical representation
figure
plot(mesh)
axis equal
view(45,45)

% Sub-meshing with small translation
delta     = 1e-2;
mesh1     = mesh.sub(mesh.col==1);
mesh2     = mesh.sub(mesh.col==2);
mesh2.vtx = mesh2.vtx + delta;
meshT     = union(mesh1,mesh2);

% Graphical representation
figure
plot(meshT)
axis equal
view(45,45)

% Clean with specified range
stp   = mesh.stp;
meshT = clean(meshT,0.5*stp(1));

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function display(tree)
% XMLTREE/DISPLAY Command window display of an XMLTree
% FORMAT display(tree)
%
% tree - XMLTree object
%__________________________________________________________________________
%
% This method is called when the semicolon is not used to terminate a
% statement which returns an XMLTree.
%__________________________________________________________________________
% Copyright (C) 2002-2011  http://www.artefact.tk/

% Guillaume Flandin
% $Id: display.m 4460 2011-09-05 14:52:16Z guillaume $

disp(' ');
disp([inputname(1),' = ']);
disp(' ');
for i=1:numel(tree)
    disp([blanks(length(inputname(1))+3) char(tree(i))]);
end
disp(' ');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [R,T] = slabVibro(f, rho, c, e)
% Copyright (c) 2018-2019, Marc Bakry, Ecole Polytechnique
% GNU General Public License v3.0.
% Computation of the coeff. of reflection / transmission for an acaustic wave
% crossing three environments with the parameters rho_ {1, 2, 3} and
% c_ {1,2, 3} for a set of frequencies f
% e.g. ~/nonRegressionTest/vibroAcoustic/transmission_1d_calcul
% To be commented...
assert(length(rho) == 3 && size(c, 1) == length(f) && size(c, 2) == 3)
omega = 2*pi*f;
R = zeros(size(f)); T = zeros(size(f));
muk = (ones(length(f), 1)*rho) .* c;
for i=1:length(f)
    tmp = muk(i, :)*omega(i);
    k   = omega(i)./c(i, :);
    b   = [-1; tmp(1); 0; 0];
    A   = [1 -1 -1 0; ...
        tmp(1) -tmp(2) tmp(2) 0; ...
        0 exp(-1i*k(2)*e) exp(1i*k(2)*e) -exp(1i*k(3)*e); ...
        0 tmp(2)*exp(-1i*k(2)*e) ...
        -tmp(2)*exp(1i*k(2)*e) tmp(3)*exp(1i*k(3)*e)];
    x = A\b;
    R(i) = x(1);
    T(i) = x(4);
end
end


text:
%%% Extraction of the final intrinsic and extrinsic paramaters:

check_active_images;

if ~exist('solution_error')
   solution_error = zeros(6*n_ima + 15,1);
end;

fc = solution(1:2);%***
cc = solution(3:4);%***
alpha_c = solution(5);%***
kc = solution(6:9);%***

fc_error = solution_error(1:2);
cc_error = solution_error(3:4);
alpha_c_error = solution_error(5);
kc_error = solution_error(6:9);

% Calibration matrix:

KK = [fc(1) fc(1)*alpha_c cc(1);0 fc(2) cc(2); 0 0 1];
inv_KK = inv(KK);

% Extract the extrinsic paramters, and recomputer the collineations

for kk = 1:n_ima,

   if active_images(kk),

      omckk = solution(15+6*(kk-1) + 1:15+6*(kk-1) + 3);%***
      Tckk = solution(15+6*(kk-1) + 4:15+6*(kk-1) + 6);%***

      omckk_error = solution_error(15+6*(kk-1) + 1:15+6*(kk-1) + 3);
      Tckk_error = solution_error(15+6*(kk-1) + 4:15+6*(kk-1) + 6);

   	Rckk = rodrigues(omckk);

   	Hkk = KK * [Rckk(:,1) Rckk(:,2) Tckk];

   	Hkk = Hkk / Hkk(3,3);

   else

      omckk = NaN*ones(3,1);
      Tckk = NaN*ones(3,1);
      Rckk = NaN*ones(3,3);
      Hkk = NaN*ones(3,3);
      omckk_error = NaN*ones(3,1);
      Tckk_error = NaN*ones(3,1);

   end;

   eval(['omc_' num2str(kk) ' = omckk;']);
   eval(['Rc_' num2str(kk) ' = Rckk;']);
   eval(['Tc_' num2str(kk) ' = Tckk;']);
   eval(['H_' num2str(kk) '= Hkk;']);
   eval(['omc_error_' num2str(kk) ' = omckk_error;']);
   eval(['Tc_error_' num2str(kk) ' = Tckk_error;']);

end;

repo: ekfmonoslam


text:
function [B,B2,dist] = parzen(data, mu, Sigma, N)
% EVAL_PDF_COND_PARZEN Evaluate the pdf of a conditional Parzen window
% function B = eval_pdf_cond_parzen(data, mu, Sigma, N)
%
% B(q,t) = Pr(data(:,t) | Q=q) = sum_{m=1}^{N(q)} w(m,q)*K(data(:,t) - mu(:,m,q); sigma)
% where K() is a Gaussian kernel with spherical variance sigma,
% and w(m,q) = 1/N(q) if m<=N(q) and = 0 otherwise
% where N(q) is the number of mxiture components for q
%
% B2(m,q,t) =  K(data(:,t) - mu(:,m,q); sigma) for m=1:max(N)

% This is like eval_pdf_cond_parzen, except mu is mu(:,m,q) instead of mu(:,q,m)
% and we use 1/N(q) instead of mixmat(q,m)

if nargout >= 2
  keep_B2 = 1;
else
  keep_B2 = 0;
end

if nargout >= 3
  keep_dist = 1;
else
  keep_dist = 0;
end

[d M Q] = size(mu);
[d T] = size(data);

M = max(N(:));

B = zeros(Q,T);
const1 = (2*pi*Sigma)^(-d/2);
const2 = -(1/(2*Sigma));
if T*Q*M>20000000 % not enough memory to call sqdist
  disp('eval parzen for loop')
  if keep_dist,
    dist = zeros(M,Q,T);
  end
  if keep_B2
    B2 = zeros(M,Q,T);
  end
  for q=1:Q
    D = sqdist(mu(:,1:N(q),q), data); % D(m,t)
    if keep_dist
      dist(:,q,:) = D;
    end
    tmp = const1 * exp(const2*D);
    if keep_B2,
      B2(:,q,:) = tmp;
    end
    if N(q) > 0
      %B(q,:) = (1/N(q)) * const1 * sum(exp(const2*D), 2);
      B(q,:) = (1/N(q)) * sum(tmp,1);
    end
  end
else
  %disp('eval parzen vectorized')
  dist = sqdist(reshape(mu(:,1:M,:), [d M*Q]), data); % D(mq,t)
  dist = reshape(dist, [M Q T]);
  B2 = const1 * exp(const2*dist); % B2(m,q,t)
  if ~keep_dist
    clear dist
  end

  % weights(m,q) is the weight of mixture component m for q
  %    = 1/N(q) if m<=N(q) and = 0 otherwise
  % e.g., N = [2   3   1], M = 3,
  % weights = [1/2 1/3 1   = 1/2 1/3 1/1      2 3 1     1 1 1
  %            1/2 1/3 0     1/2 1/3 1/1 .*   2 3 1 <=  2 2 2
  %            0   1/3 0]    1/2 1/3 1/1      2 3 1     3 3 3

  Ns = repmat(N(:)', [M 1]);
  ramp = 1:M;
  ramp = repmat(ramp(:), [1 Q]);
  n = N + (N==0); % avoid 1/0 by replacing with 0* 1/1m where 0 comes from mask
  N1 = repmat(1 ./ n(:)', [M 1]);
  mask = (ramp <= Ns);
  weights = N1 .* mask;
  B2 = B2 .* repmat(mask, [1 1 T]);

  % B(q,t) = sum_m B2(m,q,t) * P(m|q) = sum_m B2(m,q,t) * weights(m,q)
  B = squeeze(sum(B2 .* repmat(weights, [1 1 T]), 1));
  B = reshape(B, [Q T]); % undo effect of squeeze in case Q = 1
end


text:
function net = fcnInitializeModel16s(net)
%FCNINITIALIZEMODEL16S Initialize the FCN-16S model from FCN-32

%% Remove the last layer
net.removeLayer('deconv32') ;

%% Add the first Deconv layer
filters = single(bilinear_u(4, 1, 21)) ;
net.addLayer('deconv2', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 2, ...
                      'crop', [1 1 1 1], ...
                      'hasBias', false), ...
             'x38', 'x39', 'deconv1f') ;
f = net.getParamIndex('deconv1f') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

%% Add skip layers on top of pool4
net.addLayer('skip4', ...
     dagnn.Conv('size', [1 1 512 21], 'pad', 0), ...
     'x24', 'x40', {'skip4f','skip4b'});

f = net.getParamIndex('skip4f') ;
net.params(f).value = zeros(1, 1, 512, 21, 'single') ;
net.params(f).learningRate = 0.1 ;
net.params(f).weightDecay = 1 ;

f = net.getParamIndex('skip4b') ;
net.params(f).value = zeros(1, 1, 21, 'single') ;
net.params(f).learningRate = 2 ;
net.params(f).weightDecay = 1 ;

%% Add the summation layer
net.addLayer('sum1', dagnn.Sum(), {'x39', 'x40'}, 'x41') ;

%% Add deconvolutional layer implementing bilinear interpolation
filters = single(bilinear_u(32, 21, 21)) ;
net.addLayer('deconv16', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 16, ...
                      'crop', 8, ...
                      'numGroups', 21, ...
                      'hasBias', false, ...
                      'opts', net.meta.cudnnOpts), ...
             'x41', 'prediction', 'deconv16f') ;

f = net.getParamIndex('deconv16f') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

% Make sure that the output of the bilinear interpolator is not discared for
% visualization purposes
net.vars(net.getVarIndex('prediction')).precious = 1 ;

% empirical test
if 0
  figure(100) ; clf ;
  n = numel(net.vars) ;
  for i=1:n
    vl_tightsubplot(n,i) ;
    showRF(net, 'input', net.vars(i).name) ;
    title(sprintf('%s', net.vars(i).name)) ;
    drawnow ;
  end
end

author: vlfeat

repo: matconvnet-fcn


text:
function visualizeBoundary(X, y, model, varargin)
%VISUALIZEBOUNDARY plots a non-linear decision boundary learned by the SVM
%   VISUALIZEBOUNDARYLINEAR(X, y, model) plots a non-linear decision
%   boundary learned by the SVM and overlays the data on it

% Plot the training data on top of the boundary
plotData(X, y)

% Make classification predictions over a grid of values
x1plot = linspace(min(X(:,1)), max(X(:,1)), 100)';
x2plot = linspace(min(X(:,2)), max(X(:,2)), 100)';
[X1, X2] = meshgrid(x1plot, x2plot);
vals = zeros(size(X1));
for i = 1:size(X1, 2)
   this_X = [X1(:, i), X2(:, i)];
   vals(:, i) = svmPredict(model, this_X);
end

% Plot the SVM boundary
hold on
contour(X1, X2, vals, [0 0], 'Color', 'b');
hold off;

end

author: Borye


text:
function ic = check_incircle_edge(vertex, face, edge)

% check_incicle_edge - compute "empty circle" property for a set of edges
%
%   ic = check_incicle_edge(vertex,face, edge);
%
%   ic(i)==1 if edge(:,i) is delaunay valid (boundary or empty circles or non convex).
%   It thus should be flipped if ic(i)==0.
%
%   Copyright (c) 2008 Gabriel Peyre

if nargin<3
    edge = compute_edges(face);
end

n = size(edge, 2);
e2f = compute_edge_face_ring(face);

m = size(e2f,1);
f1 = full( e2f(edge(1,:)+(edge(2,:)-1)*m) );
f2 = full( e2f(edge(2,:)+(edge(1,:)-1)*m) );

if not( isempty( find(f1==0 & f2==0) ) )
    warning('Problem with triangulation');
end

ic = ones(n,1);
isconvex = ones(n,1);

I = find(f1>0 & f2>0);

f1 = f1(I); f2 = f2(I);

f1 = face(:,f1);
f2 = face(:,f2);

edge = edge(:,I);
n = length(I);

points1 = find_other(f2, edge(1,:), edge(2,:));
ic1 = check_incircle(vertex, f1, points1 );

points2 = find_other(f1, edge(1,:), edge(2,:));
ic2 = check_incircle(vertex, f2, points2 );

ic(I) = ic1 & ic2;

% check for convexity
v = vertex(:,points2) - vertex(:,points1);
u1 = vertex(:,edge(1,:)) - vertex(:,points1);
u2 = vertex(:,edge(2,:)) - vertex(:,points1);
A = zeros(2,2,n); B = zeros(2,2,n);
A(:,1,:) = reshape(v, [2 1 n]);
A(:,2,:) = reshape(u1, [2 1 n]);
B(:,1,:) = reshape(v, [2 1 n]);
B(:,2,:) = reshape(u2, [2 1 n]);
isconvex(I) = ( det3(A) .* det3(B) ) < 0;
% non convex should not be flipped
ic(isconvex==0) = 1;

%%
function f = find_other(f,a,b)

u = f - repmat(a, [3 1]); f(u==0) = 0;
u = f - repmat(b, [3 1]); f(u==0) = 0;
if sum( sum(f==0)~=2 )>0
    error('Problem with triangulation');
end
f = sum(f);

author: gpeyre

repo: matlab-toolboxes


text:
function [X, Y, Z] = faceLocations(m)
% [X, Y, Z] = faceLocations(m)
% m is a mesh type
% This function returns the X, Y, and Z
% each one is a face variable itself, and can be used for the calculation
% of face variable as a function of locations
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%  cellLocations

% Written by Ali A. Eftekhari
% See the license file
X=createFaceVariable(m, 0);
Y=createFaceVariable(m, 0);
Z=createFaceVariable(m, 0);
N=m.dims;
d = m.dimension;
switch d
case {1, 1.5, 1.8}
    	X.xvalue= m.facecenters.x;
    case {2, 2.5, 2.8}
        X.xvalue= repmat(m.facecenters.x, 1, N(2));
        X.yvalue= repmat(m.cellcenters.y', N(1)+1, 1);
        Y.xvalue= repmat(m.cellcenters.x, 1, N(2)+1);
        Y.yvalue= repmat(m.facecenters.y', N(1), 1);
    case {3, 3.2}
        X.xvalue= repmat(m.facecenters.x, 1, N(2), N(3));
        X.yvalue= repmat(m.cellcenters.y', N(1)+1, 1, N(3));
        z=zeros(1,1,N(3));
        z(1,1,:)= m.cellcenters.z;
        X.zvalue= repmat(z, N(1)+1, N(2), 1);
        % Y
        Y.xvalue= repmat(m.cellcenters.x, 1, N(2)+1, N(3));
        Y.yvalue= repmat(m.facecenters.y', N(1), 1, N(3));
        Y.zvalue= repmat(z, N(1), N(2)+1, 1);
        % Z
        z=zeros(1,1,N(3)+1);
        z(1,1,:)= m.facecenters.z;
        Z.xvalue= repmat(m.cellcenters.x, 1, N(2), N(3)+1);
        Z.yvalue= repmat(m.cellcenters.y', N(1), 1, N(3)+1);
        Z.zvalue= repmat(z, N(1), N(2), 1);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [m,v]=v_chimv(n,l,s)
%V_CHIMV approximate mean and variance of non-central chi distribution [m,v]=(n,l,s)
%
%  Inputs:	n = degrees of freedom
%           l = non-centrality parameter = sqrt(sum(mean^2)) [default 0]
%               (can be a vector or matrix to calculate many different values at once)
%           s = standard deviation of Gaussian [default 1]
%
% Outputs:  m = mean of chi distribution
%           v = variance of chi distribution
%
% If x=c+randn(n,1) is a column vector of Gaussian random numbers with mean vector c, then
% z=sqrt(x'*x) has a chi distributon with n degrees of freedom and non-centrality parameter
% l=sqrt(c'*c). The mean and variance of a chi distribution are given precisely by
%
%     m = sqrt(2)*exp(gammaln(0.5*n+0.5)-gammaln(0.5*n))*hypergeom(-0.5,0.5*n,-0.5*(l/s)^2)*s
%       = sqrt(pi/2) L(0.5,0.5*n-1,-0.5*(l/s)^2)*s
%     v = n*s^2+l^2-m^2
%
% where L(n,a,x) is the generalized Laguerre polynomial L_n^{(a)}(x) but this is very slow
% to calculate so this routine approximates these expressions.
% The expressions are from [1]  together with the following identities taken
% from [2] or [3]: 13.2.39, 13.6.19, 5.2.5, 5.4.1, 5.5.1, 5.4.6.
%
% For n=1, the accuracy is high; for n>1, accuracy improves with increasing n.
% Accuracy is worst when the non-centrality parameter, l, is close to s*sqrt(n).
% Worst case errors as a function of n are:
%                       n:    1       2      3       5      10
%   worst case error in m:  1e-15   0.007  0.004  0.0015  0.0005
%
% References:
% [1]	J. H. Park. Moments of the generalized Rayleigh distribution. Quarterly of Applied Mathematics, 19: 45–49, 1961.
% [2]	F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.
% [3]	F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Digital Library of Mathematical Functions. 2010-2014. URL http://dlmf.nist.gov/.

%      Copyright (C) Mike Brookes 2014-2020
%      Version: $Id: v_chimv.m 11260 2020-07-18 20:07:58Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent ab pp qq nab
if isempty(ab)
    nab=200; % cache a few low values of n
    pp=[ 0.595336298258636  -1.213013700592756  -0.018016200037799   1.999986150447582 0];
    qq=[ -0.161514114798972   0.368983655790737  -0.136992134476950  -0.499681107630725 2];
    ni=1./(1:nab);
    ab=[polyval(qq,ni);polyval(pp,ni)];
end
if nargin<3
    s=1;
    if nargin<2
        l=0;
    end
end
ls=l/s;
l2=(ls).^2;
s2=s^2;
if n<=nab
    if n==1
        m=l.*(1-2*normcdf(-ls))+2*s*normpdf(-ls);
    else
        m=sqrt(l2+n-1+(ab(1,n)+ab(2,n)*l2).^(-1))*s;
    end
else
    m=sqrt(l2+n-1+(polyval(qq,1/n)+polyval(pp,1/n)*l2).^(-1))*s;
end
v=(n+l2)*s2-m.^2;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function out = rval(f)
%RVAL  Obtain the value of a ONEFUN at its right endpoint.
%   RVAL(F) returns the F(1).

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

out = feval(f, 1);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function plot_path(paths, mapClass, SFCs)
% PLOT_PATH Visualize a path and sfc's ellipses
%   PLOT_PATH(path mapClass, SFCs) creates a figure showing a path through
%   the environment.
%   SFCs provide functions for drawing ellipses: drawEllipsoids();

persistent map decomps init;
if isempty(init)
    map = mapClass;
    decomps = SFCs;
    init = 1;
end

figure('Name','Animation');

hold on;
for i = 1:size(map.blocks,1)
    block = map.blocks(i, :);
    drawBlock(block);
end

path = paths{1};

% draw Inflated obstacle
for i = 1:size(map.blocks,1)
    margin = map.margin;
    block = map.blocks(i, :);
    block(1,1:3) = block(1,1:3) - margin;
    block(1,4:6) = block(1,4:6) + margin;
    block(7) = max(150, block(7));  block(8) = max(150, block(8));   block(9) = max(150, block(9));     %% set color
    drawBlock(block);
end

if size(path,1) > 0
    % pcshow(path, [0,1,0],'MarkerSize', 0.1);
    % The lower left corner position, width and height. 260 here is exactly 7cm
%     set(gcf,'position', [500 500 260 220]);
%     set(gcf,'position', [500 500 400 400]);   % big pic
%     set(gcf,'position', [500 500 300 300]);   % small pic
    set(gcf,'color','w');
    set(gca,'color','w');

    %% Display the unit of x y z
%     xlabel('x length(m)');
%     ylabel('y length(m)');
%     zlabel('z length(m)');

end

% Set the axis range
axis_add = 2*map.res(1);
axis([map.boundary.ld(1)-axis_add, map.boundary.ru(1)+axis_add,   ...
      map.boundary.ld(2)-axis_add, map.boundary.ru(2)+axis_add,   ...
      map.boundary.ld(3)-axis_add, map.boundary.ru(3)+axis_add])

hold on
for j = 1 : size(paths, 2)
    if j == 1
        % JPS's result
%         plot3(paths{j}(:,1), paths{j}(:,2), paths{j}(:,3),'-*','color','k','LineWidth',3);
    elseif j == 2
        % JPS() + simple path() result:
%         plot3(paths{j}(:,1), paths{j}(:,2), paths{j}(:,3),'-*','color','#A2142F','LineWidth',3);
        plot3(paths{j}(:,1), paths{j}(:,2), paths{j}(:,3),'-*','color','k','LineWidth',3);
    end
end
hold off

set(gca,'DataAspectRatio',[0.1 0.1 0.1]);

% drawEllipsoid = true;
if (exist('drawEllipsoid') && (drawEllipsoid == true))
    decomps{1}.drawEllipsoids();
%     decomps{2}.drawEllipsoids();
end

grid on

view(30, 10);
hold off;

end

function drawBlock(block)

    x = [ones(4,1) * block(1); ones(4,1) * block(4)];
    y = [ones(2,1) * block(5); ones(2,1) * block(2); ones(2,1) * block(5); ones(2,1) * block(2)];
    z = [block(3);block(6);block(3);block(6);block(3);block(6);block(3);block(6)];

    vert = [x, y, z];
    fac = [1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8;1 2 3 4;5 6 7 8];
    c = block(7:9)/255;
    patch('Vertices',vert,'Faces',fac,...
          'FaceVertexCData',hsv(6),'FaceColor',c,'FaceAlpha',.2);

    x = [ones(4,1) * block(1); ones(4,1) * block(4)];
    y = [block(2);block(5);block(2);block(5);block(2);block(5);block(2);block(5)];
    z = [ones(2,1) * block(3); ones(2,1) * block(6); ones(2,1) * block(3); ones(2,1) * block(6)];

    vert = [x, y, z];
    fac = [1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8;1 2 3 4;5 6 7 8];
    c = block(7:9)/255;
    patch('Vertices',vert,'Faces',fac,...
          'FaceVertexCData',hsv(6),'FaceColor',c,'FaceAlpha',.2);
end

author: LenaShengzhen

repo: AerialRobotics


text:
function db_insert_blob(tablename, fieldname, s)

% DB_INSERT_BLOB converts a Matlab variable of arbitrary type into
% a binary stream and inserts this stream into a binary blob in the
% database table.
%
% Use as
%   db_insert_blob(tablename, fieldname, s)

% Copyright (C) 2007, Robert Oostenveld
%
% $Log: db_insert_blob.m,v $
% Revision 1.1  2007/11/07 10:50:45  roboos
% created helper functions for easy access to a MySQL database table using a structure for representing the data
%

format = 'serialize';
cmd    = sprintf('insert into %s ( %s ) values ( "{M}" )', tablename, fieldname);
blob   = struct2msg(s, format);

% execute the query, this requires the mym version of the mysql mex file
mysql(str, blob);


text:
function f = p06_fun ( option, n, x )

%*****************************************************************************80
%
%% P06_FUN evaluates the integrand for problem 6.
%
%  Discussion:
%
%    The exact value is (m-1)!! * sqrt ( pi ) / sqrt ( 2**m ).
%
%    Integral ( -oo < x < +oo ) x^m exp (-x*x) dx
%
%    The parameter M is set by calling P06_PARAM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION:
%    0, integrand is f(x).
%    1, integrand is exp(-x*x) * f(x);
%    2, integrand is exp(-x*x/2) * f(x);
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real F(N), the function values.
%
  x = x ( : );
  f = zeros ( n, 1 );

  m = 0;
  m = p06_param ( 'G', 'M', m );

  f(1:n) = x(1:n).^m;

  if ( option == 0 )
    f(1:n) = f(1:n) .* exp ( - x(1:n).^2 );
  elseif ( option == 1 )

  elseif ( option == 2 )
    f(1:n) = f(1:n) .* exp ( - 0.5 * x(1:n).^2 );
  end

  return
end


text:
function str = makeValidFieldName(str)
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    isoct=exist('OCTAVE_VERSION','builtin');
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',toascii(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' ))) return;  end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos)) return; end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',toascii(str0(pos(i))))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end

author: AvaisP


text:
function x = rbmup(rbm, x)
    x = sigm(repmat(rbm.c', size(x, 1), 1) + x * rbm.W');
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function img = computeColor(u,v)

middlebury = 1;

%   computeColor color codes flow field U, V

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 21:20:30 (Wed, 31 Oct 2006) $

% Copyright 2007, Deqing Sun.
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

nanIdx = isnan(u) | isnan(v);
u(nanIdx) = 0;
v(nanIdx) = 0;

colorwheel = makeColorwheel();

if middlebury == 0
  [M N] = size(u);
  r = zeros(M,N);
  g = zeros(M,N);
  b = zeros(M,N);
  img = zeros(M,N,3);

  magnitude = sqrt(u.*u + v.*v);
  mm = max(magnitude(:));
  magnitude = magnitude/mm;

  ncols = size(colorwheel, 1);
  a = atan2(u,v)/pi; % in [-1, 1]
  fk = (a+1) /2 * (ncols-1) + 1; % in [1 ncols]

  idx = floor(fk);

  for i=1:ncols
    tmp = idx == i;
    r(tmp) = magnitude(tmp).*colorwheel(i,1);
    g(tmp) = magnitude(tmp).*colorwheel(i,2);
    b(tmp) = magnitude(tmp).*colorwheel(i,3);
  end
  img(:,:,1) = r;
  img(:,:,2) = g;
  img(:,:,3) = b;
else

  ncols = size(colorwheel, 1);

  rad = sqrt(u.^2+v.^2);

  a = atan2(-v, -u)/pi;

  fk = (a+1) /2 * (ncols-1) + 1;  % -1~1 maped to 1~ncols

  k0 = floor(fk);                 % 1, 2, ..., ncols

  k1 = k0+1;
  k1(k1==ncols+1) = 1;

  f = fk - k0;

  for i = 1:size(colorwheel,2)
    tmp = colorwheel(:,i);
    col0 = tmp(k0)/255;
    col1 = tmp(k1)/255;
    col = (1-f).*col0 + f.*col1;

    idx = rad <= 1;
    col(idx) = 1-rad(idx).*(1-col(idx));    % increase saturation with radius

    col(~idx) = col(~idx)*0.75;             % out of range

    img(:,:, i) = uint8(floor(255*col.*(1-nanIdx)));
  end;
end

%%
function colorwheel = makeColorwheel()

%   color encoding scheme

%   adapted from the color circle idea described at
%   http://members.shaw.ca/quadibloc/other/colint.htm

RY = 15;
YG = 6;
GC = 4;
CB = 11;
BM = 13;
MR = 6;

ncols = RY + YG + GC + CB + BM + MR;

colorwheel = zeros(ncols, 3); % r g b

col = 0;
%RY
colorwheel(1:RY, 1) = 255;
colorwheel(1:RY, 2) = floor(255*(0:RY-1)/RY)';
col = col+RY;

%YG
colorwheel(col+(1:YG), 1) = 255 - floor(255*(0:YG-1)/YG)';
colorwheel(col+(1:YG), 2) = 255;
col = col+YG;

%GC
colorwheel(col+(1:GC), 2) = 255;
colorwheel(col+(1:GC), 3) = floor(255*(0:GC-1)/GC)';
col = col+GC;

%CB
colorwheel(col+(1:CB), 2) = 255 - floor(255*(0:CB-1)/CB)';
colorwheel(col+(1:CB), 3) = 255;
col = col+CB;

%BM
colorwheel(col+(1:BM), 3) = 255;
colorwheel(col+(1:BM), 1) = floor(255*(0:BM-1)/BM)';
col = col+BM;

%MR
colorwheel(col+(1:MR), 3) = 255 - floor(255*(0:MR-1)/MR)';
colorwheel(col+(1:MR), 1) = 255;

author: thomas-koehler

repo: SupER