text:
% [readframe,nframes,fid,headerinfo] = get_readframe_fcn(filename)
% get_readframe_fcn inputs the name of a video, opens this video, and creates a
% function that can be used for random access to frames in the video. it can
% input files of type fmf, sbfmf, and avi. if videoio is installed on your
% machine, then get_readframe_fcn will use videoio to read all types of avi
% supported by videoio. otherwise, it uses VideoReader to open any type of avi
% supported by VideoReader. Run "help VideoReader" to determine what types of avi
% files are readable through VideoReader on your OS.
%
% Input:
% filename is the name of the video to read. The type of video is determined
% based on the extension, so be sure the extension of your movie is correct.
%
% Outputs:
% readframe is the handle to a function which inputs a frame number, reads this
% frame from filename, and outputs the corresponding frame.
% nframes is the number of frames in the video. For avis, this may be
% approximate, as it is computed from the duration and framerate.
% fid is the file identifier for filename. get_readframe_fcn will open the file
% filename. you must close fid yourself when you are done using the returned
% readframe.
% headerinfo is a struct with fields depending on the type of movie

function [readframe,nframes,fid,headerinfo] = get_readframe_fcn(filename,varargin)

i = find(strcmpi(varargin(1:2:end),'preload'));
if ~isempty(i),
  preload = varargin{i+1};
  varargin(i:i+1) = [];
else
  preload = false;
end
didpreload = false;

% allow videoio library to be used if it is installed and on the path
%CTRAX_ISVIDEOIO = exist('videoReader','file');
CTRAX_ISVIDEOIO = false;

if iscell(filename),

  readframes = cell(size(filename));
  nframes = inf;
  fid = nan(size(filename));
  headerinfo = cell(size(filename));
  for i = 1:numel(filename),
    [readframes{i},nframescurr,fid(i),headerinfo{i}] = get_readframe_fcn(filename{i},varargin{:});
    nframes = min(nframes,nframescurr);
  end

  readframe = @(f) multi_read_frame(f,readframes);
  return;

end

[~,ext] = splitext(filename);

if ispc && ~exist(filename,'file') && ~strcmpi(ext,'.seq'),
  [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
  if didfind,
    filename = actualfilename;
  end
end

if strcmpi(ext,'.fmf'),
  [header_size,version,nr,nc,bytes_per_chunk,nframes,data_format] = fmf_read_header(filename);
  fid = fopen(filename);
  readframe = @(f) fmfreadframe(fid,header_size+(f-1)*bytes_per_chunk,nr,nc,bytes_per_chunk,data_format);
  headerinfo = struct('header_size',header_size,'version',version,'nr',nr,'nc',nc,...
    'bytes_per_chunk',bytes_per_chunk,'nframes',nframes,'data_format',data_format,'type','fmf');
elseif strcmpi(ext,'.sbfmf'),
  [nr,nc,nframes,bgcenter,bgstd,frame2file] = sbfmf_read_header(filename);
  fid = fopen(filename);
  readframe = @(f) sbfmfreadframe(f,fid,frame2file,bgcenter);
  headerinfo = struct('nr',nr,'nc',nc,'nframes',nframes,'bgcenter',bgcenter,...
    'bgstd',bgstd,'frame2file',frame2file,'type','sbfmf');
elseif strcmpi(ext,'.ufmf'),
  headerinfo = ufmf_read_header(filename);
  readframe = ufmf_get_readframe_fcn(headerinfo,varargin{:});
  nframes = headerinfo.nframes;
  fid = headerinfo.fid;
elseif strcmpi(ext,'.mmf'),
  headerinfo = mmf_read_header(filename);
  readframe = mmf_get_readframe_fcn(headerinfo,varargin{:});
  nframes = headerinfo.nframes;
  fid = headerinfo.fid;
elseif strcmpi(ext,'.tif') || strcmp(ext,'.tiff'),
  info = imfinfo(filename);
  isimseq = false;
  if numel(info) == 1,
    filespec = regexprep(filename,['_\d+\.',ext,'$'],['_*.',ext]);
    imfiles = mydir(filespec);
    if numel(imfiles) > 1,
      imfiles = sort(imfiles);
      im = imread(imfiles{1});
      headerinfo = struct('nr',size(im,1),'nc',size(im,2),'ncolors',size(im,3),'nframes',numel(imfiles),...
        'type','imseq','imfiles',{imfiles});
      readframe = @(f) imseq_read_frame(f,imfiles);
      nframes = headerinfo.nframes;
      fid = -1;
      isimseq = true;
    end
  end
  if ~isimseq,
    headerinfo = struct('nr',info(1).Height,'nc',info(1).Width,'nframes',numel(info),'type','tif',...
      'bitdepth',info(1).BitDepth);
    readframe = @(f) deal(imread(filename,f),f);
    nframes = headerinfo.nframes;
    fid = -1;
  end
elseif strcmpi(ext,'.png'),
  info = imfinfo(filename);
  isimseq = false;
  filespec = regexprep(filename,'_\d+\.png$','_*.png');
  imfiles = mydir(filespec);
  if numel(imfiles) > 1,
    imfiles = sort(imfiles);
    im = imread(imfiles{1});
    headerinfo = struct('nr',size(im,1),'nc',size(im,2),'ncolors',size(im,3),'nframes',numel(imfiles),...
      'type','imseq','imfiles',{imfiles});
    readframe = @(f) imseq_read_frame(f,imfiles);
    nframes = headerinfo.nframes;
    fid = -1;
    isimseq = true;
  end
elseif strcmpi(ext,'.klb'),
  tfKLBLib = exist('readKLBslice','file')>0;
  if ~tfKLBLib
    error('get_readframe_fcn:klb','Cannot find readKLBslice function. Make sure KLB matlab library is installed.');
  end

  dim = myparse(varargin,'dim',3);
  headerinfo = readKLBheader(filename);
  readframe = klb_get_readframe_fcn(filename,varargin{:});
  nframes = headerinfo.xyzct(dim); % AL 201507: KLB doc a little unclear here but this seems right
  fid = 0;
elseif strcmpi(ext,'.mat'),
  videofiletype = load(filename,'videofiletype');
  switch videofiletype,

    case 'SingleLarvaTracker',
      videodata = load(filename);
      readframe = @(f) ReadSingleLarvaTrackerFrame(f,videodata.firstframeim,videodata.imraw,videodata.finalbbox,videodata.fps,varargin{:});
      nframes = numel(imraw);
      fid = 0;
      [nr,nc,~] = size(firstframeim);
      headerinfo = struct('nr',nr,'nc',nc,'nframes',nframes,'bgcenter',firstframeim,...
        'type','SingleLarvaTracker');
    otherwise
      error('Do not know how to parse mat file of type %s',videofiletype);
  end

elseif strcmpi(ext,'.seq'),
  [indexfilename,seqtype] = myparse(varargin,'indexfilename',0,'seqtype',0);

  if ischar(indexfilename),
    seqtype = 'shay';
  end

  % get file names
  if ispc && ~exist(filename,'file'),

    [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
    if ~didfind,
      error('Could not find movie file %s',filename);
    end
    filename0 = filename;
    % use actualfilename instead
    filename = actualfilename;

    % try the index file using the original filename
    if ~ischar(indexfilename),
      indexfilename = regexprep(filename0,'seq$','mat');
      [indexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
      if ~didfind && ~strcmp(filename0,filename),
        % try to find the index file using the source filename
        indexfilename = regexprep(filename,'seq$','mat');
        [indexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
        % if didn't find index file, then maybe this is a piotr seq file
        if ~didfind,
          if ischar(seqtype) && strcmpi(seqtype,'shay'),
            error('Could not find index file for %s',filename);
          else
            seqtype = 'piotr';
          end
        end
      end
    end

  else

    % set seqtype to piotr if indexfilename doesn't exist
    if ~ischar(seqtype) && ~ischar(indexfilename),
      indexfilename = regexprep(filename,'seq$','mat');
      if ~exist(indexfilename,'file'),
        seqtype = 'piotr';
      end
    end

  end

  % get actual filename for shortcuts
  if strcmpi(seqtype,'piotr'),
    warning('Closing file currently not implemented for Piotr''s seq files...');
    headerinfo = seqIo(filename,'getinfo');
    headerinfo.nr = headerinfo.height;
    headerinfo.nc = headerinfo.width;
    headerinfo.nframes = headerinfo.numFrames;
    headerinfo.type = 'seq_piotr';
    fid = 0;
    nframes = headerinfo.numFrames;
    sr = seqIo(filename,'r');
    readframe = @(f) seq_read_frame_piotr(f,sr);
  else

    if ischar(indexfilename),
      headerinfo = r_readseqinfo(filename,indexfilename);
    else
      headerinfo = r_readseqinfo(filename);
    end
    if numel(headerinfo.m_aiSeekPos) < headerinfo.m_iNumFrames,
      headerinfo.m_iNumFrames = numel(headerinfo.m_aiSeekPos);
    end
    headerinfo.nr = headerinfo.m_iHeight;
    headerinfo.nc = headerinfo.m_iWidth;
    headerinfo.nframes = headerinfo.m_iNumFrames;
    headerinfo.type = 'seq';
    fid = 0;
    nframes = headerinfo.m_iNumFrames;
    readframe = @(f) read_seq_frame(headerinfo,f);
  end

elseif strcmpi(ext,'.h5'),

  [dims,datasetname,fixdims,fixidx] = myparse(varargin,'slicedims',3,'datasetname',[],'fixdims',[],'fixidx',[]);
  if isempty(datasetname),
    headerinfo = h5info(filename);
    datasetname = headerinfo.Name;
    headerinfo = headerinfo.Datasets;
  else
    headerinfo = h5info(filename,datasetname);
  end
  headerinfo.datasetname = datasetname;
  headerinfo.dims = dims;
  headerinfo.fixdims = fixdims;
  headerinfo.fixidx = fixidx;
  headerinfo.sz = headerinfo.Dataspace.Size;
  headerinfo.ndims = numel(headerinfo.sz);
  headerinfo.nframes = headerinfo.sz(headerinfo.dims);

  readframe = h5_get_readframe_fcn(filename,headerinfo);
  nframes = headerinfo.nframes;
  fid = 0;

else
  fid = 0;

  isindexedmjpg = false;
  if strcmpi(ext,'.mjpg'),
    [indexfilename] = myparse(varargin,'indexfilename',[]);
    if isempty(indexfilename),
      % see if this looks like an indexed mjpg file
      [p,n] = fileparts(filename);
      indexfilename = fullfile(p,[n,'.txt']);
      if ~exist(indexfilename,'file') && ispc,
        [actualindexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
        if didfind,
          indexfilename = actualindexfilename;
        end
      end
      if exist(indexfilename,'file'),
        try
          ReadIndexedMJPGHeader(filename,indexfilename);
          isindexedmjpg = true;
        catch
          indexfilename = 0;
        end
      else
        indexfilename = 0;
      end
    elseif ischar(indexfilename),
      isindexedmjpg = true;
    end
  end

  if isindexedmjpg,
    % get file names
    if ispc && ~exist(filename,'file'),
      [actualfilename,didfind] = GetPCShortcutFileActualPath(filename);
      if ~didfind,
        error('Could not find movie file %s',filename);
      end
      filename0 = filename;
      % use actualfilename instead
      filename = actualfilename;

      [actualindexfilename,didfind] = GetPCShortcutFileActualPath(indexfilename);
      if ~didfind,
        error('Could not find index file %s',indexfilename);
      end
      indexfilename0 = indexfilename;
      % use actualfilename instead
      indexfilename = actualindexfilename;
    end

    headerinfo = ReadIndexedMJPGHeader(filename,indexfilename);
    headerinfo.fid = 0;
    nframes = headerinfo.nframes;
    readframe = @(f) read_mjpg_frame(headerinfo,f);

  else

    if CTRAX_ISVIDEOIO,
      readerobj = videoReader(filename,'preciseFrames',30,'frameTimeoutMS',5000);
      info = getinfo(readerobj);
      nframes = info.numFrames;
      seek(readerobj,0);
      seek(readerobj,1);
      readframe = @(f) videoioreadframe(readerobj,f);
      headerinfo = info;
      headerinfo.type = 'avi';
    else

      % optional argument to not use the read function of VideoReader
      [useRead,neednframes] = myparse(varargin,'useRead',true,'neednframes',true);

      try
        readerobj = VideoReader(filename);
        headerinfo = get(readerobj);

        headerinfo.type = 'avi';
        headerinfo.nr = headerinfo.Height;
        headerinfo.nc = headerinfo.Width;

        if preload,

          nframes = 0;
          preloaddata = struct;
          preloaddata.ims = {};
          preloaddata.ts = [];

          while true,

            if ~hasFrame(readerobj),
              break;
            end
            preloaddata.ts(end+1) = get(readerobj,'CurrentTime');
            preloaddata.ims{end+1} = readFrame(readerobj);
            nframes = nframes+1;

          end

          assert(numel(preloaddata.ims) == nframes);

        else
          if neednframes && useRead,
            nframes = get(readerobj,'NumberOfFrames');
          end
          if ~neednframes || ~useRead || isempty(nframes),
            % approximate nframes from duration
            nframes = round(get(readerobj,'Duration')*get(readerobj,'FrameRate'));
          end
          headerinfo.readerobj = readerobj;

        end

        headerinfo.nframes = nframes;
        frame_rate = headerinfo.FrameRate ;  % Hz
        dt = 1/frame_rate ;  % s
        headerinfo.timestamps = dt * (0:(nframes-1)) ;  % s

        if preload,

          readframe = @(f) preload_read_frame(preloaddata,f);
          didpreload = true;

        else

          if useRead,
            readframe = @(f) avi_read_frame(readerobj,headerinfo,f);
          else
            readframe = @(f) avi_read_frame_useReadFrame(readerobj,headerinfo,f);
          end

        end
      catch ME_videoreader
        error('Could not open file %s with VideoReader: %s',...
          filename,getReport(ME_videoreader));
      end
    end
  end
end

if preload && ~didpreload,

  preloaddata = struct;
  preloaddata.ims = cell(1,headerinfo.nframes);
  preloaddata.ts = nan(1,headerinfo.nframes);

  for f = 1:headerinfo.nframes,
    [preloaddata.ims{f},preloaddata.ts(f)] = readframe(f);
  end

  readframe = @(f) preload_read_frame(preloaddata,f);
  didpreload = true;

end

function varargout = multi_read_frame(f,readframes)

[im,timestamp] = readframes{1}(f);

for i = 2:numel(readframes),
  im = cat(2,im,readframes{i}(f));
end
varargout{1} = im;
if nargout >= 2,
  varargout{2} = timestamp;
end

function [im,timestamp] = seq_read_frame_piotr(f,sr)
im = [];
timestamp = [];

if ~sr.seek(f-1),
  warning('Could not seek to frame %d',f);
  return;
end

[im,timestamp] = sr.getframe();

function [im,timestamp] = avi_read_frame_useReadFrame(readerobj,headerinfo,f)

timestamp = (f-1)/headerinfo.FrameRate;
oldtimestamp = get(readerobj,'CurrentTime');
if oldtimestamp ~= timestamp,
  set(readerobj,'CurrentTime',timestamp);
end
timestamp = get(readerobj,'CurrentTime');
im = readFrame(readerobj);

function [im,timestamp] = preload_read_frame(preloaddata,f)

try
  im = preloaddata.ims{f};
  timestamp = preloaddata.ts(f);
catch ME,
  if ~isempty(preloaddata.ims),
    im = preloaddata.ims{end};
    timestamp = preloaddata.ts(end);
  end
  warning('Error reading preloaded frame %d: %s',f,getReport(ME));
end

function [im,timestamp] = avi_read_frame(readerobj,headerinfo,f)

try
  im = read(readerobj,f);
catch ME,
  warning('Error reading the first try: %s',getReport(ME));
  pause(.01);
  im = read(readerobj,f);
end
timestamp = (f-1)/headerinfo.FrameRate;

function [im,f] = imseq_read_frame(f,imfiles)

im = imread(imfiles{f});

author: kristinbranson

repo: JAABA


function [ ds ] = TrainValidateMarket(ds, learn_algs, X, X_test, X_query, idxa, idxb, flag)
%
% Input:
%
%  ds - data struct that stores the result
%  learn_algs - algorithms that are used for cross validation
%  X - input matrix, each column is an input vector [DxN*2]. N is the
%  number of pairs.
%  idxa - index of image A in X [1xN]
%  idxb - index of image B in X [1xN]
%
% Output:
%
%  ds - contains the result

c = 1;
clc;
% train
for aC=1:length(learn_algs)
    cHandle = learn_algs{aC};
    fprintf('    training %s ',upper(cHandle.type));
    if aC > 3
        s = learnPairwise(cHandle,X,idxa(1:10:end),idxb(1:10:end),logical(flag(1:10:end)));
    else
        s = learnPairwise(cHandle,X,idxa,idxb,logical(flag));
    end
    if ~isempty(fieldnames(s))
        fprintf('... done in %.4fs\n',s.t);
        ds(c).(cHandle.type) = s;
    else
        fprintf('... not available');
    end
end

% test
testID = importdata('data\testID.mat'); % identity label of all testing bboxes
testCAM = importdata('data\testCAM.mat'); % camera label of all testing bboxes
queryID = importdata('data\queryID.mat'); % identity label of all query bboxes
queryCAM = importdata('data\queryCAM.mat'); % camera label of all query bboxes
names = fieldnames(ds(c));
for nameCounter=1:length(names)
    fprintf('    evaluating %s ',upper(names{nameCounter}));
    [ds(c).(names{nameCounter}).mAP, ds(c).(names{nameCounter}).cmc] = ...
        calcmAP(ds(c).(names{nameCounter}).M, X_test, X_query, testID, testCAM, queryID, queryCAM);
    fprintf('... done \n');
end

end


text:
function [obj, p, g] = runICPInit(obj, p, g)

% Start ------------------------------------------------------------------------

% Start
msg('O', 'SetLogLevel', p.LogLevel);

% Log to file
diary off; fclose all;
p2report = fullfile(obj.OutputFolder, 'ICPLog.txt');
if exist(p2report)==2 % delete old report file if any
    try delete(p2report); end
end
if ~exist(obj.OutputFolder, 'dir')
    mkdir(obj.OutputFolder);
end
diary(p2report);

g.proc = {'GLOBALICP' 'RUNICP'};
msg('S', g.proc, 'LogLevel', 'basic');

% Initialize variables ---------------------------------------------------------

g.nPC = numel(obj.PC);
g.nItICP = 1;
obj.D.stats{g.nItICP}.normdx = Inf;

% Initialize homogeneous transformation matrices
obj.D.H = repmat({eye(4)}, 1, g.nPC);

% Estimate input parameters set to auto ----------------------------------------

if strcmpi(p.HullVoxelSize          , 'auto') || ...
   strcmpi(p.UniformSamplingDistance, 'auto') || ...
   strcmpi(p.MaxDistance            , 'auto') || ...
   strcmpi(p.MaxRoughness           , 'auto')

    msg('S', {g.proc{:} 'ESTIMATE INPUT PARAMETERS'}, 'LogLevel', 'basic');

    if strcmpi(p.HullVoxelSize, 'auto')
        firstTryHullVoxelSize = 5*p.PlaneSearchRadius;
        validVoxelSizes = [0.01 0.02 0.04 0.05 0.08 0.1 0.16 0.2 0.25 0.4 0.5 0.8 1 1.25 2 2.5 4 5 6.25 10 12.5 20 25 50 100]'; % mod(100, voxelSize) has to be zero, see pointCloud.getVoxelSize
        idx = knnsearch(validVoxelSizes, firstTryHullVoxelSize);
        p.HullVoxelSize = validVoxelSizes(idx);
        msg('V', p.HullVoxelSize, 'HullVoxelSize', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.UniformSamplingDistance, 'auto')
        p.UniformSamplingDistance = 2*p.PlaneSearchRadius;
        msg('V', p.UniformSamplingDistance, 'UniformSamplingDistance', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.MaxDistance, 'auto')
        p.MaxDistance = 3*p.PlaneSearchRadius;
        msg('V', p.MaxDistance, 'MaxDistance', 'Prec', 2, 'LogLevel', 'basic');
    end

    if strcmpi(p.MaxRoughness, 'auto')
        p.MaxRoughness = 0.2*p.PlaneSearchRadius;
        msg('V', p.MaxRoughness, 'MaxRoughness', 'Prec', 2, 'LogLevel', 'basic');
    end

    msg('E', {g.proc{:} 'ESTIMATE INPUT PARAMETERS'}, 'LogLevel', 'basic');

end

% Assign point cloud dependent parameters --------------------------------------

% Uniform sampling for each point cloud
if numel(p.UniformSamplingDistance) == 1
    p.UniformSamplingDistance = ones(1,g.nPC)*p.UniformSamplingDistance;
end

% Random subsampling for each point cloud
if numel(p.RandomSubsampling) == 1
    p.RandomSubsampling = repmat(p.RandomSubsampling, 1, g.nPC);
end

% Normal subsampling for each point cloud
if numel(p.NormalSubsampling) == 1
    p.NormalSubsampling = repmat(p.NormalSubsampling, 1, g.nPC);
end

% Max leverage subsampling for each point cloud
if numel(p.MaxLeverageSubsampling) == 1
    p.MaxLeverageSubsampling = repmat(p.MaxLeverageSubsampling, 1, g.nPC);
end

% Report input point clouds ----------------------------------------------------

msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');

msg('T', 'INPUT POINT CLOUDS:', 'LogLevel', 'basic');

msg('T', sprintf('%12s %12s %8s %8s %8s %8s    file', 'point_cloud', 'fix/loose', 'USD', 'RS', 'NS', 'MLS'), 'LogLevel', 'basic');

for i = 1:g.nPC
    if ismember(i, p.IdxFixedPointClouds), fixOrLoose = 'fix  '; else fixOrLoose = 'loose'; end
    if p.RandomSubsampling(i)     , RS  = 'true'; else RS  = 'false'; end
    if p.NormalSubsampling(i)     , NS  = 'true'; else NS  = 'false'; end
    if p.MaxLeverageSubsampling(i), MLS = 'true'; else MLS = 'false'; end

    [~, file] = fileparts(obj.PC{i});

    msg('T', sprintf('%12s %12s %8.2f %8s %8s %8s    %s', ...
        ['[' num2str(i) ']'], ...
        fixOrLoose, ...
        p.UniformSamplingDistance(i), ...
        RS, ...
        NS, ...
        MLS, ...
        file), ...
        'LogLevel', 'basic');

end

msg('T', 'where USD = UniformSamplingDistance', 'LogLevel', 'basic');
msg('T', '      RS  = RandomSubsampling'      , 'LogLevel', 'basic');
msg('T', '      NS  = NormalSubsampling'      , 'LogLevel', 'basic');
msg('T', '      MLS = MaxLeverageSubsampling' , 'LogLevel', 'basic');

msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');

end

author: pglira


text:
classdef PostProcessColumn < handle

    properties (Access = public)
        m
    end

    properties (Access = private)
        polygon
    end

    properties (Access = private)
        designVariable
        dim
        mesh
        scale
    end

    methods (Access = public)

        function obj = PostProcessColumn(cParams)
            obj.init(cParams)
        end

        function plotColumn(obj)
            obj.createPolygon();
            obj.createMesh();
            obj.plotMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.designVariable = cParams.designVariable;
            obj.dim            = cParams.dim;
            obj.mesh           = cParams.mesh;
            obj.scale          = cParams.scale;
        end

        function createPolygon(obj)
            z = obj.designVariable.getColumnArea;
            % d     = obj.dim;
            scl   = obj.scale;
            coord = obj.mesh.coord;
            nnod     = obj.mesh.nelem+1;
            dimFig   = 2;
            vertElem = 4;
            vertex = zeros(vertElem*obj.mesh.nelem+1,dimFig);
            for iNod = 1:nnod-1
                vertex(2*iNod-1,:)   = [scl*z(iNod)/2 coord(iNod)];
                vertex(2*iNod,:) = [scl*z(iNod)/2 coord(iNod+1)];
            end
            vertex = obj.flip(vertex,vertElem,dimFig);
            vertex(end,:) = vertex(1,:);
            obj.polygon = polyshape(vertex);
        end

        function createMesh(obj)
            pgon = obj.polygon;
            tr = triangulation(pgon);
            model = createpde;
            tnodes = tr.Points';
            telements = tr.ConnectivityList';
            geometryFromMesh(model,tnodes,telements);
            generateMesh(model,'GeometricOrder','linear','Hmax',0.1)
            coord = model.Mesh.Nodes';
            connec = model.Mesh.Elements';
            s.coord = coord;
            s.connec = connec;
            obj.m = Mesh(s);
        end

        function vertex = flip(obj,vertex,vertElem,dimFig)
            nElem = obj.mesh.nelem;
            nodes = dimFig*nElem+1:vertElem*nElem;
            vertex(nodes,1) = - fliplr(vertex(1:dimFig*nElem,1)')';
            vertex(dimFig*nElem+1:vertElem*nElem,2) = fliplr(vertex(1:dimFig*nElem,2)')';
        end

        function plotMesh(obj)
            figure
            clf
            obj.m.plot();
            grid on
            grid minor
            title('Optimal clamped-clamped column','Interpreter', 'latex','FontSize',20, 'fontweight','b');
            xlabel('A(x)','Interpreter', 'latex','fontsize',14,'fontweight','b');
            ylabel('x','Interpreter', 'latex','fontsize',14,'fontweight','b');
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pass = test_roots_slow( pref )
% Check that the marching squares and Bezoutian agree with each other.
%%

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3 * pref.cheb2Prefs.chebfun2eps;
j = 1;

%% (slow one)
f = chebfun2(@(x,y)exp(x-2*x.^2-y.^2).*sin(10*(x+y+x.*y.^2)));
g = chebfun2(@(x,y)exp(-x+2*y.^2+x.*y.^2).*sin(10*(x-y-2*x.*y.^2)));
r1 = roots([f;g],'ms');
r2 = roots([f;g],'resultant');
pass(j) = ( norm(sort(r1(:,1))-sort(r2(:,1))) < tol ); j = j + 1;
pass(j) = ( norm(sort(r1(:,2))-sort(r2(:,2))) < tol ); j = j + 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pass = test_inherited( )
% Test inherited commands, abs, cos, cosh, conj, etc.

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% abs
f = spherefun( @(x,y,z) -1+0*x);
h = spherefun( @(x,y,z) 1+0*x);
g = abs( f );
pass(1) = norm( g - h ) < tol;

% cos
f = spherefun( @(x,y,z) cos(x));
h = spherefun( @(x,y,z) cos(cos(x)));
g = cos( f );
pass(2) = norm( g - h ) < tol;

% cosh
f = spherefun( @(x,y,z) cos(x));
h = spherefun( @(x,y,z) cosh(cos(x)));
g = cosh( f );
pass(3) = norm( g - h ) < tol;

% conj
f = spherefun( @(x,y,z) x);
g = conj( f );
pass(4) = norm( f - g ) < tol;

% ctranspose
f = spherefun( @(x,y,z) y);
g = ctranspose( f );
pass(5) = norm( f - g ) < tol;

% diag
f = spherefun( @(lam,theta) sin(theta).*cos(lam));
h = chebfun( @(x) sin(x).*cos(x), [0 pi]);
g = diag( f );
pass(6) = norm( g - h ) < tol;

% exp
f = spherefun( @(x,y,z) cos(x)+sin(y));
h = spherefun( @(x,y,z) exp(cos(x)+sin(y)));
g = exp( f );
pass(7) = norm( g - h ) < tol;

% flipdim, fliplr, flipud
f = spherefun( @(lam, theta) sin(theta).*cos(lam));
g1 = spherefun( @(lam, theta) sin(theta).*cos(-lam));
g2 = spherefun( @(lam, theta) sin(-theta).*cos(lam));
h1 = fliplr( f );
h2 = flipud( f );
h3 = flipdim( f, 1);
h4 = flipdim( f, 2);
pass(8) = ( norm( g1 - h1 ) < tol );
pass(9) = ( norm( g2 - h2 ) < tol );
pass(10) = ( norm( g2 - h3 ) + norm( g1 - h4 ) < tol );

% imag
f = spherefun( @(x,y,z) x );
g = imag( f );
pass(11) = norm( g  ) < tol;

% isequal
f = spherefun( @(x,y,z) cos(x));
g = spherefun( @(x,y,z) cos(-x) );
pass(12) = isequal( f, g );

f = spherefun( @(x,y,z) cos(x));
g = spherefun( @(x,y,z) cos(x) + 1 );
pass(13) = ~isequal( f, g );

% isreal
f = spherefun( @(x,y,z) x);
pass(14) = isreal(f);

f = spherefun( @(x,y,z) cos(x.*y.*z));
pass(15) = isreal(f);

% iszero
f = spherefun( @(x,y,z) 0*cos(x.*y.*z));
pass(16) = iszero( f );
f = spherefun( zeros(10) );
pass(17) = iszero( f );

% length
f = spherefun( @(x,y,z) cos(z));
[m, n] = length( f );
pass(18) = (m == 1);

% log
f = spherefun( @(x,y,z) exp(z));
h = spherefun( @(x,y,z) z);
g = log( f );
pass(19) = norm( g - h ) < tol;

% rank

% sin
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) sin(cos(x.*y)));
g = sin( f );
pass(20) = norm( g - h ) < tol;

% sinh
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) sinh(cos(x.*y)));
g = sinh( f );
pass(21) = norm( g - h ) < tol;

% size
pass(22) = (size( f,1 ) == inf && size( f,2 ) == inf);

% sqrt
f = spherefun( @(x,y,z) cos(x.*y).^2);
h = spherefun( @(x,y,z) cos(x.*y));
g = sqrt( f );
pass(23) = norm( g - h ) < tol;

% tan
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tan(cos(x.*y)));
g = tan( f );
pass(24) = norm( g - h ) < 10*tol;

% tand
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tand(cos(x.*y)));
g = tand( f );
pass(25) = norm( g - h ) < tol;

% tanh
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) tanh(cos(x.*y)));
g = tanh( f );
pass(26) = norm( g - h ) < tol;

% uminus
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) -(cos(x.*y)));
g = uminus( f );
pass(27) = norm( g - h ) < tol;

% uplus
f = spherefun( @(x,y,z) cos(x.*y));
h = spherefun( @(x,y,z) +(cos(x.*y)));
g = uplus( f );
pass(28) = norm( g - h ) < tol;

end


text:
function meta=recover_meta(in, dirac)

%meta.order
%meta.scale
%meta.orientation
%meta.dirac_norm
%meta.ave

J=in{1}.meta.j;
L=max(in{2}.meta.theta);

meta.order(1)=1;
meta.scale(1)=-1;
meta.orientation(1)=0;
meta.dirac_norm(1)=norm(dirac{1}.signal{1}(:));
meta.ave(1)=mean(in{1}.signal{1}(:));
r=2;
for m=2:size(in,2)
for l=1:size(in{m}.signal,2)
meta.order(r)=m;
meta.scale(r)=simplestack(in{m}.meta.j(1:end-1,l),J,m-1);
meta.orientation(r)=simplestack(in{m}.meta.theta(:,l)-1,L,m-1);
meta.dirac_norm(r)=norm(dirac{m}.signal{l}(:));
meta.ave(r)=mean(in{m}.signal{l}(:));
r=r+1;
end
end

end

function out=simplestack(code, base, len)
out=0;
for l=1:len
out=base*out+code(l);
end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [E,J]=SynthMeasIsoGPD(d, protocol)
% Computes signals and their derivatives for isotropic free diffusion with
% diffusivity d for the protocol.
%
% [E,J]=SynthMeasIsoGPD(d, protocol)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameter d.
%
% d is the diffusivity.
%
% protocol is the object containing the acquisition protocol.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

if(strcmp(protocol.pulseseq, 'PGSE') || strcmp(protocol.pulseseq, 'STEAM'))

    GAMMA = 2.675987E8;
    modQ = GAMMA*protocol.smalldel'.*protocol.G';
    modQ_Sq = modQ.^2;
    difftime = (protocol.delta'-protocol.smalldel'/3);

    E = exp(-difftime.*modQ_Sq*d);
    if(nargout>1)
        J = -difftime.*modQ_Sq.*E;
    end
elseif(strcmp(protocol.pulseseq, 'FullSTEAM'))

    GAMMA = 2.675987E8;

    FullG = protocol.grad_dirs.*repmat(protocol.G', [1 3]);
    GdGd = sum(FullG.*FullG,2)';
    GcGc = sum(protocol.cG.*protocol.cG,2)';
    GrGr = sum(protocol.rG.*protocol.rG,2)';
    GdGc = sum(FullG.*protocol.cG,2)';
    GdGr = sum(FullG.*protocol.rG,2)';
    GcGr = sum(protocol.cG.*protocol.rG,2)';

    tdd = protocol.gap1 + protocol.gap2 + protocol.TM + 2*protocol.sdelc + 2*protocol.smalldel/3 + 2*protocol.sdelr;
    tcc = protocol.TM + 2*protocol.sdelc/3 + 2*protocol.sdelr;
    trr = protocol.TM + 2*protocol.sdelr/3;
    tdc = protocol.TM + protocol.sdelc + 2*protocol.sdelr;
    tdr = protocol.TM + protocol.sdelr;
    tcr = protocol.TM + protocol.sdelr;

    gqt = GAMMA^2*(GdGd.*protocol.smalldel.^2.*tdd + ...
    GcGc.*protocol.sdelc.^2.*tcc + ...
    GrGr.*protocol.sdelr.^2.*trr + ...
    2*GdGc.*protocol.smalldel.*protocol.sdelc.*tdc + ...
    2*GdGr.*protocol.smalldel.*protocol.sdelr.*tdr + ...
    2*GcGr.*protocol.sdelc.*protocol.sdelr.*tcr);

    E=exp(-gqt'*d);

    if(nargout>1)
        J = -gqt.*E;
    end
elseif(strcmp(protocol.pulseseq, 'DSE'))
    bValue = GetB_ValuesDSE(protocol.G', protocol.delta1', protocol.delta2', protocol.delta3', protocol.t1', protocol.t2', protocol.t3');
    E = exp(-bValue*d);
    if(nargout>1)
        J = -bValue.*E;
    end
elseif(strcmp(protocol.pulseseq, 'OGSE'))
    bValue = GetB_Values(protocol)';
    E = exp(-bValue*d);
    if(nargout>1)
        J = -bValue.*E;
    end
else
    error('Not implemented for pulse sequence: %s', protocol.pulseseq);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function varargout=randiscr(varargin)
%
% For calling details please see v_randiscr.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_randiscr.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: randiscr.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_randiscr(varargin{:});
else
    v_randiscr(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function mri = ft_defacevolume(cfg, mri)

% FT_DEFACEVOLUME allows you to de-identify an anatomical MRI by erasing specific
% regions, such as the face and ears. The graphical user interface allows you to
% position a box over the anatomical data inside which all anatomical voxel values will
% be replaced by zero. You might have to call this function multiple times when both
% face and ears need to be removed. Following defacing, you should check the result
% with FT_SOURCEPLOT.
%
% Use as
%   mri = ft_defacevolume(cfg, mri)
%
% The configuration can contain the following options
%   cfg.method     = 'box', 'spm' (default = 'box')
%
% If you specify the box method, the following options apply
%   cfg.translate  = initial position of the center of the box (default = [0 0 0])
%   cfg.scale      = initial size of the box along each dimension (default is automatic)
%   cfg.rotate     = initial rotation of the box (default = [0 0 0])
%   cfg.selection  = which voxels to keep, can be 'inside' or 'outside' (default = 'outside')
%   cfg.smooth     = 'no' or the FWHM of the gaussian kernel in voxels (default = 'no')
%   cfg.keepbrain  = 'no' or 'yes', segment and retain the brain (default = 'no')
%   cfg.feedback   = 'no' or 'yes', whether to provide graphical feedback (default = 'no')
%
% If you specify no smoothing, the selected area will be zero-masked. If you
% specify a certain amount of smoothing (in voxels FWHM), the selected area will
% be replaced by a smoothed version of the data.
%
% The spm method does not have any options, it uses SPM_DEFACE from the
% SPM12 toolbox.
%
% See also FT_ANONYMIZEDATA, FT_DEFACEMESH, FT_ANALYSISPIPELINE, FT_SOURCEPLOT

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    mri
ft_preamble provenance mri
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% set the defaults
cfg.method    = ft_getopt(cfg, 'method', 'box');
cfg.rotate    = ft_getopt(cfg, 'rotate', [0 0 0]);
cfg.scale     = ft_getopt(cfg, 'scale'); % the automatic default is determined further down
cfg.translate = ft_getopt(cfg, 'translate', [0 0 0]);
cfg.selection = ft_getopt(cfg, 'selection', 'outside');
cfg.smooth    = ft_getopt(cfg, 'smooth', 'no');
cfg.keepbrain = ft_getopt(cfg, 'keepbrain', 'no');
cfg.feedback  = ft_getopt(cfg, 'feedback', 'no');

ismri    = ft_datatype(mri, 'volume') && isfield(mri, 'anatomy');
ismesh   = isfield(mri, 'pos'); % triangles are optional

if ismri
  % check if the input data is valid for this function
  mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes');
end

switch cfg.method
  case 'spm'
    % this requires SPM12 on the path
    ft_hastoolbox('spm12', 1);

    % defacing relies on coregistration, which relies on the MRI being reasonably aligned for SPM
    mri = ft_checkdata(mri, 'hascoordsys', 'yes');

    % remember the original transformation matrix and coordinate system
    original = [];
    original.transform = mri.transform;
    original.coordsys  = mri.coordsys;
    mri = ft_convert_coordsys(mri, 'acpc');

    filename1 = {[tempname '.nii']};
    ft_write_mri(filename1{1}, mri, 'dataformat', 'nifti');

    % % apply a least squares pre-alignment step in order to make spm_deface more robust
    % % this could be done conditional on the modality/contrast, which is part of the BIDS filename
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152PD.nii'));
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152T1.nii'));
    % template = spm_vol(fullfile(spm('Dir'),'canonical','avg152T2.nii'));
    % filevol = spm_vol(filename1{1});
    % M = spm_affreg(template, filevol);
    % spm_get_space(filename1{1}, M * filevol.mat);

    filename2 = spm_deface(filename1);
    mri = ft_read_mri(filename2{1});

    % put the original transformation matrix and coordinate system back
    mri.transform = original.transform;
    mri.coordsys = original.coordsys;

    % clean up the temporary files
    delete(filename1{1});
    delete(filename2{1});

  case 'box'
    % determine the size of the "unit" sphere in the origin and the length of the axes
    switch mri.unit
      case 'mm'
        axmax = 150;
        rbol  = 5;
      case 'cm'
        axmax = 15;
        rbol  = 0.5;
      case 'm'
        axmax = 0.15;
        rbol  = 0.005;
      otherwise
        ft_error('unknown units "%s"', unit);
    end

    figHandle = figure;
    set(figHandle, 'CloseRequestFcn', @cb_close);

    % clear persistent variables to ensure fresh figure
    clear ft_plot_slice

    if ismri
      % the volumetric data needs to be interpolated onto three orthogonal planes
      % determine a resolution that is close to, or identical to the original resolution
      [corner_vox, corner_head] = cornerpoints(mri.dim, mri.transform);
      diagonal_head = norm(range(corner_head));
      diagonal_vox  = norm(range(corner_vox));
      resolution    = diagonal_head/diagonal_vox; % this is in units of "mri.unit"

      % create a contrast enhanced version of the anatomy
      mri.anatomy = double(mri.anatomy);
      dum = unique(mri.anatomy(:));
      clim(1) = dum(round(0.05*numel(dum)));
      clim(2) = dum(round(0.95*numel(dum)));
      anatomy = (mri.anatomy-clim(1))/(clim(2)-clim(1));

      ft_plot_ortho(anatomy, 'transform', mri.transform, 'unit', mri.unit, 'resolution', resolution, 'style', 'intersect');
    elseif ismesh
      if isfield(mri, 'hex')
        ft_plot_mesh(mri, 'surfaceonly', 'yes');
      else
        ft_plot_mesh(mri);
      end
    end

    axis vis3d
    view([110 36]);

    % shift the axes to the left
    ax = get(gca, 'position');
    ax(1) = 0;
    set(gca, 'position', ax);

    % get the xyz-axes
    xdat  = [-axmax 0 0; axmax 0 0];
    ydat  = [0 -axmax 0; 0 axmax 0];
    zdat  = [0 0 -axmax; 0 0 axmax];

    % get the xyz-axes dotted
    xdatdot = (-axmax:(axmax/15):axmax);
    xdatdot = xdatdot(1:floor(numel(xdatdot)/2)*2);
    xdatdot = reshape(xdatdot, [2 numel(xdatdot)/2]);
    n       = size(xdatdot,2);
    ydatdot = [zeros(2,n) xdatdot zeros(2,n)];
    zdatdot = [zeros(2,2*n) xdatdot];
    xdatdot = [xdatdot zeros(2,2*n)];

    % plot axes
    hl = line(xdat, ydat, zdat);
    set(hl(1), 'linewidth', 1, 'color', 'r');
    set(hl(2), 'linewidth', 1, 'color', 'g');
    set(hl(3), 'linewidth', 1, 'color', 'b');
    hld = line(xdatdot, ydatdot, zdatdot);
    for k = 1:n
      set(hld(k    ), 'linewidth', 3, 'color', 'r');
      set(hld(k+n*1), 'linewidth', 3, 'color', 'g');
      set(hld(k+n*2), 'linewidth', 3, 'color', 'b');
    end

    if isempty(cfg.scale)
      cfg.scale = [axmax axmax axmax]/2;
    end

    guidata(figHandle, cfg);

    % add the GUI elements
    cb_creategui(gca);
    cb_redraw(gca);

    uiwait(figHandle);
    cfg = guidata(figHandle);
    delete(figHandle);
    drawnow
    fprintf('keeping all voxels from MRI that are %s the box\n', cfg.selection)

    % the order of application is scale, rotate, translate
    S = cfg.S;
    R = cfg.R;
    T = cfg.T;

    if ismri
      % it is possible to convert the box to headcoordinates, but it is more efficient the other way around
      [X, Y, Z] = ndgrid(1:mri.dim(1), 1:mri.dim(2), 1:mri.dim(3));
      voxpos = ft_warp_apply(mri.transform, [X(:) Y(:) Z(:)]);  % voxel positions in head coordinates
      voxpos = ft_warp_apply(inv(T*R*S), voxpos);               % voxel positions in box coordinates

      remove = ...
        voxpos(:,1) > -0.5 & ...
        voxpos(:,1) < +0.5 & ...
        voxpos(:,2) > -0.5 & ...
        voxpos(:,2) < +0.5 & ...
        voxpos(:,3) > -0.5 & ...
        voxpos(:,3) < +0.5;

    elseif ismesh || issource
      meshpos = ft_warp_apply(inv(T*R*S), mri.pos);               % mesh vertex positions in box coordinates

      remove = ...
        meshpos(:,1) > -0.5 & ...
        meshpos(:,1) < +0.5 & ...
        meshpos(:,2) > -0.5 & ...
        meshpos(:,2) < +0.5 & ...
        meshpos(:,3) > -0.5 & ...
        meshpos(:,3) < +0.5;
    end

    if strcmp(cfg.selection, 'inside')
      % invert the selection, i.e. keep the voxels inside the box
      remove = ~remove;
    end

    if ismri
      if istrue(cfg.keepbrain)
        tmpcfg = [];
        tmpcfg.output = {'brain'};
        seg = ft_volumesegment(tmpcfg, mri);
        fprintf('keeping voxels in brain segmentation\n');
        % keep the tissue of the brain
        remove(seg.brain) = 0;
        clear seg
      end

      if istrue(cfg.feedback)
        tmpmri = keepfields(mri, {'anatomy', 'transform', 'coordsys', 'units', 'dim'});
        tmpmri.remove = remove;
        tmpcfg = [];
        tmpcfg.funparameter = 'remove';
        ft_sourceplot(tmpcfg, tmpmri);
      end

      if isequal(cfg.smooth, 'no')
        fprintf('zero-filling %.0f%% of the volume\n', 100*mean(remove));
        mri.anatomy(remove) = 0;
      else
        tmp = mri.anatomy;
        tmp = (1 + 0.5.*randn(size(tmp))).*tmp; % add 50% noise to each voxel
        tmp = volumesmooth(tmp, cfg.smooth, 'anatomy');
        fprintf('smoothing %.0f%% of the volume\n', 100*mean(remove));
        mri.anatomy(remove) = tmp(remove);
      end

    elseif ismesh
      % determine all fields that might need to be defaced
      fn = setdiff(fieldnames(mri), ignorefields('deface'));
      dimord = cell(size(fn));
      for i=1:numel(fn)
        dimord{i} = getdimord(mri, fn{i});
      end
      % this applies to headshapes and meshes in general
      fprintf('keeping %d and removing %d vertices in the mesh\n', sum(remove==0), sum(remove==1));
      if isfield(mri, 'tri')
        [mri.pos, mri.tri] = remove_vertices(mri.pos, mri.tri, remove);
      elseif isfield(mri, 'tet')
        [mri.pos, mri.tet] = remove_vertices(mri.pos, mri.tet, remove);
      elseif isfield(mri, 'hex')
        [mri.pos, mri.hex] = remove_vertices(mri.pos, mri.hex, remove);
      else
        mri.pos = mri.pos(~remove,1:3);
      end
      for i=1:numel(fn)
        dimtok = tokenize(dimord{i}, '_');
        % do some sanity checks
        if any(strcmp(dimtok, '{pos}'))
          ft_error('not supported');
        end
        if numel(dimtok)>5
          ft_error('too many dimensions');
        end
        % remove the same positions from each matching dimension
        if numel(dimtok)>0 && strcmp(dimtok{1}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(~remove,:,:,:,:);
        end
        if numel(dimtok)>1 && strcmp(dimtok{2}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,~remove,:,:,:);
        end
        if numel(dimtok)>2 && strcmp(dimtok{3}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,~remove,:,:);
        end
        if numel(dimtok)>3 && strcmp(dimtok{4}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,:,~remove,:);
        end
        if numel(dimtok)>4 && strcmp(dimtok{5}, 'pos')
          mri.(fn{i}) = mri.(fn{i})(:,:,:,:,~remove);
        end
      end % for fn
      mri = removefields(mri, {'dim', 'transform'}); % these fields don't apply any more
    end % ismesh

    % remove the temporary fields from the configuration, keep the rest for provenance
    cfg = removefields(cfg, {'R', 'S', 'T'});
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous mri
ft_postamble provenance mri
ft_postamble history mri
ft_postamble savevar mri

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_redraw(figHandle, varargin)
persistent p
% define the position of each GUI element
figHandle = get(figHandle, 'parent');
cfg = guidata(figHandle);

rx = str2double(get(findobj(figHandle, 'tag', 'rx'), 'string'));
ry = str2double(get(findobj(figHandle, 'tag', 'ry'), 'string'));
rz = str2double(get(findobj(figHandle, 'tag', 'rz'), 'string'));
tx = str2double(get(findobj(figHandle, 'tag', 'tx'), 'string'));
ty = str2double(get(findobj(figHandle, 'tag', 'ty'), 'string'));
tz = str2double(get(findobj(figHandle, 'tag', 'tz'), 'string'));
sx = str2double(get(findobj(figHandle, 'tag', 'sx'), 'string'));
sy = str2double(get(findobj(figHandle, 'tag', 'sy'), 'string'));
sz = str2double(get(findobj(figHandle, 'tag', 'sz'), 'string'));

% remember the user specified transformation
cfg.rotate    = [rx ry rz];
cfg.translate = [tx ty tz];
cfg.scale     = [sx sy sz];

R = rotate   (cfg.rotate);
T = translate(cfg.translate);
S = scale    (cfg.scale);

% remember the transformation matrices
cfg.R = R;
cfg.T = T;
cfg.S = S;

% start with a cube of unit dimensions
x1 = -0.5;
y1 = -0.5;
z1 = -0.5;
x2 = +0.5;
y2 = +0.5;
z2 = +0.5;

plane1 = [
  x1 y1 z1
  x2 y1 z1
  x2 y2 z1
  x1 y2 z1];

plane2 = [
  x1 y1 z2
  x2 y1 z2
  x2 y2 z2
  x1 y2 z2];

plane3 = [
  x1 y1 z1
  x1 y2 z1
  x1 y2 z2
  x1 y1 z2];

plane4 = [
  x2 y1 z1
  x2 y2 z1
  x2 y2 z2
  x2 y1 z2];

plane5 = [
  x1 y1 z1
  x2 y1 z1
  x2 y1 z2
  x1 y1 z2];

plane6 = [
  x1 y2 z1
  x2 y2 z1
  x2 y2 z2
  x1 y2 z2];

plane1 = ft_warp_apply(T*R*S, plane1);
plane2 = ft_warp_apply(T*R*S, plane2);
plane3 = ft_warp_apply(T*R*S, plane3);
plane4 = ft_warp_apply(T*R*S, plane4);
plane5 = ft_warp_apply(T*R*S, plane5);
plane6 = ft_warp_apply(T*R*S, plane6);

if all(ishandle(p))
  delete(p);
end

p(1) = patch(plane1(:,1), plane1(:,2), plane1(:,3), 'y');
p(2) = patch(plane2(:,1), plane2(:,2), plane2(:,3), 'y');
p(3) = patch(plane3(:,1), plane3(:,2), plane3(:,3), 'y');
p(4) = patch(plane4(:,1), plane4(:,2), plane4(:,3), 'y');
p(5) = patch(plane5(:,1), plane5(:,2), plane5(:,3), 'y');
p(6) = patch(plane6(:,1), plane6(:,2), plane6(:,3), 'y');
set(p, 'FaceAlpha', 0.3);

guidata(figHandle, cfg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_creategui(figHandle, varargin)
% define the position of each GUI element
figHandle = get(figHandle, 'parent');
cfg = guidata(figHandle);

% constants
CONTROL_WIDTH   = 0.05;
CONTROL_HEIGHT  = 0.08;
CONTROL_HOFFSET = 0.68;
CONTROL_VOFFSET = 0.20;

% rotateui
uicontrol('tag', 'rotateui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'rotate', 'callback', [])
uicontrol('tag', 'rx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(1)), 'callback', @cb_redraw)
uicontrol('tag', 'ry', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(2)), 'callback', @cb_redraw)
uicontrol('tag', 'rz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.rotate(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'rotateui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'rx',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'ry',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);
ft_uilayout(figHandle, 'tag', 'rz',       'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET);

% scaleui
uicontrol('tag', 'scaleui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'scale', 'callback', [])
uicontrol('tag', 'sx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(1)), 'callback', @cb_redraw)
uicontrol('tag', 'sy', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(2)), 'callback', @cb_redraw)
uicontrol('tag', 'sz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.scale(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'scaleui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sx',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sy',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'sz',      'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-CONTROL_HEIGHT);

% translateui
uicontrol('tag', 'translateui', 'parent', figHandle, 'units', 'normalized', 'style', 'text', 'string', 'translate', 'callback', [])
uicontrol('tag', 'tx', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(1)), 'callback', @cb_redraw)
uicontrol('tag', 'ty', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(2)), 'callback', @cb_redraw)
uicontrol('tag', 'tz', 'parent', figHandle, 'units', 'normalized', 'style', 'edit', 'string', num2str(cfg.translate(3)), 'callback', @cb_redraw)
ft_uilayout(figHandle, 'tag', 'translateui', 'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET,                 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'tx',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'ty',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);
ft_uilayout(figHandle, 'tag', 'tz',          'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH,   'height', CONTROL_HEIGHT/2, 'hpos', CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos', CONTROL_VOFFSET-2*CONTROL_HEIGHT);

% somehow the toolbar gets lost in 2012b
set(figHandle, 'toolbar', 'figure');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_close(figHandle, varargin)
% the figure will be closed in the main function after collecting the guidata
uiresume;

author: spm


text:
function [metabolite_structure] =createNewMetaboliteStructure(input,source,metabolite_structure_rBioNet,metab,rxnDB)
% This function creates a metabolite structure using the provided input. If
% no VMHId is provided in the header of the input, then VMHId are
% generated.
%
% INPUT
% Input                          Cell array containing the metabolites
%                                   The information provided must be as follows:
%                                   metList={
%                                   'VMH ID' 'metabolite_name' 'HMDB' 'inchistring' 'neutral_formula' 'charged_formula' 'charge'
%                                   'cot' 'Cotinine' 'HMDB0001046'  '' '' '' ''
%                                   'coto' 'Cotinine n-oxide' 'HMDB0001411'  '' '' '' ''
% source                            source of the information contained in metArray
%                                   (e.g., 'Manually assembled by IT')
% metabolite_structure_rBioNet      (optional and only necessary if the
%                                   input has no VMHId)
%
% OUTPUT
% metabolite_structure              metabolite structure
%
% Ines Thiele, 09/2021

%load rbionet data
rBioNetPath =  fileparts(which('tutorial_MetaboAnnotator'));
if exist([rBioNetPath filesep 'cache' filesep 'metab.mat'],'file')
    load([rBioNetPath filesep 'cache' filesep 'metab.mat']);
elseif exist([rBioNetPath filesep 'data' filesep 'metab.mat'],'file')
    load([rBioNetPath filesep 'data' filesep 'metab.mat']);
else
    %TODO
end

if exist([rBioNetPath filesep 'cache' filesep 'rxn.mat'],'file')
    load([rBioNetPath filesep 'cache' filesep 'rxn.mat']);
elseif  exist([rBioNetPath filesep 'data' filesep 'rxn.mat'],'file')
    load([rBioNetPath filesep 'data' filesep 'rxn.mat']);
else
    %TODO
end

RAW = input;

% load the designated field names associated with each metabolites
metaboliteStructureFieldNames;
% find VMH ID in the input file
vmh_col = find(contains(lower(RAW(1,:)),'vmh'));
% clean up some known issues:
name_col = find(contains(lower(RAW(1,:)),'name'));

% test if metabolite name column exists
if isempty(name_col)
    RAW{1} = 'metNames';
    RAW = cat(1, RAW, metList);
else
    RAW(1,name_col) = regexprep(RAW(1,name_col),RAW(1,name_col),'metNames');
end

if isempty(vmh_col) || vmh_col ==0
    if ~exist('metabolite_structure_rBioNet','var')
        load met_strc_rBioNet;
    end
end

% populate metabolite structure
for i = 2 : size(RAW,1)
    name = RAW(i,find(contains(lower(RAW(1,:)),'name')));
    if ~isempty(vmh_col) && vmh_col ~=0 % VMH ID exists
        Ori = RAW{i,vmh_col};

        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},'-','_minus_');
        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},'(','_parentO_');
        RAW{i,vmh_col} = regexprep(RAW{i,vmh_col},')','_parentC_');
        VMHId = RAW{i,vmh_col};
    else
        [VMHId] = generateVMHMetAbbr(RAW{i,name_col},metabolite_structure_rBioNet,metab,rxnDB);
        Ori = VMHId;
        vmh_col = 0;
    end
    % adjust the abbreviations to the requirements of Matlab field
    % names
    % the original designated VMH ids listed in the field VMHId

    %make a temporary VMHId name in case the existing one is an invalid field name
    VMHId_suggestion = VMHId;
    if  isvarname(strcat('VMHId_', VMHId))
        fieldname = strcat('VMH_',VMHId);
    else

        VMHId_suggestion = VMHId;
        %VMHId = datestr(now,'yyyymmddTHHMMSS');

        fieldname = strcat('VMH_',VMHId);
    end

    metabolite_structure.(fieldname) = struct();
    metabolite_structure.(fieldname).VMHId_suggestion = VMHId_suggestion;

    for j = 1:size(RAW,2)
        if j~=vmh_col
            % check whether a columnn header is part of the field2Add
            % definition
            colName1 =  regexprep(RAW(1,j),' ','_');
            colName2 =  regexprep(RAW(1,j),' ','');
            if ~isempty(find(contains(lower(field2Add),lower(colName1)))) || ~isempty(find(contains(lower(field2Add),lower(colName2))))
                clear F
                % assign correct field name
                if ~isempty(find(contains(lower(field2Add),lower(colName1))))
                    F=  field2Add(find(contains(lower(field2Add),lower(colName1))));
                elseif ~isempty(find(contains(lower(field2Add),lower(colName2))))
                    F=  field2Add(find(contains(lower(field2Add),lower(colName2))));
                end
                % make sure that the novel entry is not NaN or empty
                if isempty(find(isnan(RAW{i,j}))) && ~isempty(RAW{i,j}) && ~isempty(F)
                    metabolite_structure.(strcat('VMH_',VMHId)).(F{1}) = RAW{i,j};
                    metabolite_structure.(strcat('VMH_',VMHId)).([F{1},'_source']) = [source ': ' datestr(now)];
                else
                    metabolite_structure.(strcat('VMH_',VMHId)).(F{1}) = NaN;
                    metabolite_structure.(strcat('VMH_',VMHId)).([F{1},'_source']) = NaN;
                end
                metabolite_structure.(strcat('VMH_',VMHId)).VMHId = Ori;
            end
        end
    end
end
% add the remaining fields to the metabolite structure
metabolite_structure= addField2MetStructure(metabolite_structure);
% clean up known, potential issues with the input data
[metabolite_structure] = cleanUpMetabolite_structure(metabolite_structure);
% check for known, potentially remaining issues with the metabolite IDs
% any errors will be removed.
removeErrors = 1;
[metabolite_structure,errorFlag] = sanityCheckMetIds(metabolite_structure,removeErrors);

author: opencobra

repo: cobratoolbox


text:
%The field of definition of the quadratic function
%have been known,
%solve the range of values.
%rov(a,b,c,x1,x2),(x1<x2)
function rov(a,b,c,x1,x2)
a=sym(a);b=sym(b);c=sym(c);x1=sym(x1);x2=sym(x2);
dd=simple((-b)/(2*a));jj=simple((4*a*c-b^2)/(4*a));
fx1=simple(a*(x1^2)+b*x1+c);fx2=simple(a*(x2^2)+b*x2+c);
if double(a)>0
   if double(x2)<=double(dd)
      fmax=fx1
      fmin=fx2
   elseif double(x1)>=double(dd)
      fmax=fx2
      fmin=fx1
   elseif double(x1)<double(dd)&double(x2)>double(dd)
      fmax=big(fx1,fx2)
      fmin=jj
   end
else
   if double(x2)<=double(dd)
      fmax=fx2
      fmin=fx1
   elseif double(x1)>=double(dd)
      fmax=fx1
      fmin=fx2
   elseif double(x1)<double(dd)&double(x2)>double(dd)
      fmax=jj
      fmin=small(fx1,fx2)
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function r8vec_midspace_test ( )

%*****************************************************************************80
%
%% R8VEC_MIDSPACE_TEST tests R8VEC_MIDSPACE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_MIDSPACE_TEST\n' );
  fprintf ( 1, '  For a R8VEC:\n' );
  fprintf ( 1, '  R8VEC_MIDSPACE: evenly spaced midpoints between A and B\n' );

  a = 10.0;
  b = 20.0;

  x = r8vec_midspace ( n, a, b );
  r8vec_print ( n, x, '  r8vec_midspace ( 5, 10, 20 )' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function output = DesignBlitEqualizer(raw_BLIT_response)
% DesignBlitEqualizer -- Equalize raw BLIT (Band Limitted Impulse Train)
% to have flat gain using FIR equalizer.
% This equalizer assumes Nuttall-based BLIT with zero at 0.55 times
% sampling frequency.
% This is a supporting function for  fractional pitch allocation of
% excitation pulse
%
% output = DesignBlitEqualizer(raw_BLIT_response)
%
% Input argument
%   raw_BLIT_response: impulse response of antialiasing BLIT filter
%
% Output
%   output : FIR equalizer to flatten frequency response

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

fftl = 8192; % Finally the central 121 points are used.
% The truncation length 121 is decided for the response to decrease
% relatively small enough.
amp = abs(fft(raw_BLIT_response / sum(raw_BLIT_response), fftl));
eqAmp = 1.0 ./ (amp + 10.0 ^ (-45 / 20));
output = real(fftshift(ifft(eqAmp)));
output = output(fftl / 2 + 1 + (-60:60));
end

author: google


text:
function blas1_d_test12 ( )

%*****************************************************************************80
%
%% TEST12 tests IDAMAX, DAXPY, DDOT, DSCAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST12\n' );
  fprintf ( 1, '  Use IDAMAX, DAXPY, DDOT and DSCAL\n' );
  fprintf ( 1, '  in a Gauss elimination routine.\n' );
  fprintf ( 1, '\n' );
%
%  Set the matrix.
%
  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 2.0;
      elseif ( i == j + 1 )
        a(i,j) = - 1.0;
      elseif ( i == j - 1 )
        a(i,j) = - 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end
%
%  Set the right hand side.
%
  b(1:n-1) = 0.0;
  b(n) = n + 1.0;

  info = 0;

  for k = 1 : n - 1

    l = idamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;

    if ( a(l,k) == 0.0 )

      info = k;

    else

      if ( l ~= k )
        t = a(l,k);
        a(l,k) = a(k,k);
        a(k,k) = t;
      end

      t = -1.0 / a(k,k);
      a(k+1:n,k) = dscal ( n-k, t, a(k+1:n,k), 1 );

      for j = k+1 : n
        t = a(l,j);
        if ( l ~= k )
          a(l,j) = a(k,j);
          a(k,j) = t;
        end
        a(k+1:n,j) = daxpy ( n-k, t, a(k+1:n,k), 1, a(k+1:n,j), 1 );
      end

    end

  end

  ipvt(n) = n;
  if ( a(n,n) == 0.0 )
    info = n;
  end

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The matrix is singular.\n' );
    return
  end

  for k = 1 : n-1
    l = ipvt(k);
    t = b(l);
    if ( l ~= k )
      b(l) = b(k);
      b(k) = t;
    end
    b(k+1:n) = daxpy ( n-k, t, a(k+1:n,k)', 1, b(k+1:n), 1 );
  end

  for k = n : -1 : 1
    b(k) = b(k) / a(k,k);
    t = - b(k);
    b(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, b(1:k-1), 1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution (should be (1,2,3,...,n)):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12f\n', b(i) );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Trdata,Vldata] = getTrVlData(TrsRGB_dir, TrXYZ_dir, ...
    VlsRGB_dir, VlXYZ_dir, aug, Tr_ImageNum, Vl_ImageNum, imgSize)

%get training imgs
Tr_in_images = imageDatastore(TrsRGB_dir, 'ReadFcn', @inRead);

%get training imgs GT
Tr_gt_images = imageDatastore(TrsRGB_dir);

Tr_gt_images.ReadFcn = @(filename)gtRead(filename, ...
    TrsRGB_dir, TrXYZ_dir);

%get random images instead of using the entire dataset
if Tr_ImageNum ~=0  && Tr_ImageNum < length(Tr_in_images.Files)
    inds = randperm(Tr_ImageNum);
    Tr_in_images.Files = Tr_in_images.Files(inds(1:Tr_ImageNum));
    Tr_gt_images.Files = Tr_gt_images.Files(inds(1:Tr_ImageNum));
end

%get validation imgs
Vl_in_images = imageDatastore(VlsRGB_dir, 'ReadFcn', @inRead);

%get validation imgs GT
Vl_gt_images = imageDatastore(VlsRGB_dir);

Vl_gt_images.ReadFcn = @(filename)gtRead(filename, ...
    VlsRGB_dir, VlXYZ_dir);

%get random images instead of using the entire dataset
if Vl_ImageNum ~=0  && Vl_ImageNum < length(Vl_in_images.Files)
    inds = randperm(Vl_ImageNum);
    Vl_in_images.Files = Vl_in_images.Files(inds(1:Vl_ImageNum));
    Vl_gt_images.Files = Vl_gt_images.Files(inds(1:Vl_ImageNum));
end

if aug == 1
    augmenter = imageDataAugmenter( ...
        'RandScale',[1 1.2], ... %scaling
        'RandXReflection',true, ... %X reflection
        'RandYReflection',true); %Y reflection

    Trdata = randomPatchExtractionDatastore(Tr_in_images,Tr_gt_images, ...
        imgSize,'DataAugmentation',augmenter,'PatchesPerImage',...
        1);
elseif aug == 0
    Trdata = randomPatchExtractionDatastore(Tr_in_images,Tr_gt_images, ...
        imgSize,'PatchesPerImage', 1);
end

Vldata = randomPatchExtractionDatastore(Vl_in_images,Vl_gt_images, ...
        imgSize,'PatchesPerImage', 2);

author: mahmoudnafifi


text:
function x = VBA_spm_invBcdf(F,v,w,tol)
% Inverse Cumulative Distribution Function (CDF) of Beta distribution
% FORMAT x = spm_invBcdf(F,v,w,tol)
%
% F   - CDF (lower tail p-value)
% v   - Shape parameter (v>0)
% w   - Shape parameter (w>0)
% x   - Beta ordinates at which CDF F(x)=F
% tol - tolerance [default 10^-6]
%__________________________________________________________________________
%
% spm_invBcdf implements the inverse Cumulative Distribution Function
% for Beta distributions.
%
% Returns the Beta-variate x such that Pr{X<x} = F for X a Beta
% random variable with shape parameters v and w.
%
% Definition:
%--------------------------------------------------------------------------
% The Beta distribution has two shape parameters, v and w, and is
% defined for v>0 & w>0 and for x in [0,1] (See Evans et al., Ch5).
% The Cumulative Distribution Function (CDF) F(x) is the probability
% that a realisation of a Beta random variable X has value less than
% x. F(x)=Pr{X<x}: This function is usually known as the incomplete Beta
% function. See Abramowitz & Stegun, 26.5; Press et al., Sec6.4 for
% definitions of the incomplete beta function.
%
% Variate relationships:
%--------------------------------------------------------------------------
% Many: See Evans et al., Ch5
%
% Algorithm:
%--------------------------------------------------------------------------
% Interval bisection to find the inverse CDF to within tol.
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invBcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Parameters
%--------------------------------------------------------------------------
Dtol  = 10^-8;
maxIt = 10000;

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<4, tol = Dtol; end
if nargin<3, error('Insufficient arguments'), end

ad = [ndims(F);ndims(v);ndims(w)];
rd = max(ad);
as = [[size(F),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation - Undefined and special cases
%--------------------------------------------------------------------------
%-Initialise result to zeros
x = zeros(rs);

%-Only defined for F in [0,1] & strictly positive v & w.
% Return NaN if undefined.
md = ( F>=0  &  F<=1  &  v>0  &  w>0 );
if any(~md(:))
    x(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special cases: x=0 when F=0, x=1 when F=1
x(md & F==1) = 1;

%-Compute where defined & not special case
%--------------------------------------------------------------------------
Q  = find( md  &  F>0  &  F<1 );
if isempty(Q), return, end
if xa(1), FQ=F(Q); FQ=FQ(:); else FQ=F*ones(length(Q),1); end
if xa(2), vQ=v(Q); vQ=vQ(:); else vQ=v*ones(length(Q),1); end
if xa(3), wQ=w(Q); wQ=wQ(:); else wQ=w*ones(length(Q),1); end

%-Interval bisection
%--------------------------------------------------------------------------
a  = zeros(length(Q),1); fa = a-FQ;
b  =  ones(length(Q),1); fb = b-FQ;
i  = 0;
xQ = a+1/2;
QQ = 1:length(Q);

while ~isempty(QQ) &&  i<maxIt
    fxQQ        = betainc(xQ(QQ),vQ(QQ),wQ(QQ))-FQ(QQ);
    mQQ         = fa(QQ).*fxQQ > 0;
    a(QQ(mQQ))  = xQ(QQ(mQQ));   fa(QQ(mQQ))  = fxQQ(mQQ);
    b(QQ(~mQQ)) = xQ(QQ(~mQQ));  fb(QQ(~mQQ)) = fxQQ(~mQQ);
    xQ(QQ)      = a(QQ) + (b(QQ)-a(QQ))/2;
    QQ          = QQ( (b(QQ)-a(QQ))>tol );
    i           = i+1;
end

if i==maxIt
    warning('convergence criteria not reached - maxIt reached');
end

x(Q) = xQ;

author: MBB-team

repo: VBA-toolbox


text:
% DEMSWISSROLLLLE4 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('swissRoll');

options = lleOptions(32);
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model, 2);

lvmScatterPlotColor(model, model.Y(:, 2));

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, model.Y(:, 2), 'SwissRoll', 4, true);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function plot_comparison(P1,P2,fig,ms1,ms2)

% Plots two point clouds "P1" and "P2" so that those points of "P2" which are
% not in "P1" are plotted in red whereas the common points are plotted in
% blue. "fig" and "ms1" and "ms2" are the figure number and marker sizes.

if nargin == 3
    ms1 = 3;
    ms2 = 3;
elseif nargin == 4
    ms2 = 3;
end

if ms1 == 0
    ms1 = 3;
end
if ms2 == 0
    ms2 = 3;
end

P2 = setdiff(P2,P1,'rows');

figure(fig)
if size(P1,2) == 3
    plot3(P1(:,1),P1(:,2),P1(:,3),'.b','Markersize',ms1)
    hold on
    plot3(P2(:,1),P2(:,2),P2(:,3),'.r','Markersize',ms2)
elseif size(P1,2) == 2
    plot(P1(:,1),P1(:,2),'.b','Markersize',ms1)
    hold on
    plot(P2(:,1),P2(:,2),'.r','Markersize',ms2)
end
hold off
axis equal

author: InverseTampere

repo: TreeQSM


text:
function b = diag(a,k)
% function B = diag(A,K)
%
% DESCRIPTION
%   Diagonal polynomial matrices and diagonals of polynomial matrices.
%
% INPUTS
%   A: polynomial
%   K: integer
%
% OUTPUTS
%   B: polynomial
%
% SYNTAX
%   B=diag(A,K)
%     If A is a polynomial vector of length N then B is an M-by-M
%     square matrix with M=N+ABS(K) and the elements of A on the K^th
%     diagonal.  K=0 is the main diagonal, K>0 is above the main diagonal
%     and K<0 is below the main diagonal
%   B=diag(A,K)
%     If A is a polynomial matrix then B is a column vector whose elements
%     are the entries of the K^th diagonal of A.
%   B=diag(A)
%     This is the same as B=diag(A,0).
%
% See also triu, tril, blkdiag

% 12/3/2009: PJS  Initial Coding

if nargin==1
    k=0;
end

sza = size(a);
if any(sza==1)
    % a is a vector --> b is a matrix

    % Use double/diag to find right dimensions and indices
    la = length(a);
    bmat = diag(1:la,k);
    idx = find(bmat);

    % Assign values
    b = polynomial( zeros(size(bmat)) );
    L.type = '()';
    L.subs = {idx};
    b = subsasgn(b,L,a);
else
    % a is a matrix --> b is a vector
    if k>=sza(2) || k <= -sza(1)
        b = polynomial(zeros(0,1));
        return
    end

    % Use double/diag to find right indices
    amat = reshape( 1:(sza(1)*sza(2)), sza );
    idx = diag(amat,k);

    % Assign values
    L.type = '()';
    L.subs = {idx};
    b = subsref(a,L);
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function g = lfmKernGradient(kern, t, varargin)

% LFMKERNGRADIENT Gradient of LFM kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the latent force
% model kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in t.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG t1 : the input locations associated with the rows of the
% kernel matrix.
% ARG t2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as t1 and the same number of columns
% as t2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO lfmKernParamInit, kernGradient, lfmKernDiagGradient, kernGradX,
% lfmXlfmKernGradient
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

if length(varargin)<2
  t2 = t;
else
  t2 = varargin{1};
end

[g1, g2] = lfmXlfmKernGradient(kern, kern, t, t2, varargin{end});

g = real(g1 + g2);

%g = real(g1);


text:
classdef TestAgastFeatureDetector
    %TestAgastFeatureDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            obj = cv.AgastFeatureDetector('Threshold',10, 'NonmaxSuppression',true);
            assert(obj.Threshold == 10);
            assert(obj.NonmaxSuppression == true);
            typename = obj.typeid();

            img = imread(TestAgastFeatureDetector.im);
            kpts = obj.detect(img);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpts))));
        end

        function test_detect_imgset
            img = imread(TestAgastFeatureDetector.im);
            imgs = {img, img};
            obj = cv.AgastFeatureDetector();

            kpts = obj.detect(imgs);
            validateattributes(kpts, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, {'vector'}), kpts);
            cellfun(@(kpt) assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpt)))), kpts);
        end

        function test_detect_mask
            img = imread(TestAgastFeatureDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.AgastFeatureDetector();
            kpts = obj.detect(img, 'Mask',mask);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestAgastFeatureDetector.kfields, fieldnames(kpts))));
            xy = cat(1, kpts.pt);
            assert(all(xy(:,1) <= ceil(size(img,2)/2)));
        end

        function test_detect_types
            types = {'AGAST_5_8', 'AGAST_7_12d', 'AGAST_7_12s', 'OAST_9_16'};
            img = imread(TestAgastFeatureDetector.im);
            obj = cv.AgastFeatureDetector();
            for i=1:numel(types)
                obj.Type = types{i};
                kpts = obj.detect(img);
                validateattributes(kpts, {'struct'}, {'vector'});
            end
        end

        function test_error_1
            try
                cv.AgastFeatureDetector('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Construction Data_MI
% MI: 50 trials
% 4-class (Cylindrical/Spherical/LeftTwist/RightTwist)
%% Initializing
clc; close all; clear all;
%% file
dd='C:\Users\Doyeunlee\Desktop\Journal_dylee\1_ConvertedData\';
fsldx=2;
fs={'100','250','1000'};
ref='FCz';
channelMatrix={'F3','F1','Fz','F2','F4';
    'FC3','FC1','FCz','FC2','FC4';
    'C3','C1', 'Cz', 'C2', 'C4';
    'CP3','CP1','CPz','CP2','CP4';
    'P3','P1','Pz','P2','P4'};

trainRatio=0.8;
modldx=1;
modality={'MI','realMove'};

filelist={'sub01','sub02','sub03','sub04','sub05'};

%%
for sub = 1:length(filelist)
    ival = [0 3000];

    % Data load - MI
    [cntGrasp,mrkGrasp,mntGrasp]=eegfile_loadMatlab([dd 'realMove' '\' fs{fsldx} '\' filelist{sub} '_multigrasp_' 'realMove']);
    [cntTwist,mrkTwist,mntTwist]=eegfile_loadMatlab([dd 'realMove' '\' fs{fsldx} '\' filelist{sub} '_twist_' 'realMove']);

    % theta, mu, and beta band
    filterBank = {'[4 8]','[8 13]','[13 40]'};

    for filt = 1:length(filterBank)
        filtBank = {filt};

        % Butterworth
        cntGrasp=proc_filtButter(cntGrasp,5,[4 40]);
        cntTwist=proc_filtButter(cntTwist,5,[4 40]);

        % CAR
        cntGrasp = proc_commonAverageReference(cntGrasp);
        cntTwist = proc_commonAverageReference(cntTwist);

        epoGrasp=cntToEpo(cntGrasp,mrkGrasp,ival);
        epoTwist=cntToEpo(cntTwist,mrkTwist,ival);

        epoGrasp=proc_selectChannels(epoGrasp,{'F3','F1','Fz','F2','F4',...
            'FC3','FC1','FCz','FC2','FC4',...
            'C3','C1', 'Cz', 'C2', 'C4', ...
            'CP3','CP1','CPz','CP2','CP4',...
            'P3','P1','Pz','P2','P4'});
        epoTwist=proc_selectChannels(epoTwist,{'F3','F1','Fz','F2','F4',...
            'FC3','FC1','FCz','FC2','FC4',...
            'C3','C1', 'Cz', 'C2', 'C4', ...
            'CP3','CP1','CPz','CP2','CP4',...
            'P3','P1','Pz','P2','P4'});

        % Retrieve trial amounts
        %% 4-class: Forward, Cylindrical, LeftTwist, Rest
        epoCylindrical=proc_selectClasses(epoGrasp,{'Cylindrical'});
        epoSpherical=proc_selectClasses(epoGrasp,{'Spherical'});
        epoLeftTwist=proc_selectClasses(epoTwist,{'LeftTwist'});
        epoRightTwist=proc_selectClasses(epoTwist,{'RightTwist'});

        % append
        epo = proc_appendEpochs(epoCylindrical,epoSpherical);
        epo = proc_appendEpochs(epo,epoLeftTwist);
        epo = proc_appendEpochs(epo,epoRightTwist);

        % Data shuffling
        epoCylindrical.x=datasample(epoCylindrical.x,size(epoCylindrical.x,3),3,'Replace',false);
        epoSpherical.x=datasample(epoSpherical.x,size(epoSpherical.x,3),3,'Replace',false);
        epoLeftTwist.x=datasample(epoLeftTwist.x,size(epoLeftTwist.x,3),3,'Replace',false);
        epoRightTwist.x=datasample(epoRightTwist.x,size(epoRightTwist.x,3),3,'Replace',false);

        % CSP
        [csp_fv,csp_w,csp_eig]=proc_multicsp(epo,3);
        proc=struct('memo','csp_w');

        proc.train=['[fv,csp_w]= proc_multicsp(fv, 3);' ...
            'fv=proc_variance(fv);' ...
            'fv=proc_logarithm(fv);'];

        proc.apply=['fv=proc_linearDerivation(fv, csp_w);','fv=proc_variance(fv);','fv=proc_logarithm(fv);'];

        % Setting training Data
        minTrial=min([size(epoCylindrical.x,3), ...
            size(epoSpherical.x,3), ...
            size(epoLeftTwist.x,3), size(epoRightTwist.x,3)]);
        trainNum=minTrial*trainRatio;

        % training
        epoCylindricalTrain=proc_selectEpochs(epoCylindrical,1:trainNum);
        epoSphericalTrain=proc_selectEpochs(epoSpherical,1:trainNum);
        epoLeftTwistTrain=proc_selectEpochs(epoLeftTwist,1:trainNum);
        epoRightTwistTrain=proc_selectEpochs(epoRightTwist,1:trainNum);

        % test
        epoCylindricalTest=proc_selectEpochs(epoCylindrical,trainNum+1:size(epoCylindrical.x,3));
        epoSphericalTest=proc_selectEpochs(epoSpherical,trainNum+1:size(epoSpherical.x,3));
        epoLeftTwistTest=proc_selectEpochs(epoLeftTwist,trainNum+1:size(epoLeftTwist.x,3));
        epoRightTwistTest=proc_selectEpochs(epoRightTwist,trainNum+1:size(epoRightTwist.x,3));

        % Concat
        epoTrain=proc_appendEpochs(epoCylindricalTrain,epoSphericalTrain);
        epoTrain=proc_appendEpochs(epoTrain,epoLeftTwistTrain);
        epoTrain=proc_appendEpochs(epoTrain,epoRightTwistTrain);
        Train=epoTrain;

        epoTest=proc_appendEpochs(epoCylindricalTest,epoSphericalTest);
        epoTest=proc_appendEpochs(epoTest,epoLeftTwistTest);
        epoTest=proc_appendEpochs(epoTest,epoRightTwistTest);
        Test=epoTest;

        %%
        trainX=epoTrain.x;
        testX=epoTest.x;
        trainY=epoTrain.y;testY=epoTest.y;
        trainY=vec2ind(trainY)-1;
        testY=vec2ind(testY)-1;
        str= char('_train');
        str2=char('_test');
        str3=char('.mat');
%         str4=char();
        save(['C:\Users\Doyeunlee\Desktop\AAAI_20190717\Dataconstruction_filter\realMove_train/',num2str(sub),filterBank{filt},str3],'trainX','trainY');
        save(['C:\Users\Doyeunlee\Desktop\AAAI_20190717\Dataconstruction_filter\realMove_test/',num2str(sub),filterBank{filt},str3],'testX','testY');
    end % filterBank
end % subject

disp('Done');

author: PatternRecognition

repo: OpenBMI


text:
function []=Main(result,speed,options,problem,sgd)
global mones  mzeros convert  in_size gate_size numMmcell globalData globalMask usegpu...
    out_size share_size share_size2  in ingate cellstate cells outgate...
    node_outgateInit cellinInit node_cellbiasInit delta_outInit cellstatusInit signum

 signum = 3.5;
if problem.batchsize>problem.numsamples
    problem.batchsize=problem.numsamples;
end

[mones,mzeros,convert,usegpu]=gputype(speed.usegpu);
[data,mask,test,masktest ]=feval(['gen' problem.name],problem);
problem.numsamples=size(data,1);

[in_size,gate_size,out_size,share_size,share_size2,numMmcell,W,in,ingate,cellstate,cells,outgate,...
    node_outgateInit,cellinInit,node_cellbiasInit,delta_outInit,cellstatusInit ]=netInit(problem);
if problem.continualPredict
    Bpfunc=@batch_equal_nomask_lstm;
else
    Bpfunc=@batch_cell_lstm;
end

if speed.gradientchecking

    for i=1:fix(problem.numsamples/problem.batchsize)

        bashdata{i} =  data(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
    end

    globalData=bashdata{1};
    globalMask=maskdata{1};
    [ err,dw,inLL,right ]=Bpfunc(W);

    [numgrad]=computeNumericalGradient(Bpfunc,W);
    disp( [numgrad dw])
    diff = norm(numgrad-dw)/norm(numgrad+dw);
    disp(diff);
    fprintf('Norm of the difference between numerical and analytical gradient (should be < 1e-9)\n\n');
    return
end

if speed.usecluster
    serveroptions.slavecount=speed.numcluster;
    serveroptions.slavedir='/usr/local/anew-lstm';
    server=Server(serveroptions);

    markers = floor(linspace(1, size(data,2), server.slaveCount+1));
    markers(end) = size(data,2);
    markers(1) = 0;
    trainStarts = cell(1,server.slaveCount);
    trainEnds =cell(1,server.slaveCount);
    for i=1:server.slaveCount
        trainStarts{i} = markers(i)+1;
        trainEnds{i} = markers(i+1);
        fprintf('i = %d, trainStart = %d, trainEnd = %d\n', i, trainStarts{i}, trainEnds{i});
    end

    options.Method = 'lbfgs';

    for i=1:fix(problem.numsamples/problem.batchsize/server.slaveCount )
        bashdata{i} =  data(  (i-1)*problem.batchsize*server.slaveCount  + 1:...
            i*problem.batchsize*server.slaveCount  , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize*server.slaveCount + 1:...
            i*problem.batchsize*server.slaveCount  , : , :   );
    end
    clear data

    record=1;
    timer=tic;
    for epoch =1:1
        for i=1:fix(problem.numsamples/problem.batchsize/server.slaveCount )
            [Winit]=  server.rpc('clientLoadDataMinibatchNomask_ref',...
                trainStarts, trainEnds, problem.batchsize  ,...
                problem,speed,bashdata{i},maskdata{i},{1:server.slaveCount});
            if i==1&&epoch==1
                W=Winit{1};
            end
            for repeat = 1:32
                [W, cost] = minFunc( @server_batch_cell_lstm,W,options,server);
                if  toc(timer)>1
                    timearray{1}( record )= toc(timer)
                    [errorarray{1}(record),~,~,rightarray{1}(record) ]=testmodel(W,test,masktest ,problem);
                    timer=tic;
                    record=record+1;
                end
            end
        end
    end
    figure
    for method = 1:1
        plot(cumsum(timearray{method}), rightarray{method} , result.linspec{method}); hold on
    end
    legend('server  ' );
    xlabel('time(seconds)');
    ylabel('test obj');

    saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest) problem.name '.fig']);
    saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps'],'epsc');
    eps2pdf([pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps']);

    save('server.mat','timearray','errorarray');

elseif ~speed.usecluster

    for i=1:fix(problem.numsamples/problem.batchsize)

        bashdata{i} =  data(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
        maskdata{i}= mask(  (i-1)*problem.batchsize + 1:i*problem.batchsize , : , :   );
    end

    Winit=W;

   timer=tic;
   record=1;
    for i=1:fix(problem.numsamples/problem.batchsize)
        for inrepeat=1:20

            globalData =  bashdata{i};
            globalMask= maskdata{i};

            options.Method = 'lbfgs';

            [W, cost] = minFunc( Bpfunc,W,options);
            if toc(timer)>2

                timearray{1}( record  )= toc(timer);
                [errorarray{1}(record),~,~,rightarray{1}(record) ]=testmodel(W,test,masktest ,problem);
                disp([ 'lbfgs error '  num2str(errorarray{1}(record) )]  )
                record=record+1;
                timer=tic;
            end
        end
    end

        W=Winit;
        momentum=sgd.momentum;
        alpha=sgd.alpha;
        oldGradient=0;

        record =1;
        timer=tic;
        for repeat =1 :5
            for inner =1:50

                for  i=1:fix(problem.numsamples/problem.batchsize)

                    globalData =  bashdata{i};
                    globalMask= maskdata{i};

                    [~,dw,~,~]=Bpfunc(W);
                    oldGradient = alpha*dw + momentum* oldGradient;
                    W= W - alpha* oldGradient;
                end
                timearray{2}(record  )= toc(timer);
                [errorarray{2}(record),~,~, rightarray{2}(record) ]=testmodel(W,test,masktest ,problem);
                record=record+1;
                timer=tic;
            end
            disp([ 'sgd error '  num2str(errorarray{2}(record-1) )]  )
        end

    figure
    for method = 1:2
        plot(cumsum(timearray{method}), rightarray{method} , result.linespec{method}); hold on
    end
    legend('lbfgs','sgd');
    xlabel('time(seconds)');
    ylabel('test obj');

    figure
    for method = 1:2
        plot(cumsum(timearray{method}),errorarray{method} , result.linespec{method}); hold on
    end
    legend('lbfgs','sgd');
    xlabel('time(seconds)');
    ylabel('test e');

    if result.savepic
        saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest) problem.name '.fig']);
        saveas(gcf,[ pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps'],'epsc');
     %   eps2pdf([pwd '/picture/' num2str(problem.numsamples) '_' num2str(problem.Ttest)  problem.name '.eps']);
    end

end

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function colorDetectHSV(fileName, hsvVal, tol)

%
% function colorDetectHSV(fileName, hsvVal, tol)
%
% This function is used for detecting a specified hsv value in images.
%
% ARGUMENTS:
% fileName: the name of the jpg file to be loaded
% hsvVal: 3x1 array containing the HSV values to be detected
% tol: 1x1 or 2x1 or 3x1 array containing the tolerance (i.e. the maximum
% distance - in each hsv coefficient - of each pixel from hsvVal).
%
% Example:
% colorDetectHSV('train/face07.jpg', median(HSV), [0.05 0.05 0.1]);
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theodoros Giannakopoulos - January 2008
% www.di.uoa.gr/~tyiannak
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RGB = imread(fileName);

HSV = rgb2hsv(RGB);

% find the difference between required and real H value:
diffH = abs(HSV(:,:,1) - hsvVal(1));

[M,N,t] = size(RGB);
I1 = zeros(M,N); I2 = zeros(M,N); I3 = zeros(M,N);

T1 = tol(1);

I1( find(diffH < T1) ) = 1;

if (length(tol)>1)
    % find the difference between required and real S value:
    diffS = abs(HSV(:,:,2) - hsvVal(2));
    T2 = tol(2);
    I2( find(diffS < T2) ) = 1;
    if (length(tol)>2)
        % find the difference between required and real V value:
        difV = HSV(:,:,3) - hsvVal(3);
        T3 = tol(3);
        I3( find(diffS < T3) ) = 1;
        I = I1.*I2.*I3;
    else
        I = I1.*I2;
    end
else
    I = I1;
end

subplot(1,2,1),imshow(RGB); title('Original Image');
subplot(1,2,2),imshow(I,[]); title('Detected Areas');


text:
function x = solvedbi_sm(ah, rho, b, c)

% solvedbi_sm -- Solve a diagonal block linear system with a scaled identity
%                term using the Sherman-Morrison equation
%
%         The solution is obtained by independently solving a set of linear
%         systems of the form (see wohlberg-2016-efficient)
%
%                  (a a^H + rho I) x = b
%
%         In this equation inner products and matrix products are taken along
%         the 3rd dimension of the corresponding multi-dimensional arrays; the
%         solutions are independent over the 1st and 2nd (and 4th, if
%         non-singleton) dimensions.
%
% Usage:
%       x = solvedbi_sm(ah, rho, b, c);
%
% Input:
%       ah          Multi-dimensional array containing a^H
%       rho         Scalar rho
%       b           Multi-dimensional array containing b
%       c           Multi-dimensional array containing pre-computed quantities
%                   a^H / (a^H a + rho)
%
% Output:
%       x           Multi-dimensional array containing linear system solution
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2014-12-18
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

a = conj(ah);

if nargin < 4 || isempty(c),
  c = bsxfun(@rdivide, ah, sum(ah.*a, 3) + rho);
end

cb = sum(bsxfun(@times, c, b), 3);
clear ah c;
cba = bsxfun(@times, cb, a);
clear a cb;
x = (b - cba) / rho;

return

author: thfylsty


text:
%% flandmark_learn.m
% TODO - add description
%
% 14-07-10 Michal Uricar
% 17-05-11 Michal Uricar
% 02-08-11 Michal Uricar, flandmark2 - detection of corners of eyes and mouth
% 21-03-12 Michal Uricar, LFW annotation in one file
% 13-09-12 Michal Uricar, refactoring

clc;
close all; clearvars;

%% User defined options

% setup paths to face database (and benchmarks)
ROOT_FACE_FOLDER = '../LFW/';                           % path to folder with database images
DB_ANN = '../data/LFW_annotation.txt';                  % annotation file for the database
%TRN_EXAMPLES = '../benchmarks/LFW_TRN_annotation.txt';  % benchmark db division, leave empty to create different division
%VAL_EXAMPLES = '../benchmarks/LFW_VAL_annotation.txt';  % benchmark db division, leave empty to create different division
%TST_EXAMPLES = '../benchmarks/LFW_TST_annotation.txt';  % benchmark db division, leave empty to create different division
TRN_EXAMPLES = [];
VAL_EXAMPLES = [];
TST_EXAMPLES = [];

% options configuration
% order of components in LFW database
options.compnames = { 'center', 'canthus-rl', 'canthus-lr', 'mouth-corner-r', 'mouth-corner-l', ...
                      'canthus-rr', 'canthus-ll', 'nose' };
options.bw = [40; 40];
options.components = [20 10 10 10 10 10 10 10;
                      20 10 10 10 10 10 10 10];
options.comselect = [10, 4, 5, 8, 9, 3, 6, 10];
% number of components
options.M = size(options.components, 2);
% base window margin (percents)
options.bw_margin = [20; 20];
% path to LFW database (where images and anotations are stored) must contain subdir /mgt_v2 with xml annotations
options.image_path = ROOT_FACE_FOLDER;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Add paths

	addpath('./code/');

%% Check project dir structure

    fprintf('Checking project dir structure...\n');

    DIRS = {'./code/img/', './code/results/', './code/MAT'};

    for i = 1 : numel(DIRS)
        fprintf('%s ... ', DIRS{i});
        if (~exist(DIRS{i}, 'dir'))
            fprintf('MISSING\n');
            fprintf('Creating dir %s...\n', DIRS{i});
            mkdir(DIRS{i});
        else
            fprintf('OK\n');
        end;
    end;

    DIRS = {'./bmrm', './libqp', './code/Functions', ...
            ROOT_FACE_FOLDER};

    for i = 1 : numel(DIRS)
        fprintf('%s ... ', DIRS{i});
        if (~exist(DIRS{i}, 'dir'))
            fprintf('MISSING\n');
            error('??? Corrupted project structure.');
        else
            fprintf('OK\n');
        end;
    end;

%% Save above specified user data

    options.image_path = ['../' ROOT_FACE_FOLDER];

    save('./code/MAT/paths.mat', 'ROOT_FACE_FOLDER', 'TRN_EXAMPLES', 'VAL_EXAMPLES', 'TST_EXAMPLES', 'DB_ANN');
    save('./code/MAT/options.mat', 'options');

%% Compile mex files if neccessary

    fprintf('Checking for neccessary mex-files...\n');

    if (~exist(['../matlab_toolbox/mex/lbp2sparse.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/argmax_mex.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/getNormalizedFrame.' mexext], 'file') || ...
        ~exist(['../matlab_toolbox/mex/flandmark_detector.' mexext], 'file'))
        error(['??? Mex-files not found! Please compile the flandmark with BUILD_MATLAB_BINDINGS turned on'...
               ' or run /matlab_toolbox/flandmark_compilemex.m to create them.\n']);
    end

    fprintf('Done.\n');

%% Collect statistics about components positions

    % check if statistics already exist
    if (exist('./code/results/statistics_orig.mat', 'file'))
        fprintf('Statistics found. Skipping...\n');
    else
        fprintf('Collecting statistics about components...\n');

        run './code/runstats.m';
    end;

%% Division of database into 3 parts (60% traning, 20% validating, 20% testing)

    % check if division already done.
    if (exist('./code/MAT/TRN.mat', 'file') && ...
            exist('./code/MAT/VAL.mat', 'file') && ...
            exist('./code/MAT/TST.mat', 'file'))
        fprintf('Division already done. Skipping...\n');
    else
        fprintf('Dividing database into 3 parts...\n');

        run './code/dbDivision.m';
    end;

%% Create data for BMRM solver (TRN, VAL, TST)

    %% Functions selection

    fprintf('Prepare functions configurations\n');
    f_argmax = @argmax2;
    f_risk = @risk2;
%     f_risk = @risk2_unnormalized;
    f_buildPsi = @buildPsi2;
    f_createMapTable = @createMapTable2;
    displacement = true;
    save('./code/MAT/functions_gdisp.mat', 'f_argmax', 'f_risk', ...
        'f_buildPsi', 'f_createMapTable', 'displacement');

    %% Configuration - nImages, height_of_pyramid
    fprintf('Looking for config file...\n');
    if (~exist('./code/MAT/config.mat', 'file'))
        fprintf('Not found. Creating...\nconfig = \n');
        config.nImages = 'all';
        config.heightOfPyramid = 4;
        disp(config);
        save('./code/MAT/config.mat', 'config');
    else
        fprintf('File %s found.\n', './code/MAT/config.mat');
    end;

    fprintf('Prepare data for LBP\n');
    run './code/prepDataForLBP.m';

%% Call BMRM solver

    % prepare lamdas
    fprintf('Looking for file with lambda values...\n');
    if (~exist('./code/MAT/lambdas.mat', 'file'))
        lambdas = [1e-3, 1e-2, 1e-1, 1];
        fprintf('Not found. Creating...\n');
        save('./code/MAT/lambdas.mat', 'lambdas');
        disp(lambdas);
    else
        fprintf('Lambda values found:');
        load('./code/MAT/lambdas.mat');
        disp(lambdas);
    end;

    % prepare normalization coefficients for each image in database
    if (~exist('./code/MAT/kappas.mat', 'file'))
	run('./code/computeKappas.m');
    end;

    fprintf('Calling BMRM solver...\n');
    run './code/test_bmrm.m';

%% Validation

    if (exist('./code/results/W.mat', 'file'))
        fprintf('Best vector W already found. Skipping...\n');
    else
        fprintf('Running validaton...\n');

        run './code/testValidation.m';
    end;

%% Test

    if (exist('./code/results/errors.mat', 'file'))
        fprintf('Test already passed. Skipping...\n');
    else
        fprintf('Running test...\n');

        run './code/test1.m';
    end;

    fprintf('Show histograms.\n');
    run './code/test2.m';

%% Generate model.mat (for use with flandmark_load_model - binary dat file creation)

    run './code/detectorModelBuild.m';

%% Timestamp

 fprintf(1,'Finished on %s.\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function vl_demo_kdtree
% VL_DEMO_KDTREE
%   Demonstrates the construction of a kd-tree for uniformly
%   distributed 2-D data and data distributed on a circle. Both the
%   median and mean threshold selection methods are illustrated.

  randn('state',0) ;
  rand('state',0) ;

  % uniform 2D points
  Q = single(rand(2,1)) ;
  X = single(rand(2, 100)) ;

  % kd-tree with median and mean thresholds
  figure(1) ; clf ; do('thresholdmethod', 'mean') ;
  figure(2) ; clf ; do('thresholdmethod', 'median') ;

  % 2D points on a circle
  thr = 2*pi*rand(1,100) ;
  X = single(.4*[cos(thr); sin(thr)] + .5) ;
  Q = single(Q) ;

  % kd-tree with median and mean thresholds
  figure(3) ; clf ; do('thresholdmethod', 'mean') ;
  figure(4) ; clf ; do('thresholdmethod', 'median') ;

  vl_demo_print(1, 'kdtree_uniform_mean') ;
  vl_demo_print(2, 'kdtree_uniform_median') ;
  vl_demo_print(3, 'kdtree_circle_mean') ;
  vl_demo_print(4, 'kdtree_circle_median') ;

  function do(varargin)
    % build the forest
    kdforest = vl_kdtreebuild (X, 'verbose', varargin{:}) ;

    % query the kd-tree
    [i, d] = vl_kdtreequery (kdforest, X, Q, 'numneighbors', 10, 'verbose') ;
    vl_plotframe(X, 'ro') ;
    hold on ;
    xl = [0, +1] ;
    yl = [0, +1] ;
    axis equal ;
    xlim(xl) ;
    ylim(yl) ;
    vl_demo_kdtree_plot(kdforest.trees(1), 1, xl, yl) ;

    % plot query point and returns
    vl_plotframe(Q,'b*') ;
    for k=1:length(i)
      vl_plotframe([Q ; sqrt(d(k))],'b:','linewidth',1) ;
      vl_plotframe(X(:, i(k)), 'bx') ;
    end
    set(gca,'xtick',[],'ytick',[]) ;
  end

end

author: Cloud-CV

repo: object-proposals