text:
function [grid, cfg] = ft_prepare_sourcemodel(cfg, vol, sens)

% FT_PREPARE_SOURCEMODEL helps to make a source model that can be
% used for source reconstruction, beamformer scanning, linear
% estimation and MEG interpolation.
%
% Use as
%   [grid, cfg] = prepare_dipole_grid(cfg)
% where the configuration structure contains the details on how the source
% model should be constructed.
%
% A source model can be constructed based on
%   - regular 3D grid with explicit specification
%   - regular 3D grid with specification of the resolution
%   - regular 3D grid, based on segmented MRI, restricted to gray matter
%   - surface grid based on brain surface from volume conductor
%   - surface grid based on head surface from external file
%   - using user-supplied grid positions, which can be regular or irregular
%   - cortical sheet that was created in MNE or Freesurfer
% The approach that will be used depends on the configuration options that
% you specify.
%
% Configuration options for generating a regular 3-D grid
%   cfg.grid.xgrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.grid.ygrid      = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.grid.zgrid      = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.grid.resolution = number (e.g. 1 cm) for automatic grid generation
%
% Configuration options for a predefined grid
%   cfg.grid.pos        = N*3 matrix with position of each source
%   cfg.grid.dim        = [Nx Ny Nz] vector with dimensions in case of 3-D grid (optional)
%   cfg.grid.inside     = vector with indices of the sources inside the brain (optional)
%   cfg.grid.outside    = vector with indices of the sources outside the brain (optional)
% The following fields are not used in this function, but will be copied along to the output
%   cfg.grid.leadfield
%   cfg.grid.filter or alternatively cfg.grid.avg.filter
%   cfg.grid.subspace
%   cfg.grid.lbex
%
% Configuration options for cortex segmentation, i.e. for placing dipoles in grey matter
%   cfg.mri           = can be filename, MRI structure or segmented MRI structure
%   cfg.sourceunits   = 'mm' or 'cm' (default is 'cm')
%   cfg.mriunits      = 'mm' or 'cm' (default is 'mm')
%   cfg.threshold     = 0.1, relative to the maximum value in the segmentation
%   cfg.smooth        = 5, smoothing in voxels
%
% Configuration options for reading a cortical sheet from file
%   cfg.headshape     = string, should be a *.fif file
%
% Other configuration options
%   cfg.vol          = volume conduction model
%   cfg.grad         = gradiometer definition
%   cfg.elec         = electrode definition
%   cfg.grid.tight   = 'yes' or 'no' (default is automatic)
%   cfg.inwardshift  = depth of the bounding layer for the source space, relative to the head model surface (default = 0)
%   cfg.symmetry     = 'x', 'y' or 'z' symmetry for two dipoles, can be empty (default = [])
%   cfg.headshape    = a filename containing headshape, a structure containing a
%                      single triangulated boundary, or a Nx3 matrix with surface
%                      points
%
% See also FT_SOURCEANALYSIS, FT_MEGREALIGN, FT_DIPOLEFITTING, FT_PREPARE_LEADFIELD

% Searching through the code, it seems that the following cfg fields are being used
% cfg.grid
% cfg.mri
% cfg.headshape
% cfg.tightgrid
% cfg.symmetry
% cfg.smooth
% cfg.threshold
% cfg.spheremesh
% cfg.inwardshift
% cfg.mriunits
% cfg.sourceunits

% Copyright (C) 2004-2011, Robert Oostenveld
%
% %Log$

% set the defaults
if ~isfield(cfg, 'symmetry'),         cfg.symmetry    = [];       end
if ~isfield(cfg, 'grid'),             cfg.grid        = [];       end

if ~isfield(cfg, 'vol') && nargin>1
  % put it in the configuration structure
  % this is for backward compatibility, 13 Januari 2011
  cfg.vol = vol;
end

if ~isfield(cfg, 'sens') && nargin>2
  % put it in the configuration structure
  % this is for backward compatibility, 13 Januari 2011
  cfg.grad = sens;
end

% for backward compatibility
if isfield(cfg, 'tightgrid')
  cfg.grid.tight = cfg.tightgrid;
  cfg = rmfield(cfg, 'tightgrid');
end

if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'xgrid') && ~ischar(cfg.grid.xgrid)
  error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.xgrid simultaneously');
end
if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'ygrid') && ~ischar(cfg.grid.ygrid)
  error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.ygrid simultaneously');
end
if isfield(cfg.grid, 'resolution') && isfield(cfg.grid, 'zgrid') && ~ischar(cfg.grid.zgrid)
  error('You cannot specify cfg.grid.resolution and an explicit cfg.grid.zgrid simultaneously');
end

% a grid can be constructed based on a number of ways
basedonauto   = isfield(cfg.grid, 'resolution');                        % regular 3D grid with specification of the resolution
basedongrid   = isfield(cfg.grid, 'xgrid') && ~ischar(cfg.grid.xgrid);  % regular 3D grid with explicit specification
basedonpos    = isfield(cfg.grid, 'pos');                               % using user-supplied grid positions, which can be regular or irregular
basedonshape  = isfield(cfg, 'headshape') && ~isempty(cfg.headshape);   % surface grid based on inward shifted head surface from external file
basedonmri    = isfield(cfg, 'mri');                                    % regular 3D grid, based on segmented MRI, restricted to gray matter
basedonvol    = false;                                                  % surface grid based on inward shifted brain surface from volume conductor
basedoncortex = isfield(cfg, 'headshape') && (iscell(cfg.headshape) || ft_filetype(cfg.headshape, 'neuromag_fif') || ft_filetype(cfg.headshape, 'freesurfer_triangle_binary')); % cortical sheet from MNE or Freesurfer, also in case of multiple files/hemispheres

if basedonshape && basedoncortex
  % treating it as cortical sheet has preference
  basedonshape = false;
end

if basedongrid && basedonpos
  % fall back to default behaviour, in which the pos overrides the grid
  basedongrid = false;
end

if ~any([basedonauto basedongrid basedonpos basedonshape basedonmri basedoncortex]) && nargin>1 && ~isempty(cfg.vol)
  % fall back to default behaviour, which is to create a surface grid (e.g. used in MEGRELIGN)
  basedonvol = 1;
end

% these are mutually exclusive, but printing all requested methods here
% facilitates debugging of weird configs. Also specify the defaults here to
% keep the overview
if basedonauto
  fprintf('creating dipole grid based on automatic 3D grid with specified resolution\n');
  if ~isfield(cfg.grid, 'xgrid'),       cfg.grid.xgrid = 'auto';    end
  if ~isfield(cfg.grid, 'ygrid'),       cfg.grid.ygrid = 'auto';    end
  if ~isfield(cfg.grid, 'zgrid'),       cfg.grid.zgrid = 'auto';    end
  if ~isfield(cfg, 'inwardshift'),      cfg.inwardshift = 0;        end % in this case for inside detection, % FIXME move to cfg.grid
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'yes';     end
end

if basedongrid
  fprintf('creating dipole grid based on user specified 3D grid\n');
  if ~isfield(cfg, 'inwardshift'),      cfg.inwardshift = 0;        end % in this case for inside detection, % FIXME move to cfg.grid
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'no';      end
end

if basedonpos
  fprintf('creating dipole grid based on user specified dipole positions\n');
  if ~isfield(cfg, 'inwardshift'),      cfg.inwardshift = 0;        end % in this case for inside detection, % FIXME move to cfg.grid
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'no';      end
end

if basedonshape
  fprintf('creating dipole grid based on inward-shifted head shape\n');
  if ~isfield(cfg, 'spheremesh'),       cfg.spheremesh = 642;  end  % FIXME move spheremesh to cfg.grid
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'no'; end
end

if basedoncortex
  if ~isfield(cfg, 'sourceunits');      cfg.sourceunits = 'cm';     end
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'no'; end
end

if basedonmri
  fprintf('creating dipole grid based on grey matter from segmented MRI\n');
  if ~isfield(cfg, 'threshold'),        cfg.threshold = 0.1;        end % relative
  if ~isfield(cfg, 'smooth');           cfg.smooth    = 5;          end % in voxels
  if ~isfield(cfg, 'sourceunits');      cfg.sourceunits = 'cm';     end
  if ~isfield(cfg, 'mriunits');         cfg.mriunits = 'mm';        end
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'yes';     end
end

if basedonvol
  fprintf('creating dipole grid based on inward-shifted brain surface from volume conductor model\n');
  if ~isfield(cfg, 'spheremesh'),       cfg.spheremesh = 642;  end   % FIXME move to cfg.grid
  if ~isfield(cfg.grid, 'tight'),       cfg.grid.tight = 'no'; end
end

% these are mutually exclusive
if sum([basedonauto basedongrid basedonpos basedonshape basedonmri basedonvol basedoncortex])~=1
  error('incorrect cfg specification for constructing a dipole grid');
end

needspm = isfield(cfg, 'smooth') && ~strcmp(cfg.smooth, 'no');
if needspm
  % check if SPM is in path, do not add any of them to the path yet
  hasspm2 = ft_hastoolbox('SPM2');
  hasspm8 = ft_hastoolbox('SPM8');

  if ~hasspm2 && ~hasspm8
    % is neither is on the path, try adding SPM8
    try, hasspm8 = ft_hastoolbox('SPM8', 1); end
  end

  if ~hasspm8
    % if SPM8 could not be added, try adding SPM2
    try, ft_hastoolbox('SPM2', 1); end
  end
end

% start with an empty grid
grid = [];

% copy the volume conductor and sensor array out of the cfg
if isfield(cfg, 'vol')
  vol = cfg.vol;
else
  vol = [];
end

% these are mutually exclusive
if isfield(cfg, 'grad')
  sens = cfg.grad;
elseif isfield(cfg, 'elec')
  sens = cfg.elec;
else
  sens = [];
end

if basedonauto
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct a regular 3D grid that spans a box encompassing all electrode
  % or gradiometer coils, this will typically also cover the complete brain
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ischar(cfg.grid.xgrid) && strcmp(cfg.grid.xgrid, 'auto')
    grid.xgrid = floor(min(sens.pnt(:,1))):cfg.grid.resolution:ceil(max(sens.pnt(:,1)));
  end
  if ischar(cfg.grid.ygrid) && strcmp(cfg.grid.ygrid, 'auto')
    grid.ygrid = floor(min(sens.pnt(:,2))):cfg.grid.resolution:ceil(max(sens.pnt(:,2)));
  end
  if ischar(cfg.grid.zgrid) && strcmp(cfg.grid.zgrid, 'auto')
    grid.zgrid = floor(min(sens.pnt(:,3))):cfg.grid.resolution:ceil(max(sens.pnt(:,3)));
  end
  grid.dim   = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  [X, Y, Z]  = ndgrid(grid.xgrid, grid.ygrid, grid.zgrid);
  grid.pos   = [X(:) Y(:) Z(:)];
end

if basedongrid
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % a detailled xgrid/ygrid/zgrid has been specified, the other details
  % still need to be determined
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  grid.xgrid = cfg.grid.xgrid;
  grid.ygrid = cfg.grid.ygrid;
  grid.zgrid = cfg.grid.zgrid;
  grid.dim   = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  [X, Y, Z]  = ndgrid(grid.xgrid, grid.ygrid, grid.zgrid);
  grid.pos   = [X(:) Y(:) Z(:)];
end

if basedonpos
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % a grid is already specified in the configuration, reuse as much of the
  % prespecified grid as possible (but only known objects)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if isfield(cfg.grid, 'pos')
    grid.pos = cfg.grid.pos;
  end
  if isfield(cfg.grid, 'mom')
    grid.mom = cfg.grid.mom;
  end
  if isfield(cfg.grid, 'dim')
    grid.dim = cfg.grid.dim;
  end
  if isfield(cfg.grid, 'xgrid')
    % FIXME is it desirable to have this in the grid?
    grid.xgrid = cfg.grid.xgrid;
  end
  if isfield(cfg.grid, 'ygrid')
    % FIXME is it desirable to have this in the grid?
    grid.ygrid = cfg.grid.ygrid;
  end
  if isfield(cfg.grid, 'zgrid')
    % FIXME is it desirable to have this in the grid?
    grid.zgrid = cfg.grid.zgrid;
  end
  if isfield(cfg.grid, 'dim')
    grid.dim = cfg.grid.dim;
  elseif isfield(grid, 'xgrid') && isfield(grid, 'ygrid') && isfield(grid, 'zgrid')
    grid.dim = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  end
  if isfield(cfg.grid, 'inside')
    grid.inside = cfg.grid.inside;
  end
  if isfield(cfg.grid, 'outside')
    grid.outside = cfg.grid.outside;
  end
  if isfield(cfg.grid, 'lbex')
    grid.lbex = cfg.grid.lbex;
  end
  if isfield(cfg.grid, 'subspace')
    grid.subspace = cfg.grid.subspace;
  end
  if isfield(cfg.grid, 'leadfield')
    grid.leadfield = cfg.grid.leadfield;
  end
  if isfield(cfg.grid, 'filter')
    grid.filter = cfg.grid.filter;
  end
  if isfield(cfg.grid, 'tight')
    grid.tight = cfg.grid.tight;
  end
  % this is not supported any more
  if isfield(cfg.grid, 'avg') && isfield(cfg.grid.avg, 'filter')
    error('please put your filters in cfg.grid instead of cfg.grid.avg');
  end
end

if basedonmri
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct a grid based on the segmented MRI that is provided in the
  % configuration, only voxels in gray matter will be used
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % convert the source/functional data into the same units as the anatomical MRI
  scale = 1;
  switch cfg.sourceunits
    case 'mm'
      scale = scale / 1000;
    case 'cm'
      scale = scale / 100;
    case 'dm'
      scale = scale / 10;
    case 'm'
      scale = scale / 1;
    otherwise
      error('unknown physical dimension in cfg.sourceunits');
  end
  switch cfg.mriunits
    case 'mm'
      scale = scale * 1000;
    case 'cm'
      scale = scale * 100;
    case 'dm'
      scale = scale * 10;
    case 'm'
      scale = scale * 1;
    otherwise
      error('unknown physical dimension in cfg.mriunits');
  end

  if ~isfield(cfg.grid, 'resolution')
    switch cfg.sourceunits
      case 'mm'
        cfg.grid.resolution = 10;
      case 'cm'
        cfg.grid.resolution = 1;
      case 'dm'
        cfg.grid.resolution = 0.1;
      case 'm'
        cfg.grid.resolution = 0.01;
    end
  end

  if ischar(cfg.mri)
    % read the segmentation from file
    mri      = ft_read_mri(cfg.mri);
    mri.gray = double(mri.anatomy);
  elseif isstruct(cfg.mri) && ~isfield(cfg.mri, 'gray')
    % looks like a segmentation that has already been loaded in memory
    mri      = cfg.mri;
    mri.gray = double(mri.anatomy);
  elseif isstruct(cfg.mri) && isfield(cfg.mri, 'gray')
    % looks like a complete segmentation from VOLUMESEGMENT
    mri      = cfg.mri;
    mri.gray = double(mri.gray);
  else
    error('cannot determine the format of the segmentation in cfg.mri');
  end

  % apply a smoothing of a certain amount of voxels
  if ~strcmp(cfg.smooth, 'no');
    fprintf('smoothing gray matter segmentation with %d voxels\n', cfg.smooth);
    spm_smooth(mri.gray, mri.gray, cfg.smooth);
  end

  % determine for each voxel whether it belongs to the cortex
  if isfield(cfg, 'threshold')
    fprintf('thresholding gray matter segmentation at a relative value of %f\n', cfg.threshold);
    head = mri.gray./max(mri.gray(:)) > cfg.threshold;
  else
    error('you must specify cfg.threshold for cortex segmentation');
  end

  ind                 = find(head(:));
  fprintf('%d from %d voxels in the segmentation are marked as cortex (%.0f%%)\n', length(ind), prod(size(head)), 100*length(ind)/prod(size(head)));
  [X,Y,Z]             = ndgrid(1:mri.dim(1), 1:mri.dim(2), 1:mri.dim(3));   % create the grid in MRI-coordinates
  posmri              = [X(ind) Y(ind) Z(ind) ones(length(ind),1)];         % take only the inside voxels
  poshead             = mri.transform * posmri';                            % transform to head coordinates
  poshead             = poshead(1:3,:)';
  posmri              = posmri(:,1:3);
  resolution          = cfg.grid.resolution*scale;                                        % source and mri can be expressed in different units (e.g. cm and mm)
  xgrid               = floor(min(poshead(:,1))):resolution:ceil(max(poshead(:,1)));      % create the grid in head-coordinates
  ygrid               = floor(min(poshead(:,2))):resolution:ceil(max(poshead(:,2)));      % with 'consistent' x,y,z definitions
  zgrid               = floor(min(poshead(:,3))):resolution:ceil(max(poshead(:,3)));
  [X,Y,Z]             = ndgrid(xgrid,ygrid,zgrid);
  pos2head            = [X(:) Y(:) Z(:) ones(length(X(:)),1)]';
  pos2mri             = mri.transform \ pos2head;                                         % transform to MRI-coordinates
  pos2mri             = round(pos2mri(1:3,:))';
  pos2head            = pos2head(1:3,:)';
  pos2mri             = pos2mri(:,1:3);
  % it might be that the box with the points does not completely fit into the MRI
  sel = find(pos2mri(:,1)<1 |  pos2mri(:,1)>size(head,1) | ...
    pos2mri(:,2)<1 |  pos2mri(:,2)>size(head,2) | ...
    pos2mri(:,3)<1 |  pos2mri(:,3)>size(head,3));
  if isempty(sel)
    % use the efficient implementation
    inside = head(sub2ind(mri.dim, pos2mri(:,1), pos2mri(:,2), pos2mri(:,3)));
  else
    % only loop over the points that can be dealt with
    inside = zeros(length(xgrid)*length(ygrid)*length(zgrid), 1);
    for i=setdiff(1:size(pos2mri,1), sel(:)')
      inside(i) = head(pos2mri(i,1), pos2mri(i,2), pos2mri(i,3));
    end
  end
  inside = find(inside);

  grid.pos            = pos2head/scale;                                     % convert to source units
  grid.xgrid          = xgrid/scale;                                        % convert to source units
  grid.ygrid          = ygrid/scale;                                        % convert to source units
  grid.zgrid          = zgrid/scale;                                        % convert to source units
  grid.dim            = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
  grid.inside         = inside(:);
  grid.outside        = setdiff(1:size(grid.pos,1),grid.inside)';

  fprintf('the regular 3D grid encompassing the cortex contains %d grid points\n', size(grid.pos,1));
  fprintf('%d grid points inside gray matter\n', length(grid.inside));
  fprintf('%d grid points outside gray matter\n', length(grid.outside));
end

if basedoncortex
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read it from a *.fif file that was created using Freesurfer and MNE
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if iscell(cfg.headshape)
    % FIXME loop over all files, this should be two hemispheres
    keyboard
  end
  switch ft_filetype(cfg.headshape)
    case 'freesurfer_triangle_binary'
      % it contains a cortical sheet which was created by the Freesurfer software
      shape = ft_read_headshape(cfg.headshape, 'format', 'freesurfer_triangle_binary');
    case 'neuromag_fif'
      % fif files can contain a variety of objects
      % here we assume that it contains a cortical sheet which was created by the MNE software
      shape = ft_read_headshape(cfg.headshape, 'format', 'mne_source');
    otherwise
      % use autodetection
      shape = ft_read_headshape(cfg.headshape);
  end
  grid.pos = shape.pnt;
  grid.tri = shape.tri;
  grid = ft_convert_units(grid, cfg.sourceunits);
end

if basedonshape
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the headshape  to make a superficial dipole layer (e.g.
  % for megrealign). Assume that all points are inside the volume.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get the surface describing the head shape
  if isstruct(cfg.headshape) && isfield(cfg.headshape, 'pnt')
    % use the headshape surface specified in the configuration
    headshape = cfg.headshape;
  elseif isnumeric(cfg.headshape) && size(cfg.headshape,2)==3
    % use the headshape points specified in the configuration
    headshape.pnt = cfg.headshape;
  elseif ischar(cfg.headshape)
    % read the headshape from file
    headshape = ft_read_headshape(cfg.headshape);
  else
    error('cfg.headshape is not specified correctly')
  end
  if ~isfield(headshape, 'tri')
    % generate a closed triangulation from the surface points
    headshape.pnt = unique(headshape.pnt, 'rows');
    headshape.tri = projecttri(headshape.pnt);
  end
  grid.pos = headsurface([], [], 'headshape', headshape, 'inwardshift', cfg.inwardshift, 'npnt', cfg.spheremesh);
  grid.inside  = 1:size(grid.pos,1);
  grid.outside = [];
end

if basedonvol
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the volume conduction model to make a superficial dipole layer (e.g.
  % for megrealign). Assume that all points are inside the volume.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  grid.pos = headsurface(vol, sens, 'inwardshift', cfg.inwardshift, 'npnt', cfg.spheremesh);
  grid.inside  = 1:size(grid.pos,1);
  grid.outside = [];
end

% FIXME use inside_vol instead of this replication of code
% determine the dipole locations inside the brain volume
if ~isfield(grid, 'inside') && ~isfield(grid, 'outside')
  if ft_voltype(vol, 'infinite') || ft_voltype(vol, 'halfspace')
    % an empty vol in combination with gradiometers indicates a magnetic dipole
    % in an infinite vacuum, i.e. all dipoles can be considered to be inside
    grid.inside = 1:size(grid.pos,1);
    grid.outside = [];
    outside = zeros(1,size(grid.pos,1));
    for i =1:size(grid.pos,1);
      invacuum = false;
      dip1 = grid.pos(i,:);
      % condition of dipoles falling in the non conductive halfspace
      invacuum = acos(dot(vol.ori,(dip1-vol.pnt)./norm(dip1-vol.pnt))) < pi/2;
      if invacuum
        outside(i) = 1;
      end
    end
    grid.outside = find(outside);
    grid.inside  = find(~outside);
  else
    if isfield(sens, 'ori') && isfield(sens, 'pnt') && isfield(sens, 'tra')
      % in case of MEG, make a triangulation of the outermost surface
      if isfield(cfg, 'headshape')
        % use the specified headshape to construct the bounding triangulation
        [pnt, tri] = headsurface(vol, sens, 'headshape', cfg.headshape, 'inwardshift', cfg.inwardshift, 'surface', 'skin');
      else
        % use the volume conductor model to construct the bounding triangulation
        [pnt, tri] = headsurface(vol, sens, 'inwardshift', cfg.inwardshift, 'surface', 'skin');
      end
    else
      % in case of EEG, make a triangulation of the innermost surface
      [pnt, tri] = headsurface(vol, sens, 'inwardshift', cfg.inwardshift, 'surface', 'brain');
    end
    % determine the dipole positions that are inside the triangulated surface
    tmp = bounding_mesh(grid.pos, pnt, tri);
    grid.inside  = find(tmp==1);
    grid.outside = find(tmp==0);
  end
elseif ~isfield(grid, 'inside')
  grid.inside = setdiff(1:size(grid.pos,1), grid.outside);
elseif ~isfield(grid, 'outside')
  grid.outside = setdiff(1:size(grid.pos,1), grid.inside);
end

if strcmp(cfg.grid.tight, 'yes')
  fprintf('%d dipoles inside, %d dipoles outside brain\n', length(grid.inside), length(grid.outside));
  fprintf('making tight grid\n');
  xmin = min(grid.pos(grid.inside,1));
  ymin = min(grid.pos(grid.inside,2));
  zmin = min(grid.pos(grid.inside,3));
  xmax = max(grid.pos(grid.inside,1));
  ymax = max(grid.pos(grid.inside,2));
  zmax = max(grid.pos(grid.inside,3));
  xmin_indx = find(grid.xgrid==xmin);
  ymin_indx = find(grid.ygrid==ymin);
  zmin_indx = find(grid.zgrid==zmin);
  xmax_indx = find(grid.xgrid==xmax);
  ymax_indx = find(grid.ygrid==ymax);
  zmax_indx = find(grid.zgrid==zmax);
  sel =       (grid.pos(:,1)>=xmin & grid.pos(:,1)<=xmax); % select all grid positions inside the tight box
  sel = sel & (grid.pos(:,2)>=ymin & grid.pos(:,2)<=ymax); % select all grid positions inside the tight box
  sel = sel & (grid.pos(:,3)>=zmin & grid.pos(:,3)<=zmax); % select all grid positions inside the tight box
  grid.pos   = grid.pos(sel,:);
  % update the inside and outside vector
  tmp = zeros(1,prod(grid.dim));
  tmp(grid.inside)  = 1;        % these are originally inside the brain
  tmp(grid.outside) = 0;        % these are originally outside the brain
  tmp               = tmp(sel); % within the tight box, these are inside the brain
  grid.inside  = find(tmp);
  grid.outside = find(~tmp);
  grid.xgrid   = grid.xgrid(xmin_indx:xmax_indx);
  grid.ygrid   = grid.ygrid(ymin_indx:ymax_indx);
  grid.zgrid   = grid.zgrid(zmin_indx:zmax_indx);
  grid.dim     = [length(grid.xgrid) length(grid.ygrid) length(grid.zgrid)];
end
fprintf('%d dipoles inside, %d dipoles outside brain\n', length(grid.inside), length(grid.outside));

% apply the symmetry constraint, i.e. add a symmetric dipole for each location defined sofar
% set up the symmetry constraints
if ~isempty(cfg.symmetry)
  if strcmp(cfg.symmetry, 'x')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 -1 1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 -x1 y
  elseif strcmp(cfg.symmetry, 'y')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 1 -1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 -y
  elseif strcmp(cfg.symmetry, 'z')
    reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    mirror = [1 1 1 1 1 -1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 y1
  else
    error('unrecognized symmetry constraint');
  end
  fprintf('each source describes two dipoles with symmetry along %s axis\n', cfg.symmetry);
  % expand the number of parameters from one (3) to two dipoles (6)
  grid.pos = grid.pos(:,expand) .* repmat(mirror, size(grid.pos,1), 1);
end

% FIXME should the cfg be added to the output grid?

author: PatternRecognition

repo: OpenBMI


text:
function displaytable(data,headings,labels,cellsize,decimalsize,permutation,cellindices,cellcolors)
%
% Prints out a condensed, legible N-dimensional data table with dimension
% headings and variable cell width and precision.
%
% Entries can be highlighted with color using the cprintf function available
% on Matlab Central.
%
% data is an n-dimensional array
% headings {'Row' 'Column'} is a cell array containing the name of each dimension.
% labels {{1 2 3} {4 5 6}} is a cell array containing lists of values for each dimension.
%   each element of labels can be either:
%     cell array of strings;
%     cell array of numbers;
%     array of numbers.
%
% cellsize (optional, default 6) is the width of each displayed cell
% decimalsize (optional, default 2) is the number of decimal places
%   displayed in each cell
%
% permutation [3 4 1 2] describes the ordering of dimensions in the final
%   displayed table. Data, headings, and labels are all permuted accordingly
%
% cellindices [m x ndim]: m n-dimensional tuplets containing the indices of
%   entries to be highlighted
% cellcolors [m x 3]: [R,G,B] color tuplets for entries to be highlighted
%
% Matt Caywood
%
%   Versions:
% 0.1:  2010.06.16 Initial version
% 0.2:  2010.06.29 Fixed degenerate tables with dimension values = 1
% 0.3:  2010.07.15 Added convenient table permutations
% 0.31: 2010.07.16 Fixed dimension error checking
% 0.4:  2010.07.20 Added color highlighting of entries using cprintf
% 0.41: 2010.07.21 Fixed header spacing
% 0.42: 2011.03.08 Worked around -0 issue in printf

% displaytable(rand(5,5),{'1' '2'},{{'C1' 'C2' 'C3' 'I1' 'I2'},{'C1' 'C2' 'C3' 'I1' 'I2'}})

global cellwidth decimalplaces cellformat myprintf;

%% check input valid
if ~isnumeric(data), error('displaytable:NotNumeric','Data must be numeric.'); end
if ~isreal(data), error('displaytable:NotReal','Data must contain only real elements.'); end

% use headings as canonical length, because matlab ignores singleton trailing dimensions in data
nd = length(headings);
dims = size(data);

if any(dims == 0), fprintf('[]\n'); return; end % some dimension is zero
if ndims(data) > nd, error('displaytable:headings','Not enough headings for data.'); end
if length(labels) ~= nd, error('displaytable:headings','Headings and labels are different lengths.'); end
for i = 1:nd
    % accept singleton trailing dimensions
    if ~( (i > length(dims) && 1 == length(labels{i})) || (dims(i) == length(labels{i})) )
       error('displaytable:labels','Label list #%d (%s) is wrong length for data.',i,headings{i});
    end
end

%% create cell format string
if ~exist('cellsize','var'), cellwidth = 6; else cellwidth = cellsize; end
if ~exist('decimalsize', 'var'), decimalplaces = 2; else decimalplaces = decimalsize; end

assert(cellwidth > decimalplaces+1);
cellformat = ['%' int2str(cellwidth-(decimalplaces+1)) '.' int2str(decimalplaces) 'f'];

%% normalize labels to cell array of strings
for i = 1:nd
    if ~iscell(labels{i}), labels{i} = num2cell(labels{i}); end
    for j = 1:length(labels{i})
        if isnumeric(labels{i}{j}), labels{i}{j} = num2str(labels{i}{j}); end
    end
end

%% Use cprintf if available and cellindices specified, otherwise builtin fprintf
myprintf = @(col,fmt,varargin) fprintf(fmt,varargin{:});

if ~exist('cellindices','var')
    cellindices = [];
    cellcolors = [];
else
    assert(size(cellindices,2) == nd,'Cellindices have wrong dimensions.');
    if ~exist('cprintf','file')
        fprintf('WARNING: cprintf not found, colors will not be shown!\nPlease download cprintf from Matlab Central:\nhttp://www.mathworks.com/matlabcentral/fileexchange/24093\n');
    else
        myprintf = @cprintf;
        if ~exist('cellcolors','var')
            cellcolors = repmat([1 0 0],[size(cellindices,1) 1]); % default red
        else
            assert(size(cellcolors,2) == 3,'Colors must be RGB triplets');
            assert(size(cellcolors,1) == size(cellindices,1),'Colors must be equal length to indices');
        end
    end
end

%% permute table
if exist('permutation','var') && ~isempty(permutation)
    assert(length(permutation) == nd,'Wrong permutation length.');
    data = permute(data,permutation);
    headings = {headings{permutation}};
    labels = {labels{permutation}};

    if ~isempty(cellindices)
        for m = 1:size(cellindices,1)
            cellindices(m,:) = cellindices(m,permutation);
        end
    end
end

% recurse to print out all tables
displaytable_helper(data,headings,labels,cellindices,cellcolors);

end

%% recursion helper
function displaytable_helper(data,headings,labels,cellindices,cellcolors)

global myprintf;

nd = ndims(data); % nd always >= 2

if (nd == 2) % 1D row or column vector, or 2D matrix
    display2dtable(data,headings,labels,cellindices,cellcolors);

else % nd > 2
    title = headings{1};
    sublabels = labels{1};
    labelchars = maxchars(sublabels);

    nidx = size(data,1);
    remainingdims = size(data); remainingdims = remainingdims(2:nd);
    for subscr = 1:nidx

        whichidx = cell(1,nd);
        whichidx{1} = subscr;
        for i = 2:nd; whichidx{i} = ':'; end
        subs.type = '()';
        subs.subs = whichidx;
        subtab = subsref(data,subs);
        subtab = reshape(subtab,remainingdims);

        myprintf('text','%s = %s\t\n',title,truncpadstring(sublabels{subscr},labelchars));
        if isempty(cellindices), passindices = [];
        else passindices = (cellindices(:,1) == subscr);
        end
        displaytable_helper(subtab,headings(2:end),labels(2:end),cellindices(passindices,2:end),cellcolors(passindices,:));
    end
end
end

%% base case
function display2dtable(data,headings,labels,cellindices,cellcolors)

assert(~isempty(headings));
assert(nargin >= 3,'Need row, column labels');
assert(nargin == ndims(data) + 3,'Wrong # arguments');

[nrow,ncol] = size(data);

global cellformat cellwidth decimalplaces myprintf;

if length(headings) == 1
    % handle degenerate 1 x n, n x 1 cases
    if (ncol == 1)
        colheading = '';
        collabels = cell(size(data));
        rowheading = headings{1};
        rowlabels = labels{1};
    elseif (nrow == 1)
        colheading = headings{1};
        collabels = labels{1};
        rowheading = '';
        rowlabels = cell(size(data));
    end
else
    rowheading = headings{1};
    rowlabels = labels{1};
    colheading = headings{2};
    collabels = labels{2};
end

rowchars = max(3,maxchars(rowlabels));

% squeeze dimension label into row label space
dimchars = max(1,floor(rowchars/2));
dim = [truncpadstring(rowheading,dimchars) '/' truncpadstring(colheading,dimchars)];

% display header
myprintf('text',dim);
for i = 1:ncol
    myprintf('text',' %s',truncpadstring(collabels{i},cellwidth));
end
myprintf('text','\n');

% display rows
for j = 1:nrow
    myprintf('text',truncpadstring(rowlabels{j},rowchars));
    for i = 1:ncol

        if isempty(cellindices)
            idx = [];
        else
            idx = find(cellindices(:,1) == j & cellindices(:,2) == i);
            if ~isempty(idx), idx = idx(1); end % use first only
        end

        % clean up the rounding a bit so "-0.00" never appears, only "0.00"
        roundeddata = fround(data(j,i),decimalplaces);
        if (roundeddata == 0), roundeddata = 0; end % this gets rid of floating point -0

        if isempty(idx)
            myprintf('text',' %s',truncpadstring(sprintf(cellformat,roundeddata),cellwidth));
        else
            myprintf(cellcolors(idx,:),' %s',truncpadstring(sprintf(cellformat,roundeddata),cellwidth));
        end
    end
    myprintf('text','\n');
end

end

%% -----------------
function n = maxchars(c)
% find maximum length in characters of strings in cell array

n = 0;
for i = 1:length(c)
    n = max(n,length(c{i}));
end

end

%% -----------------
function ps = truncpadstring(s,n)
% fit string s into n characters

if length(s) > n
    ps = s(1:n);
else
    ps = sprintf(['%' num2str(n) 's'],s);
end

end

%% -----------------
function b = fround(a,n)
% fix-round

if ~exist('n','var'), n = 0; end

b = round(a*10^n)/10^n;

end

author: Sable

repo: mcbench-benchmarks


text:
function [resid] = residual(Xw,Xd,aXw,aXd,sXw,sXd,Y,D)
%
%  'w' refers to time series on wet days
%  'd' refers to time series on dry days
%   X refers to observed time series
%   aX refers to average values of periodic time series (calculated by Fourier analysis)
%   sX refers to standard deviation of periodic time series (Fourier analysis)
%
residw=[];
residd=[];
%
for i = 1:Y
   resw=zeros(1,D);
   resd=resw;
   jw=find(Xw(i,:)~=-999 & Xw(i,:)~=0);
   jd=find(Xd(i,:)~=-999 & Xd(i,:)~=0);
   %
   resw(jw)=(Xw(i,jw)-aXw(jw))./sXw(jw);
   resd(jd)=(Xd(i,jd)-aXd(jd))./sXd(jd);
   %
   residw=[residw resw];
   residd=[residd resd];
   %
   clear resw;
   clear resd;
   %
end
resid=residw+residd;


text:
 function xs = eql_sps_os(x, Gb, yi, ci, ri, R, ...
		niter, pixmax, curv, relax0, chat)
%function xs = eql_sps_os(x, Gb, yi, ci, ri, R, ...
%		niter, pixmax, curv, relax0, chat)
%
% E-QPL-SPS-OS algorithm for emission Poisson problem
% (ordered subsets separable paraboloidal surrogates)
% in
%	x	[np,1]		initial estimate
%	Gb	[nd,np]		Gblock object (see eql_osps_test.m)
%	yi,ci,ri [nb,na]	see em_fbp.m (for model too)
%	R			penalty structure (see Robject.m)
%					(or empty for ML case)
%	niter			# iterations
%	pixmax			upper constraint for pixel values
%	curv		'oc' for erdogan's optimal curvatures
%			'pc' for erdogan's fast precomputed curvatures,
%				which usually provides faster convergence,
%				but can be nonmonotone.
%	relax0	[1] or [2]	relax0 or (relax0, relax_rate)
% out
%	xs [np,niter]	updated image vectors each iteration
%
% Copyright Mar 2000, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

Gb = block_ob(Gb, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Gb, 'n');
starts = subset_start(nblock);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end

if ~isvar('R'), R = []; end
if ~isvar('niter')	| isempty(niter),	niter = 1;	end
if ~isvar('pixmax')	| isempty(pixmax),	pixmax = inf;	end
if ~isvar('curv')	| isempty(curv),	curv = 'oc';	end
if ~isvar('chat')	| isempty(chat),	chat = false;	end

if ~isvar('relax0')	| isempty(relax0),	relax0 = 1;	end
if length(relax0) == 1
	relax_rate = 0;
elseif length(relax0) == 2
	relax_rate = relax0(2);
	relax0 = relax0(1);
else
	error relax
end

eml_check(yi, ci, ri, 'os', nblock);

[nb, na] = size(yi);

gi = sum(Gb')';		% g_i = sum_j g_ij
gi = reshape(gi, nb, na);

%
% precomputed curvatures
%
if streq(curv, 'pc')
	ni = ci.^2 ./ max(yi,1);	% precomputed
	ni = eml_curvature(yi, ci, ri, [], [], 'pc');

	% efficient single denominator consistent with aspire
	denom = Gb' * col(gi .* ni);
elseif ~streq(curv, 'oc')
	error 'curv not implemented'
end

%
% loop over iterations
%
xs = zeros(numel(x), niter);
x = max(x,0);
x = min(x,pixmax);
xs(:,1) = x;

for iter = 2:niter
	if chat, printf('E-QL-SPS-OS iteration %d', iter-1), end

	relax = relax0 / (1 + relax_rate * (iter-2));

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

		li = Gb{iblock} * x;			% l=G*x "line integrals"
		li = reshape(li, nb, length(ia));
		yb = ci(:,ia) .* li + ri(:,ia);		% predicted meas. means

		% fix: need to be careful here with 0/0 -> 0
		dothi = ci(:,ia) .* (yi(:,ia) ./ yb - 1);
		grad = Gb{iblock}' * dothi(:);

		if streq(curv, 'oc')
			% optimal curvatures (for monotone increase)
			ni = eml_curvature(yi(:,ia), ci(:,ia), ri(:,ia), ...
					li, yb, 'oc');
			denom = Gb{iblock}' * col(gi(:,ia) .* ni);
			denom = nblock * denom;
		end

		if isempty(R)
			num = nblock * grad;
			den = denom;
		else
			num = nblock * grad - R.cgrad(R, x);
			den = denom + R.denom(R, x);
		end

		x = x + relax * num ./ den;	% relaxed update
		x = max(x,0);			% lower bound
		x = min(x,pixmax);		% upper bound
	end

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

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function s = freeSpace()

% FREESPACE Check for free space in Map.
%   FREESPACE() returns the number of free states in the map Map.
%
%   See also NEWRANGE, BLOCKRANGE, RANGES.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.
%   Copyright 2015-     Joan Sola @ IRI-CSIC-UPC

global Map

s = sum(~Map.used);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
classdef select_pairs< dagnn.ElementWise
%   In the first form of label_type, C has dimension H x W x 2 x N and specifies a
%   2 categorical label for each spatial location.   C(:,:,1,:) == ones(H,W,N). C(:,:,2,:)== zeros(H,W,N);
%   In the second form of label_type, C has dimension H x W x 1 x N and specifies a
%   categorical label for each spatial location.

  properties
    label_type = 'form2'
  end

  methods
      % --------------------------------------------------------------------
      function y = zerosLike(obj,x)
          % --------------------------------------------------------------------
          if isa(x,'gpuArray')
              y = gpuArray.zeros(size(x),classUnderlying(x)) ;
          else
              y = zeros(size(x),'like',x) ;
          end
      end
%
      function y = onesLike(obj, x)
          % --------------------------------------------------------------------
          if isa(x,'gpuArray')
              y = gpuArray.ones(size(x),classUnderlying(x)) ;
          else
              y = ones(size(x),'like',x) ;
          end
      end

    function outputs = forward(obj, inputs, params)
        score = inputs{1};
        label = inputs{2};
        ind_pos = label>0;
        ind_neg = label<0;
        ind_pos = ind_pos(:);
        ind_neg = ind_neg(:);
        n_pos = sum(ind_pos);
        n_neg = sum(ind_neg);
        ind_select = ones(n_pos,n_neg)>0;
        pos_use = sum(ind_select,2);
        neg_use = sum(ind_select,1);
%         instance_weight = 1./(pos_use*neg_use+eps);
        instance_weight = ones(n_pos,n_neg);
        n_pairs = sum(ind_select(:));
        sz_in = size(score);
        try
            batch_size = sz_in(4);
        catch
            batch_size = 1;
        end
%         out = 0;
        if isa(score,'gpuArray')
              in0 = gpuArray.zeros([n_pairs,2,batch_size],classUnderlying(score)) ;
              instance_weight =gpuArray(instance_weight );
          else
              in0 = zeros([n_pairs,2,batch_size],'single');
        end

        for i = 1:batch_size
            tmp_s = score(:,:,:,i);
            pos = tmp_s(ind_pos);
            neg = tmp_s(ind_neg);
%             ef1 = exp(pos); ef2 = exp(neg);
            pos_m = repmat(pos,1,n_neg).*instance_weight;
            neg_m = repmat(neg',n_pos,1).*instance_weight;
            pos_v = pos_m(ind_select);
            neg_v = neg_m(ind_select);
            if size(pos_v,1)==1
                pos_v = pos_v'; neg_v = neg_v';
            end
            in0(:,:,i) = [pos_v,neg_v];
        end
        in0 = reshape(in0,[1,n_pairs,2,batch_size]);
%         switch obj.label_type
%             case 'form1'
%         label0 = obj.zerosLike(in0);
%         label0(:,:,1,:) = 1;
%             case 'form2'
%                 label0 = ones(1,n_pairs,1,batch_size,'single');
%             otherwise
%
%         end

        outputs{1} = in0;
%         outputs{2} = label0;
%         sz_in = size(inputs{1});
%         n = prod(sz_in(1:3));
%         ef1 = exp(inputs{1});
%         ef2 = exp(inputs{2});
%         out = 2*(ef2./(ef1+ef2+eps)).^2 ;
%       outputs{1} = sum(out(:))/n;
%       n = obj.numAveraged ;
%       m = n + size(inputs{1},4) ;
%       obj.average = (n * obj.average + gather(outputs{1})) / m ;
%       obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
        score = inputs{1};
        label = inputs{2};
        ind_pos = label>0;
        ind_neg = label<0;
        ind_pos = ind_pos(:);
        ind_neg = ind_neg(:);
        n_pos = sum(ind_pos);
        n_neg = sum(ind_neg);
        ind_select0 = ones(n_pos,n_neg,'single');
        ind_select = ind_select0>0;
        pos_use = sum(ind_select,2);
        neg_use = sum(ind_select,1);
%         instance_weight = 1./(pos_use*neg_use+eps);
        instance_weight = ones(n_pos,n_neg);
        if isa(score,'gpuArray')
            ind_select0 = gpuArray(ind_select0);
            instance_weight =gpuArray(instance_weight );
        end
%         n_pairs = sum(ind_select(:));
        sz_in = size(score);
        try
            batch_size = sz_in(4);
        catch
            batch_size = 1;
        end
        der1 = squeeze(derOutputs{1});

%         out = 0;
        derIn0 = obj.zerosLike(score);
%         derOutputs{1} = derOutputs{1}/n_pairs;
        for i = 1:batch_size
            tmp = obj.zerosLike(ind_select0);
            tmp(ind_select) = der1(:,1,i);
            tmp = tmp .*instance_weight;
            tmp_s = score(:,:,:,i);
            der_tmp = obj.zerosLike(tmp_s);
            der_tmp(ind_pos) = sum(tmp,2);

            tmp = obj.zerosLike(ind_select0);
            tmp(ind_select) = der1(:,2,i);
            tmp = tmp .*instance_weight;
%             der_tmp = obj.zerosLike(tmp_s);
            der_tmp(ind_neg) = sum(tmp,1);
            derIn0(:,:,1,i) = reshape(der_tmp,size(tmp_s));

%             pos = tmp_s(ind_pos);
%             neg = tmp_s(ind_neg);
%             ef1 = exp(pos); ef2 = exp(neg);
%             ef1_m = repmat(ef1,1,n_neg).*ind_select;
%             ef2_m = repmat(ef2',n_pos,1).*ind_select;
%             derIn1 =  -4.*ef2_m.^2./((ef1_m+ef2_m).^3+eps).*derOutputs{1};
%             der_tmp(ind_pos) = sum(derIn1,2);
%             derIn2 = derIn1 *(-1).* ef1_m;
%             der_tmp(ind_neg) = sum(derIn2);
%             derIn0(:,:,1,i) = reshape(der_tmp,size(tmp_s));
        end
        derInputs{1} = derIn0;
        derInputs{2} = [];
%         sz_in = size(inputs{1});
%         n = prod(sz_in(1:3));
%         derOutputs{1} = derOutputs{1}/n;
%         ef1 = exp(inputs{1});
%         ef2 = exp(inputs{2});
%       derInputs{1} = -4.*ef2.^2./((ef1+ef2).^3+eps).*derOutputs{1};
%       derInputs{2} = derInputs{1}*(-1).* ef1;
      derParams = {} ;
    end

%     function reset(obj)
%       obj.average = 0 ;
%       obj.numAveraged = 0 ;
%     end

    function outputSizes = getOutputSizes(obj, inputSizes, paramSizes)
      outputSizes{1} = [1 1 1 inputSizes{1}(4)] ;
    end

    function rfs = getReceptiveFields(obj)
      % the receptive field depends on the dimension of the variables
      % which is not known until the network is run
      rfs(1,1).size = [NaN NaN] ;
      rfs(1,1).stride = [NaN NaN] ;
      rfs(1,1).offset = [NaN NaN] ;
      rfs(2,1) = rfs(1,1) ;
    end

    function obj = select_pairs(varargin)
      obj.load(varargin) ;
    end
  end
end

author: shenjianbing

repo: TripletTracking


text:
function newmodel = addGenes(model,geneIDs,varargin)
% Add Genes (or a single gene) to the model
%
% USAGE:
%    model = addGenes(model,geneIDs,varargin)
%
% INPUT:
%    model:             The model to add the Metabolite batch to.
%    geneIDs:           The IDs of the genes that shall be added.
%    varargin:          fieldName, Value pairs. The given fields will be set
%                       according to the values. Only defined COBRA fields may be
%                       used. The following fields will be ignored (as they
%                       are dependent on the existing model structure):
%
%                        * rxnGeneMat - The size will be extended, but no associations will be set.
% OUTPUT:
%    newmodel:     The model structure with the additional reactions.
%
% EXAMPLE:
%    %To add genes with specific fields use:
%    `model = addGenes(model, {'G1', 'Gene2', 'InterestingGene'}, 'proteins', {'Protein1','Protein B','Protein Alpha'}, 'geneField2', {'D','E','F'})`
%
% NOTE:
%    All fields (geneField1/geneField2) have to be present in
%    the model, or defined in the field definitions, otherwise they are
%    ignored.
%

if (isfield(model,'genes') && ~isempty(model.genes) && any(ismember(model.genes,geneIDs))) || numel(unique(geneIDs)) < numel(geneIDs)
    error('Duplicate Reaction ID detected.');
end

if ~isfield(model,'genes') % Should not happen but if, lets add it.
    model.genes = cell(0,1);
end

nGenes = numel(model.genes);

%We have make sure, that the new fields are in sync, so we create those
%first.
fieldDefs = getDefinedFieldProperties();
fieldDefs = fieldDefs(cellfun(@(x) strcmp(x,'genes'), fieldDefs(:,2)) | cellfun(@(x) strcmp(x,'genes'), fieldDefs(:,3)));
modelGeneFields = getModelFieldsForType(model,'genes');
model.genes = [model.genes;columnVector(geneIDs)];

for field = 1:2:numel(varargin)
    %If the field does not exist and is not defined, we will not add that
    %information, as we don't know how the field should look.
    cfield = varargin{field};
    if strcmp('rxnGeneMat',cfield) || (~any(ismember(fieldDefs(:,1),cfield)) && ~any(ismember(modelGeneFields,cfield)))
        warning('Field %s is excluded.',cfield);
        continue;
    end

    if ~any(size(varargin{field+1}) == numel(geneIDs)) %something must fit
        error('Size of field %s does not fit to the rxnList size', varargin{field});
    end
    %Now add the field data.
    if ~isfield(model,cfield)
        %If its not yet in the model, create it
        %according to its defaults and replace the final elements.
        model = createEmptyFields(model,cfield);
        model.(cfield)((end-numel(varargin{field+1})+1):end) = columnVector(varargin{field+1});
    else
        %Or just extend it with the supplied data.
        model.(cfield) = [model.(cfield);columnVector(varargin{field+1})];
    end

end
%Extend the remaining fields, to keep the fields in sync
newmodel = extendModelFieldsForType(model,'genes','originalSize',nGenes);

author: opencobra

repo: cobratoolbox


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

syms q1 q2 q3
% matrices DH
A01 = dh_sym(q1, 0, 0.5, 0);
A12 = dh_sym(q2, 0, 0.5, 0);
A23 = dh_sym(q3, 0, 0.5, 0);

A02 = A01*A12;
A03 = A01*A12*A23;
A03 = simplify(A03);

z0 = [0 0 1]';

p03=A03(1:3,4);
p13=A03(1:3,4)-A01(1:3,4);
p23=A03(1:3,4)-A02(1:3,4);

Jv = [cross(z0, p03) cross(z0, p13) cross(z0, p23)];
Jv = simplify(Jv)

Jw = [z0 z0 z0];
J = [Jv; Jw]
% transform to valid dimensions
J = [J(1:2,:); J(6,:)]

J = simplify(J)

singularities = det(J)
singularities = simplify(singularities)

%Redefine J
J = [J(1:2,:)]
% compute moore-penrose pseudo inverse
Jp = J'*inv(J*J')
Jp = simplify(Jp)

% compute projector P
P = eye(3)-Jp*J

function A = dh_sym(theta, d, a, alpha)
syms q1 q2 q3 q4
% avoid almost zero elements in cos(alpha) and sin(alpha)
ca = cos(alpha);
sa = sin(alpha);
if abs(ca) < 1e-6
    ca = 0;
end
if abs(sa) < 1e-6
    sa = 0;
end

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

repo: ARTE


text:
% The COBRAToolbox: testMatrixCoherence.m
%
% Purpose:
%     - testMatriceCoherence tests the functionality of matrixCoherence.
%
% Authors:
%     - Sylvain Arreckx March 2017
%
% Exemple comes from http://stemblab.github.io/mutual-coherence/

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testMatrixCoherence'));
cd(fileDir);

% define the tolerance
tol = 1e-5;

A = [1, -1, 1;
     1,  2, 4];

[mu, Q] = matrixCoherence(A);

assert(abs(mu - 0.8575) < tol);

% change the directory
cd(currentDir)


text:
function [minVal, minPos] = min(f)
%MIN   Global minimum of a SINGFUN on [-1,1].
%   MINVAL = MIN(F) returns the global minimum of the SINGFUN F on [-1,1].
%
%   [MINVAL, MINPOS] = MIN(F) returns also a value such that MINVAL = F(MINPOS).
%
% See also MAX, MINANDMAX.

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

% To avoid code duplication, we simply call MINANDMAX():
[minVal, minPos] = minandmax(f);
minVal = minVal(1,:);
minPos = minPos(1,:);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function u = tan(a)
%TAN          Slope tangent tan(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = tan(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  Xxs = hull(a.r(:,indexc),a.r(:,indexr));
  u.s = a.s ./ sqr(cos(Xxs));

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function qr = plane_normal_xyz_to_qr ( pp, normal, pq, pr, n, xyz )

%*****************************************************************************80
%
%% PLANE_NORMAL_XYZ_TO_QR: XYZ to QR coordinates for a normal form plane.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      NORMAL is a normal vector to the plane.
%
%    Two vectors PQ and PR can be computed with the properties that
%    * NORMAL, PQ and PR are pairwise orthogonal;
%    * PQ and PR have unit length;
%    * every point P in the plane has a "QR" representation
%      as P = PP + q * PQ + r * PR.
%
%    This function is given the XYZ coordinates of a set of points on the
%    plane, and returns the QR coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.
%
%    Input, real NORMAL(3), a normal vector N to the plane.  The
%    vector must not have zero length, but it is not necessary for N
%    to have unit length.
%
%    Input, real PQ(3), a vector of unit length,
%    perpendicular to the vector N and the vector PR.
%
%    Input, real PR(3), a vector of unit length,
%    perpendicular to the vector N and the vector PQ.
%
%    Input, integer N, the number of points on the plane.
%
%    Input, real XYZ(3,N), the XYZ coordinates of the points.
%
%    Output, real QR(2,N), the QR coordinates of the points.
%
  qr = [ pq'; pr' ] * ( xyz(1:3,1:n) - repmat ( pp, 1, n ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% EX_LAPLACE_BELTRAMI: solve the Laplace-Beltrami equation on a NURBS surface, discretized with B-splines (non-isoparametric approach).

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS structure from the NURBS toolbox
problem_data.geo_name = 'geo_roof.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4];

% Physical parameters
problem_data.c_diff  = @(x, y, z) ones(size(x));

% Source and boundary terms
alp = 5; bet = 1; L = 1;
g1 = @(x,y) (1-cos(atan(x./y))).*(1-sin(atan(x./y)));
g2 = @(x,y) cos(atan(x./y)) + sin(atan(x./y)) - 4*cos(atan(x./y)).*sin(atan(x./y));
gz = @(z) sin(pi*alp/L*z);
dg1 = @(x,y) sin(atan(y./x)).*(1-sin(atan(y./x))) - cos(atan(y./x)).*(1-cos(atan(y./x)));

problem_data.f = @(x, y, z) bet*((alp*pi/L)^2 * g1(x,y) - g2(x,y)).*gz(z);
problem_data.h = @(x, y, z, ind) bet * g1(x,y).*gz(z);

% Exact solution (optional)
problem_data.uex     = @(x, y, z) bet * g1(x,y).*gz(z);
problem_data.graduex = @(x, y, z) cat (1, ...
                reshape (-1./sqrt(x.^2+y.^2).*sin(atan(y./x))*bet.*gz(z).*dg1(x,y), [1, size(x)]), ...
                reshape ( 1./sqrt(x.^2+y.^2).*cos(atan(y./x))*bet.*gz(z).*dg1(x,y), [1, size(x)]), ...
                reshape ( alp*pi/L * cos(alp*pi/L*z) .* bet .* g1(x,y), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];       % Degree of the splines
method_data.regularity = [2 2];       % Regularity of the splines
method_data.nsub       = [9 9];       % Number of subdivisions
method_data.nquad      = [4 4];       % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_laplace (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW

output_file = 'Laplace-Beltrami.vts';

vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) PLOT IN MATLAB. COMPARISON WITH THE EXACT SOLUTION

[eu, F] = sp_eval (u, space, geometry, vtk_pts);
[X, Y, Z]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)), squeeze(F(3,:,:)));
subplot (1,2,1)
surf (X, Y, Z, eu)
title ('Numerical solution'), axis tight
subplot (1,2,2)
surf (X, Y, Z, problem_data.uex (X,Y,Z))
title ('Exact solution'), axis tight

% Display errors of the computed solution in the L2 and H1 norm
[error_h1, error_l2] = ...
           sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex)

%!demo
%! ex_laplace_beltrami

%!test
%! problem_data.geo_name = 'geo_roof.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4];
%! problem_data.c_diff  = @(x, y, z) ones(size(x));
%! alp = 5; bet = 1; L = 1;
%! g1 = @(x,y) (1-cos(atan(x./y))).*(1-sin(atan(x./y)));
%! g2 = @(x,y) cos(atan(x./y)) + sin(atan(x./y)) - 4*cos(atan(x./y)).*sin(atan(x./y));
%! gz = @(z) sin(pi*alp/L*z);
%! dg1 = @(x,y) sin(atan(y./x)).*(1-sin(atan(y./x))) - cos(atan(y./x)).*(1-cos(atan(y./x)));
%! problem_data.f = @(x, y, z) bet*((alp*pi/L)^2 * g1(x,y) - g2(x,y)).*gz(z);
%! problem_data.h = @(x, y, z, ind) bet * g1(x,y).*gz(z);
%! problem_data.uex     = @(x, y, z) bet * g1(x,y).*gz(z);
%! problem_data.graduex = @(x, y, z) cat (1, ...
%!                 reshape (-1./sqrt(x.^2+y.^2).*sin(atan(y./x))*bet.*gz(z).*dg1(x,y), [1, size(x)]), ...
%!                 reshape ( 1./sqrt(x.^2+y.^2).*cos(atan(y./x))*bet.*gz(z).*dg1(x,y), [1, size(x)]), ...
%!                 reshape ( alp*pi/L * cos(alp*pi/L*z) .* bet .* g1(x,y), [1, size(x)]));
%! method_data.degree     = [3 3];       % Degree of the splines
%! method_data.regularity = [2 2];       % Regularity of the splines
%! method_data.nsub       = [9 9];       % Number of subdivisions
%! method_data.nquad      = [4 4];       % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_laplace (problem_data, method_data);
%! [error_h1, error_l2] = ...
%!            sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex);
%! assert (msh.nel, 81)
%! assert (space.ndof, 144)
%! assert (error_h1, 0.0407258650845198, 1e-15)
%! assert (error_l2, 0.00101798638120372, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function christmas

% Anselm Ivanovas, anselm.ivanovas@student.unisg.ch

%Basically just a nice plot for some christmas fun.
%3D Plot of a hhristmas tree with some presents and snow

%% setup
snow=800;     % number of snow flakes [0 .. 5000]

%% draw tree
h=0:0.2:25; %vertical grid
[X,Y,Z] = cylinder(tree(h)); %produce a tree formed cylinder
Z=Z*25; %scale to the right heigth

%Add some diffusion to the surface of the tree to make it look more real

treeDiffusion=rand(126,21)-0.5;%some horizontal diffusion data

%add diffusion to the grid points
for cnt1=1:21

    for cnt2=16:126%starting above the trunk
        %get the angle to always diffuse in direction of the radius
        angle=atan(Y(cnt2,cnt1)/X(cnt2,cnt1));
        %split the diffusion in the two coordinates, depending on the angle
        X(cnt2,cnt1)=X(cnt2,cnt1)+cos(angle)*treeDiffusion(cnt2,cnt1);
        Y(cnt2,cnt1)=Y(cnt2,cnt1)+sin(angle)*treeDiffusion(cnt2,cnt1);
        %some Vertical diffusion for each point
        Z(cnt2,cnt1)=Z(cnt2,cnt1)+(rand-0.5)*0.5;
    end

end
%draw the tree
surfl(X,Y,Z,'light')

%% View and format

%Use as nice green color map (darker at the bottom, lighter at the top)
r=(0.0430:(0.2061/50):0.2491)';%red component
g=(0.2969:(0.4012/50):0.6981)';%green component
b=(0.0625:(0.2696/50):0.3321)';%blue component
map=[r,g,b];%join in a map
for cnt=1:6
    %change the lower part to brown for the trunk
    map(cnt,:)=[77,63,5]/265;
end

colormap(map)%set the map
view([-37.5,4])%Change the view to see a little more of the Actual 3D tree
lighting phong %some nice lighting
shading interp %remove grid and smoothen the surface color
axis equal %takes care of display in the right proportion
axis([-10 10 -10 10 0 30]) %give some more axis space (for the snow later)
axis off %but don't show axis
hold on %to draw the rest
title('HAPPY HOLIDAYS')%self explaining

%% Presents
%Draw some presents around the tree (each with random color)
drawPresent(2,-4,0,3,3,2);
drawPresent(-4,3,0,2,3,1.5);
drawPresent(5,3,0,4,3,3);
drawPresent(-14,-5,0,6,3,1);
drawPresent(-9,-10,0,2,2,2);
drawPresent(0,4,0,4,3,3);
drawPresent(-6,-13,0,3,3,3);

%% Snow

%create some random 3D coordinates for the snow (amount as in setup above)
snowX=(rand(snow,1)*25-12.5);
snowY=(rand(snow,1)*25-12.5);
snowZ=(rand(snow,1)*27);
%Note:Some flakes will end up IN the tree but just can't be seen then
plot3(snowX,snowY,snowZ,'w*')%plot coordinates as white snow flakes
hold off%Done
end % of function

%% ============= private functions

function r=tree(h)%Gives a profile for the tree
for cnt=1:length(h)

    if(h(cnt)==0)%no Width at the bottom. Ensures a "closed" trunk
        r(cnt)=0;
    end
    %smaller radius for the trunk
    if (h(cnt)>0 && h(cnt)<=3)
        r(cnt)=1.5;
    end

    %reduce radius gradually from 8 to 0. Note: will only work with a trunk heigth
    %of 3 and a whole tree heigth of 25. Scale the height of the tree in
    %the "draw tree" section, since the cylinder command will return a 1
    %unit high cylinder anyway
    if(h(cnt)>3)
        r(cnt)=8-(h(cnt)-3)*0.3636;
    end

end

end % of function

%Draws a present with the given coordinate + size in a random color
%Note:Given coordinates apply to the lower front + left corner of the
%present (the one closest to the viewer) as seen in the plot
function drawPresent(dx,dy,dz,scalex,scaley,scalez)

%the standard present coordinates
presentX=[0.5 0.5 0.5 0.5 0.5; 0 1 1 0 0; 0 1 1 0 0; 0 1 1 0 0; 0.5 0.5 0.5 0.5 0.5];
presentY=[0.5 0.5 0.5 0.5 0.5; 0 0 1 1 0; 0 0 1 1 0; 0 0 1 1 0; 0.5 0.5 0.5 0.5 0.5];
presentZ=[0 0 0 0 0; 0 0 0 0 0; 0.5 0.5 0.5 0.5 0.5; 1 1 1 1 1; 1 1 1 1 1];

%draw some presents with random colors
%scale present and move it to the right place and get the plot handle
myHandle=surf((presentX*scalex+dx),(presentY*scaley+dy), (presentZ*scalez+dz));
%some random color map
randColorMap(:,:,1)=repmat(rand,[5,5]);%r component
randColorMap(:,:,2)=repmat(rand,[5,5]);%g component
randColorMap(:,:,3)=repmat(rand,[5,5]);%b component
%Assign colormap just to the plot handle object of the present, so the tree
%does not change color
set(myHandle,'CData',randColorMap)
shading interp %Nice shding + without grid

end % of function


text:
classdef SVM < handle

  properties

    svm;

  end

  methods

    function a = SVM()
      %Wrapper on libsvm.

    end

    function train(a, inputs, outputs)

      a.svm = svmtrain(outputs, inputs);

    end

    function outputs = test(a, inputs)

      outputs = svmpredict(ones(size(inputs, 1), 1), inputs, a.svm);

    end

    function margins = margins(a, inputs)

      [guess acc margins] = svmpredict(ones(size(inputs, 1), 1), inputs, a.svm);

    end

    function [weights theta] = weights(a)
        weights = (a.svm.sv_coef' * full(a.svm.SVs));
        theta = a.svm.rho;
    end

  end

end

author: epfl-lasa


text:
function [coregisteredImage, affineCoregistrationGeometry, outputOtherImages] =...
    coregister_to(this, stationaryImage, varargin)
% Coregister this MrImage to another MrImage using SPM's coregister.
%
%   Y = MrImage()
%   cY = Y.coregister_to(stationaryImage)
%
% This is a method of class MrImage.
%
% IN
%   stationaryImage     cell {nImages,n} of MrImages or single MrImage that
%                       serves as stationary" or reference image to which
%                       this image is coregistered to
%
%  optional parameter name/value pairs:
%
%  transformation estimation and application
%       applyTransformation
%                   'geometry'      MrImageGeometry is updated,
%                                   MrImage.data remains untouched
%                   'data'          MrImage.data is resliced to new
%                                   geometry
%                                   NOTE: An existing
%                                   transformation in MrImageGeometry will
%                                   also be applied to MrImage, combined
%                                   with the calculated one for
%                                   coregistration
%                   'none'          transformation matrix is
%                                   computed, but not applied to geometry
%                                   of data of this image
%       trafoParameters             'translation', 'rigid', 'affine',
%                                   'rigidscaled' or
%                                   [1,1-12] vector of starting parameters
%                                   for transformation estimation number of
%                                   elements decides whether
%                                   translation only (1-3)
%                                   rigid (4-6)
%                                   rigid and scaling (7-9)
%                                   affine (10-12)
%                                   transformation is performed
%                                   default: 'rigid' (as in SPM)
% SPM input parameters:
%          separation           optimisation sampling steps (mm)
%                               default: [4 2]
%          objectiveFunction    cost function string:
%                               'mi'  - Mutual Information
%                               'nmi' - Normalised Mutual Information
%                               'ecc' - Entropy Correlation Coefficient
%                               'ncc' - Normalised Cross Correlation
%                               default: 'nmi'
%          tolerances           tolerances for accuracy of each param
%                               default: [0.02 0.02 0.02 0.001 0.001 0.001]
%          histSmoothingFwhm    smoothing to apply to 256x256 joint histogram
%                               default: [7 7]
%          otherImages          cell(nImages,1) of other images (either
%                               file names or MrImages) that should undergo
%                               the same trafo as this images due to coreg
%                               default: {}
%           doPlot              set to true for graphical output and PS file creation
%                               in SPM graphics window
%                               default: false
%
%   Parameters for high-dim application:
%       representationIndexArray:   either an MrImageObject or a selection
%                                   (e.g. {'echo', 1} which is then applied to
%                                   obtain one 4D image
%                                   default representationIndexArray: first
%                                   index of all extra (non-4D) dimensions
%       applicationIndexArray:      a selection which defines one or multiple
%                                   4D images on which the estimated parameters
%                                   are applied
%                                   default applicationIndexArray: all non-4D
%                                   dimensions
%       splitComplex                'ri' or 'mp'
%                                   If the data are complex numbers, real and imaginary or
%                                   magnitude and phase are realigned separately.
%                                   default: ri (real and imaginary)
%
% OUT
%
% EXAMPLE
%   coregister_to
%
%   See also MrImageSpm4D.coregister_to

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-11-28
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% SPM parameters
defaults.doPlot = false; % set to true for graphical output and PS file creation
defaults.otherImages = {};
defaults.objectiveFunction = 'nmi';
defaults.separation = [4 2 1];
defaults.tolerances = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
defaults.histSmoothingFwhm = [7 7];
defaults.trafoParameters = 'rigid';
% other coreg_to parameters of MrImageSpm4D
defaults.applyTransformation = 'data';
defaults.otherImages = {};
args = tapas_uniqc_propval(varargin, defaults);
% tapas_uniqc_strip_fields(args);

methodParameters = {args};

%% Distinguish cases

%% 1) this: 3D; stationary: 3D; other: nD
%  -> trivial case, directly mimicking SPM's coregister and passing to
%  MrImageSpm4D

thisNonSDims = this.dimInfo.get_non_singleton_dimensions;
thisIs3D = numel(thisNonSDims) == 3;
thisIsReal = isreal(this);
thisIsReal3D = thisIs3D && thisIsReal;

stationaryNonSDims = stationaryImage.dimInfo.get_non_singleton_dimensions;
stationaryIs3D = numel(stationaryNonSDims) == 3;
stationaryIsReal = isreal(stationaryImage);
stationaryIsReal3D = stationaryIs3D && stationaryIsReal;

if thisIsReal3D &&  stationaryIsReal3D % just do as SPM does!

    % for now, the stationary is 3D
    stationaryImage = stationaryImage.remove_dims();

    xyzDimLabels = {'x','y','z'};
    splitDimIndexArray = setdiff(1:this.dimInfo.nDims, ...
        this.dimInfo.get_dim_index(xyzDimLabels));
    splitDimLabels = this.dimInfo.dimLabels(splitDimIndexArray);

    [coregisteredImage, affineCoregistrationGeometry, outputOtherImages] = ...
        this.apply_spm_method_per_4d_split(...
        @(x, y) coregister_to(x, stationaryImage, y), ...
        'methodParameters', methodParameters, 'splitDimLabels', splitDimLabels);

end

end
%% 2) this: nD; stationary: 3D; other: cell(nSplits,nOtherImages) of 3D images
%% 2a) representation: 3D; application: nD
%  -> one 3D part represents this for the coregistration to stationary, and
%  the estimated coreg parameters are then applied to all images in the
%  application selections.
%% 2b) representation: nD; application: nD
%  -> each 3D part of "this" will be individually coregistered to stationary;
%     and kth coreg is applied to coresponding k-th cell of other images
%     Note: cell of nSplits of 3D images can be created by
%     other.split('splitDimLabels', {'t','echo'}) or similar
%% 3) this: 3D; stationary: cell(nStationary,1) of 3D images; other: cell (nStationary, nOtherImages)
%  -> "this" is coregistered to each of the stationaries (therefore 3+x D
%  output), and parameters of iStationary's coeregistration are applied to corresponding other images in cell {iStationary,:}
%% 4) this: nD; stationary: cell(nStationary,1) of 3D images; other cell(nStationary, nOtherImages)
%     Note: nStationary == nSplits of "this" here
%  -> Each 3D part of "this" is coregistered to the corresponding iSplit
%  (iStationary) entry of cell stationary, and the coreg parameters are
%  applied to all entries in cell {iStationary,:} of other images.

author: translationalneuromodeling

repo: tapas


text:
function [L,U,P] = lu_rightp (A)
%LU_RIGHTP right-looking LU factorization, with partial pivoting.
%
% Example:
%   [L,U,P] = lu_rightp (A)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

n = size (A,1) ;
P = eye (n) ;
for k = 1:n
    [x,i] = max (abs (A (k:n,k))) ;                           % partial pivoting
    i = i+k-1 ;
    P ([k i],:) = P ([i k], :) ;
    A ([k i],:) = A ([i k], :) ;                              % (6.10), (6.11)
    A (k+1:n,k) = A (k+1:n,k) / A (k,k) ;                               % (6.12)
    A (k+1:n,k+1:n) = A (k+1:n,k+1:n) - A (k+1:n,k) * A (k,k+1:n) ;     % (6.9)
end
L = tril (A,-1) + eye (n) ;
U = triu (A) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Compile MatConvNet
% TODO

%% Setup classifier
setup_path = 'D:\Users\sim0629\matconvnet-1.0-beta20\matlab';
% need to be called
run(fullfile(setup_path, 'vl_setupnn.m'));
load('imdb.mat');
load('net-epoch-57.mat');

%% Detection and Classification
close all;
[filename, pathname, ~] = uigetfile({'*.jpg;*.tif;*.png;*.gif','All Image Files'});
im = imread(fullfile(pathname, filename));

[ellipses, resized_im] = detectEllipses(im, true);
num_of_ellipses = length(ellipses);
features = zeros(50 * num_of_ellipses, 100, 3, 'uint8');
scoreses = zeros(num_of_ellipses, 8);

for i = 1 : num_of_ellipses
    feature = extractFeatureImage(resized_im, ellipses, i, false);
    scoreses(num_of_ellipses+1-i,:) = classifyFeature(feature, images, net)';
    features((num_of_ellipses-i)*50+1 : (num_of_ellipses+1-i)*50, :, :) = feature(52:101,1:100,:);
end

figure;

subplot('position', [0, 0.3, 0.3, 0.7]);
imshow(features);

subplot('position', [0.3, 0, 0.7, 0.3]);
imshow('labels.jpg');

subplot('position', [0.3, 0.3, 0.7, 0.7]);
colormap('jet');
imagesc(scoreses);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef globalICP < handle
    % GLOBALICP Class for global ICP problems.

    properties

        % Corresponding points
        CP

        % Point clouds
        PC

        % Path to folder for output data
        OutputFolder

        % Path to folder for temporary data
        TempFolder

        % Various data
        D

    end

    methods

        function obj = globalICP(varargin)
        % GLOBALICP Constructor method for global ICP class.
        % ----------------------------------------------------------------------
        % DESCRIPTION/NOTES
        % With the globalICP class the alignment of two or more point clouds can
        % be refined. A prerequisite for this is an approximate alignment of the
        % point clouds.
        % ----------------------------------------------------------------------
        % INPUT
        % 1 ['OutputFolder', OutputFolder]
        %   Path to directory in which output files are stored. If this option
        %   is omitted, the path given by the command 'cd' is used as directory.
        %
        % 2 ['TempFolder', TempFolder]
        %   Folder in which temporary files are saved, e.g. imported point
        %   clouds. If this option is omitted, the path given by the command
        %   'tempdir' is used as directory.
        % ----------------------------------------------------------------------
        % OUTPUT
        % [obj]
        % Object instance of class globalICP.
        % ----------------------------------------------------------------------
        % EXAMPLES
        % 1 Minimal working example with 6 point clouds.
        %
        %   % Create globalICP object
        %   icp = globalICP('OutputFolder', 'D:\temp');
        %
        %   % Add point clouds to object from plain text files
        %   % (Added point clouds are saved as mat files, e.g. LionScan1Approx.mat)
        %   icp.addPC('LionScan1Approx.xyz');
        %   icp.addPC('LionScan2Approx.xyz');
        %   icp.addPC('LionScan3Approx.xyz');
        %   icp.addPC('LionScan4Approx.xyz');
        %   icp.addPC('LionScan5Approx.xyz');
        %   icp.addPC('LionScan6Approx.xyz');
        %
        %   % Plot all point clouds BEFORE ICP (each in a different random color)
        %   figure; icp.plot('Color', 'random');
        %   title('BEFORE ICP'); view(0,0);
        %
        %   % Run ICP!
        %   icp.runICP('PlaneSearchRadius', 2);
        %
        %   % Plot all point clouds AFTER ICP
        %   figure; icp.plot('Color', 'random');
        %   title('AFTER ICP'); view(0,0);
        %
        % 2 Continued: Demo of the methods 'loadPC' and 'exportPC'.
        %
        %   % Load fifth point cloud to workspace and plot
        %   pc = icp.loadPC(5);
        %   figure; pc.plot;
        %
        %   % Export final point clouds
        %   icp.exportPC(1, 'LionScan1.xyz');
        %   icp.exportPC(2, 'LionScan2.xyz');
        %   icp.exportPC(3, 'LionScan3.xyz');
        %   icp.exportPC(4, 'LionScan4.xyz');
        %   icp.exportPC(5, 'LionScan5.xyz');
        %   icp.exportPC(6, 'LionScan6.xyz');
        % ----------------------------------------------------------------------
        % philipp.glira@gmail.com
        % ----------------------------------------------------------------------

        % Input parsing --------------------------------------------------------

        p = inputParser;
        p.addParameter('OutputFolder', cd     , @ischar);
        p.addParameter('TempFolder'  , tempdir, @ischar);
        p.parse(varargin{:});
        p = p.Results;

        % Directories ----------------------------------------------------------

        if ~exist(p.OutputFolder, 'dir'), mkdir(p.OutputFolder), end
        obj.OutputFolder = p.OutputFolder;

        if ~exist(p.TempFolder, 'dir'), mkdir(p.TempFolder), end
        obj.TempFolder = p.TempFolder;

        end

    end

end

author: pglira


text:
function varargout = rank(varargin)
%RANK      Rank of a SPHEREFUN.
%   RANK(F) produces an estimate of the rank of the approximant F.
%
%   RANK(F, TOL) is the number of singular values of F greater than TOL/N, where
%   N is the first singular value of F.
%
% See also SPHEREFUN/LENGTH.

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

[varargout{1:nargout}] = rank@separableApprox(varargin{:});

end


text:
function Population = EnvironmentalSelection(Population,N)
% The environmental selection of CMOPSO

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,N);
    Next = false(1,length(FrontNo));
    Next(FrontNo<MaxFNo) = true;

    PopObj = Population.objs;
    fmax   = max(PopObj(FrontNo==1,:),[],1);
    fmin   = min(PopObj(FrontNo==1,:),[],1);
    PopObj = (PopObj-repmat(fmin,size(PopObj,1),1))./repmat(fmax-fmin,size(PopObj,1),1);

    %% Select the solutions in the last front
    Last = find(FrontNo==MaxFNo);
    del  = Truncation(PopObj(Last,:),length(Last)-N+sum(Next));
    Next(Last(~del)) = true;
    % Population for next generation
    Population = Population(Next);
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    N = size(PopObj,1);

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,N);
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ f,qualMeasOut] = PCSD(proj,geo,angles,maxiter,varargin)
%PCSD solves the reconstruction problem using projection-controlled steepest descent method
%
%   PCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using
%   the projection data PROJ taken over ALPHA angles, corresponding to the
%   geometry described in GEO, using NITER iterations.
%
%   PCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter for the SART iterations.
%                  Default is 1
%
%   'lambdared':   Reduction of lambda.Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%       'init':    Describes diferent initialization techniques.
%                   •  'none'     : Initializes the image to zeros (default)
%                   •  'FDK'      : intializes image to FDK reconstrucition
%
%   'TViter':      Defines the amount of TV iterations performed per SART
%                  iteration. Default is 20
%
%   'maxL2err'     Maximum L2 error to accept an image as valid. This
%                  parameter is crucial for the algorithm, determines at
%                  what point an image should not be updated further.
%                  Default is 20% of the FDK L2 norm.
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri and Manasavee Lohvithee
%--------------------------------------------------------------------------

%% parse inputs
[beta,beta_red,f,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);

measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<2 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

% does detector rotation exists?
if ~isfield(geo,'rotDetector')
    geo.rotDetector=[0;0;0];
end

%% Create weigthing matrices for the SART step
% the reason we do this, instead of calling the SART fucntion is not to
% recompute the weigths every AwASD-POCS iteration, thus effectively doubling
% the computational time
% Projection weigth, W

geoaux=geo;
geoaux.sVoxel([1 2])=geo.sVoxel([1 2])*1.1; % a Bit bigger, to avoid numerical division by zero (small number)
geoaux.sVoxel(3)=max(geo.sDetector(2),geo.sVoxel(3)); % make sure lines are not cropped. One is for when image is bigger than detector and viceversa
geoaux.nVoxel=[2,2,2]'; % accurate enough?
geoaux.dVoxel=geoaux.sVoxel./geoaux.nVoxel;
W=Ax(ones(geoaux.nVoxel','single'),geoaux,angles,'Siddon','gpuids',gpuids);
W(W<min(geo.dVoxel)/4)=Inf;
W=1./W;

% Compute V
V=computeV(geo,angles,num2cell(angles),num2cell(1:length(angles)),'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end

%Initialize image.
%f=zeros(geo.nVoxel','single');

iter=0;
offOrigin=geo.offOrigin;
offDetector=geo.offDetector;
rotDetector=geo.rotDetector;
stop_criteria=0;
DSD=geo.DSD;
DSO=geo.DSO;
%%
while ~stop_criteria %POCS
    % If quality is going to be measured, then we need to save previous image
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = f; % only store if necesary
    end
    if (iter==0 && verbose==1);tic;end
    iter=iter+1;

    %Estimation error in the projection domain
    est_proj=Ax(f,geo,angles,'interpolated','gpuids',gpuids);
    delta_p=im3Dnorm(est_proj-proj,'L2');

    %Enforcing ART along all projections if squared delta_p > epsilon
    if (delta_p^2)>epsilon
        for jj=1:size(angles,2)
            if size(offOrigin,2)==size(angles,2)
                geo.offOrigin=offOrigin(:,jj);
            end
            if size(offDetector,2)==size(angles,2)
                geo.offDetector=offDetector(:,jj);
            end
            if size(rotDetector,2)==size(angles,2)
                geo.rotDetector=rotDetector(:,jj);
            end
            if size(DSD,2)==size(angles,2)
                geo.DSD=DSD(jj);
            end
            if size(DSO,2)==size(angles,2)
                geo.DSO=DSO(jj);
            end
            f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,jj).*(proj(:,:,jj)-Ax(f,geo,angles(:,jj),'gpuids',gpuids)),geo,angles(:,jj),'gpuids',gpuids));

        end
    end

    %Non-negativity projection on all pixels
    if nonneg
        f=max(f,0);
    end

    geo.offDetector=offDetector;
    geo.offOrigin=offOrigin;
    geo.DSD=DSD;
    geo.DSO=DSO;
    geo.rotDetector=rotDetector;
    if measurequality
        qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts);
    end

    % Compute L2 error of actual image. Ax-b
    dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2');
    % Compute change in the image after last SART iteration
    dp_vec=(f-f0);

    if iter==1
        step=1;
    else
        step=delta_p/delta_p_first;
    end
    f0=f;
    %  TV MINIMIZATION
    % =========================================================================
    %  Call GPU to minimize TV
    f=minimizeTV(f0,step,ng,'gpuids',gpuids);    %   This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays
    %                                             for ii=1:ng
    %                                                 %delta=-0.00038 for thorax phantom
    %                                                 df=weighted_gradientTVnorm(f,delta);
    %                                                 df=df./im3Dnorm(df,'L2');
    %                                                 f=f-(step.*df);
    %                                             end

    % Compute change by TV min
    dg_vec=(f-f0);

    if iter==1
        delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2');
    end

    % Reduce SART step
    beta=beta*beta_red;

    % Check convergence criteria
    % ==========================================================================

    %Define c_alpha as in equation 21 in the journal
    c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2));
    %This c is examined to see if it is close to -1.0

    if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter
        if verbose
            disp(['Stopping criteria met']);
            disp(['   c    = ' num2str(c)]);
            disp(['   beta = ' num2str(beta)]);
            disp(['   iter = ' num2str(iter)]);
        end
        stop_criteria=true;
    end

    if (iter==1 && verbose==1)
        expected_time=toc*maxiter;
        disp('PCSD');
        disp(['Expected duration  :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end

end
end

function [beta,beta_red,f0,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin)
opts=     {'lambda','lambda_red','init','tviter','verbose','maxl2err','qualmeas','nonneg','gpuids','redundancy_weighting'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('CBCT:PCSD:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('CBCT:PCSD:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option isnot default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('CBCT:PCSD:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end
    % parse inputs
    switch opt
        % Verbose
        %  =========================================================================
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
        % Lambda
        %  =========================================================================
        case 'lambda'
            if default
                beta=1;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:PCSD:InvalidInput','Invalid lambda')
                end
                beta=val;
            end
        % Lambda reduction
        %  =========================================================================
        case 'lambda_red'
            if default
                beta_red=0.99;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:PCSD:InvalidInput','Invalid lambda')
                end
                beta_red=val;
            end
        % Initial image
        %  =========================================================================
        case 'init'
            if default || strcmp(val,'none')
                f0=zeros(geo.nVoxel','single');

            else
                if strcmp(val,'FDK')
                    f0=FDK(proj, geo, angles);
                else
                    error('TIGRE:PCSD:InvalidInput','Invalid init')
                end
            end
        % Number of iterations of TV
        %  =========================================================================
        case 'tviter'
            if default
                ng=20;
            else
                ng=val;
            end
        %  Maximum L2 error to have a "good image"
        %  =========================================================================
        case 'maxl2err'
            if default
                epsilon=im3Dnorm(FDK(proj,geo,angles))*0.2; %heuristic
            else
                epsilon=val;
            end
        %Image Quality Measure
        %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:PCSD:InvalidInput','Invalid quality measurement parameters');
                end
            end
        %  Non negative
        %  =========================================================================
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
        %  GPU Ids
        %  =========================================================================
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        otherwise
            error('TIGRE:PCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in PCSD()']);

    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
 function dot = ir_dot_double(a, b)
%function dot = ir_dot_double(a, b)
%| compute dot product between vectors a and b in double precision
%| dot = sum(a(:) .* b(:), 'double'); % double accumulate
%| user must apply conj() to a or b before calling for complex case

if nargin < 2, ir_usage, end

dot = sum(a(:) .* b(:), 'double'); % double accumulate


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = putpayoff(S,Strike)
% payoff of a European put option
    y =  max(Strike - S, 0);
end


text:
% Test file for trigtech/prolong.m

function pass = test_prolong(pref)

if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

F = @(x) exp(sin(pi*x));
f = testclass.make(F, [], pref);

% Odd prolongation.
k = 101;
go = prolong(f, k);
x = testclass.trigpts(k);
pass(1) = size(go,1) == k  && norm(go.values - F(x), inf) < 10*vscale(go)*eps;

% Even prolongation.
k = 100;
ge = prolong(f, k);
x = testclass.trigpts(k);
pass(2) = size(ge,1) == k  && norm(ge.values - F(x), inf) < 10*vscale(ge)*eps;

% Odd restriction from odd length f.
k = 89;
g = prolong(go, k);
x = testclass.trigpts(k);
pass(3) = size(g,1) == k  && norm(g.values - F(x), inf) < 10*vscale(g)*eps;

% Even restriction from odd length f.
k = 88;
g = prolong(go, k);
x = testclass.trigpts(k);
pass(4) = size(g,1) == k  && norm(g.values - F(x), inf) < 10*vscale(g)*eps;

% Odd restriction from even length f.
k = 89;
g = prolong(ge, k);
x = testclass.trigpts(k);
pass(5) = size(g,1) == k  && norm(g.values - F(x), inf) < 10*vscale(g)*eps;

% Even restriction from even length f.
k = 88;
g = prolong(ge, k);
x = testclass.trigpts(k);
pass(6) = size(g,1) == k  && norm(g.values - F(x), inf) < 10*vscale(g)*eps;

% Restriction to length 1.
g = prolong(f, 1);
pass(7) = size(g,1) == 1;

% Array-valued tests.
F = @(x) exp(sin(pi*x));
f = testclass.make(@(x) [F(x), -F(x)], [], pref);
k = 101;
g = prolong(f, k);
x = testclass.trigpts(k);
values = g.coeffs2vals(g.coeffs);
pass(8) = size(g,1) == k && norm(values - [F(x), -F(x)], inf) < ...
    100*max(vscale(g)*eps);

g = prolong(f, 1);
pass(9) = size(g,1) == 1;

% Same length test
g = prolong(f, size(f,1));
fvalues = f.coeffs2vals(f.coeffs);
gvalues = g.coeffs2vals(g.coeffs);
pass(10) = all(fvalues(:) == gvalues(:));

% Values test.
v = [1 2 3];
f = testclass.make(v, [], pref);
g = prolong(f, 5);
values = g.coeffs2vals(g.coeffs);
pass(11) = norm(values - repmat([1 2 3], 5, 1), inf) < ...
   10*max(vscale(g)*eps);

end