text:
function [t,sts] = cfg_getfile(varargin)
% File selector
% FORMAT [t,sts] = cfg_getfile(n,typ,mesg,sel,wd,filt,frames)
%     n    - Number of files
%            A single value or a range.  e.g.
%            1       - Select one file
%            Inf     - Select any number of files
%            [1 Inf] - Select 1 to Inf files
%            [0 1]   - select 0 or 1 files
%            [10 12] - select from 10 to 12 files
%     typ  - file type
%           'any'   - all files
%           'batch' - SPM batch files (.m, .mat and XML)
%           'dir'   - select a directory
%           'image' - Image files (".img" and ".nii")
%                     Note that it gives the option to select
%                     individual volumes of the images.
%           'mat'   - Matlab .mat files or .txt files (assumed to contain
%                     ASCII representation of a 2D-numeric array)
%           'mesh'  - Mesh files (".gii" and ".mat")
%           'nifti' - NIfTI files without the option to select frames
%           'xml'   - XML files
%           Other strings act as a filter to regexp.  This means
%           that e.g. DCM*.mat files should have a typ of '^DCM.*\.mat$'
%      mesg - a prompt (default 'Select files...')
%      sel  - list of already selected files
%      wd   - Directory to start off in
%      filt - value for user-editable filter (default '.*')
%      frames - Image frame numbers to include (default '1')
%
%      t    - selected files
%      sts  - status (1 means OK, 0 means window quit)
%
% FORMAT [t,ind] = cfg_getfile('Filter',files,typ,filt,frames)
% filter the list of files (cell array) in the same way as the
% GUI would do. There is an additional typ 'extimage' which will match
% images with frame specifications, too. The 'frames' argument
% is currently ignored, i.e. image files will not be filtered out if
% their frame numbers do not match.
% When filtering directory names, the filt argument will be applied to the
% last directory in a path only.
% t returns the filtered list (cell array), ind an index array, such that
% t = files(ind).
%
% FORMAT cpath = cfg_getfile('CPath',path,cwd)
% function to canonicalise paths: Prepends cwd to relative paths, processes
% '..' & '.' directories embedded in path.
% path     - string matrix containing path name
% cwd      - current working directory [default '.']
% cpath    - conditioned paths, in same format as input path argument
%
% FORMAT [files,dirs]=cfg_getfile('List',direc,filt)
% Returns files matching the filter (filt) and directories within dire
% direc    - directory to search
% filt     - filter to select files with (see regexp) e.g. '^w.*\.img$'
% files    - files matching 'filt' in directory 'direc'
% dirs     - subdirectories of 'direc'
% FORMAT [files,dirs]=cfg_getfile('ExtList',direc,filt,frames)
% As above, but for selecting frames of 4D NIfTI files
% frames   - vector of frames to select (defaults to 1, if not
%            specified). If the frame number is Inf, all frames for the
%            matching images are listed.
% FORMAT [files,dirs]=cfg_getfile('FPList',direc,filt)
% FORMAT [files,dirs]=cfg_getfile('ExtFPList',direc,filt,frames)
% As above, but returns files with full paths (i.e. prefixes direc to each)
% FORMAT [files,dirs]=cfg_getfile('FPListRec',direc,filt)
% FORMAT [files,dirs]=cfg_getfile('ExtFPListRec',direc,filt,frames)
% As above, but returns files with full paths (i.e. prefixes direc to
% each) and searches through sub directories recursively.
%
% FORMAT cfg_getfile('prevdirs',dir)
% Add directory dir to list of previous directories.
% FORMAT dirs=cfg_getfile('prevdirs')
% Retrieve list of previous directories.
%
% This code is based on the file selection dialog in SPM5, with virtual
% file handling turned off.
%____________________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% John Ashburner and Volkmar Glauche
% $Id: cfg_getfile.m 6251 2014-10-30 13:58:44Z volkmar $

t = {};
sts = false;
if nargin > 0 && ischar(varargin{1})
    switch lower(varargin{1})
        case {'addvfiles', 'clearvfiles', 'vfiles'}
            cfg_message('matlabbatch:deprecated:vfiles', ...
                        'Trying to use deprecated ''%s'' call.', ...
                        lower(varargin{1}));
        case 'cpath'
            cfg_message(nargchk(2,Inf,nargin,'struct'));
            t = cpath(varargin{2:end});
            sts = true;
        case 'filter'
            filt    = mk_filter(varargin{3:end});
            t       = varargin{2};
            if numel(t) == 1 && isempty(t{1})
                sts = 1;
                return;
            end;
            t1 = cell(size(t));
            if any(strcmpi(varargin{3},{'dir','extdir'}))
                % only filter last directory in path
                for k = 1:numel(t)
                    t{k} = cpath(t{k});
                    if t{k}(end) == filesep
                        [p n] = fileparts(t{k}(1:end-1));
                    else
                        [p n] = fileparts(t{k});
                    end
                    if strcmpi(varargin{3},'extdir')
                        t1{k} = [n filesep];
                    else
                        t1{k} = n;
                    end
                end
            else
                % only filter filenames, not paths
                for k = 1:numel(t)
                    [p n e] = fileparts(t{k});
                    t1{k} = [n e];
                end
            end
            [t1,sts1] = do_filter(t1,filt.ext);
            [t1,sts2] = do_filter(t1,filt.filt);
            sts = sts1(sts2);
            t = t(sts);
        case {'list', 'fplist', 'extlist', 'extfplist'}
            if nargin > 3
                frames = varargin{4};
            else
                frames = 1; % (ignored in listfiles if typ==any)
            end;
            if regexpi(varargin{1}, 'ext') % use frames descriptor
                typ = 'extimage';
            else
                typ = 'any';
            end
            filt    = mk_filter(typ, varargin{3}, frames);
            [t sts] = listfiles(varargin{2}, filt); % (sts is subdirs here)
            sts = sts(~(strcmp(sts,'.')|strcmp(sts,'..'))); % remove '.' and '..' entries
            if regexpi(varargin{1}, 'fplist') % return full pathnames
                direc = cfg_getfile('cpath', varargin{2});
                % remove trailing path separator if present
                direc = regexprep(direc, [filesep '$'], '');
                if ~isempty(t)
                    t = strcat(direc, filesep, t);
                end
                if nargout > 1
                    % subdirs too
                    sts = cellfun(@(sts1)cpath(sts1, direc), sts, 'UniformOutput',false);
                end
            end
        case {'fplistrec', 'extfplistrec'}
            % list directory
            [f1 d1] = cfg_getfile(varargin{1}(1:end-3),varargin{2:end});
            f2 = cell(size(d1));
            d2 = cell(size(d1));
            for k = 1:numel(d1)
                % recurse into sub directories
                [f2{k} d2{k}] = cfg_getfile(varargin{1}, d1{k}, ...
                                            varargin{3:end});
            end
            t = vertcat(f1, f2{:});
            if nargout > 1
                sts = vertcat(d1, d2{:});
            end
        case 'prevdirs',
            if nargin > 1
                prevdirs(varargin{2});
            end;
            if nargout > 0 || nargin == 1
                t = prevdirs;
                sts = true;
            end;
        otherwise
            cfg_message('matlabbatch:usage','Inappropriate usage.');
    end
else
    [t,sts] = selector(varargin{:});
end
%=======================================================================

%=======================================================================
function [t,ok] = selector(n,typ,mesg,already,wd,filt,frames,varargin)
if nargin<1 || ~isnumeric(n) || numel(n) > 2
    n = [0 Inf];
else
    if numel(n)==1,   n    = [n n];    end;
    if n(1)>n(2),     n    = n([2 1]); end;
    if ~isfinite(n(1)), n(1) = 0;        end;
end

if nargin<2 || ~ischar(typ), typ     = 'any';   end;

if nargin<3 || ~(ischar(mesg) || iscellstr(mesg))
    mesg    = 'Select files...';
elseif iscellstr(mesg)
    mesg = char(mesg);
end

if nargin<4 || isempty(already) || (iscell(already) && isempty(already{1}))
    already = {};
else
    % Add folders of already selected files to prevdirs list
    pd1 = cellfun(@(a1)strcat(fileparts(a1),filesep), already, ...
                  'UniformOutput',false);
    prevdirs(pd1);
end

if nargin<5 || isempty(wd) || ~ischar(wd)
    if isempty(already)
        wd      = pwd;
    else
        wd = fileparts(already{1});
        if isempty(wd)
            wd = pwd;
        end
    end;
end

if nargin<6 || ~ischar(filt), filt    = '.*';    end;

if nargin<7 || ~(isnumeric(frames) || ischar(frames))
    frames  = '1';
elseif isnumeric(frames)
    frames = char(gencode_rvalue(frames(:)'));
elseif ischar(frames)
    try
        ev = eval(frames);
        if ~isnumeric(ev)
            frames = '1';
        end
    catch
        frames = '1';
    end
end
ok  = 0;

t = '';
sfilt = mk_filter(typ,filt,eval(frames));

[col1,col2,col3,lf,bf] = colours;

% delete old selector, if any
fg = findobj(0,'Tag',mfilename);
if ~isempty(fg)
    delete(fg);
end
% create figure
fg = figure('IntegerHandle','off',...
    'Tag',mfilename,...
    'Name',mesg,...
    'NumberTitle','off',...
    'Units','Pixels',...
    'MenuBar','none',...
    'DefaultTextInterpreter','none',...
    'DefaultUicontrolInterruptible','on',...
    'Visible','off');

cfg_onscreen(fg);
set(fg,'Visible','on');

sellines = min([max([n(2) numel(already)]), 4]);
[pselp pcntp pfdp pdirp] = panelpositions(fg, sellines+1);

uicontrol(fg,...
          'style','text',...
          'units','normalized',...
          'position',posinpanel([0 sellines/(sellines+1) 1 1/(sellines+1)],pselp),...
          lf,...
          'BackgroundColor',get(fg,'Color'),...
          'ForegroundColor',col3,...
          'HorizontalAlignment','left',...
          'string',mesg,...
          'tag','msg');

% Selected Files
sel = uicontrol(fg,...
    'style','listbox',...
    'units','normalized',...
    'Position',posinpanel([0 0 1 sellines/(sellines+1)],pselp),...
    lf,...
    'Callback',@unselect,...
    'tag','selected',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'Max',10000,...
    'Min',0,...
    'String',already,...
    'Value',1);
c0 = uicontextmenu('Parent',fg);
set(sel,'uicontextmenu',c0);
uimenu('Label','Unselect All', 'Parent',c0,'Callback',@unselect_all);

% get cwidth for buttons
tmp=uicontrol('style','text','string',repmat('X',[1,50]),bf,...
    'units','normalized','visible','off');
fnp = get(tmp,'extent');
delete(tmp);
cw = 3*fnp(3)/50;

if strcmpi(typ,'image'),
    uicontrol(fg,...
        'style','edit',...
        'units','normalized',...
        'Position',posinpanel([0.61 0 0.37 .45],pcntp),...
        'Callback',@update_frames,...
        'tag','frame',...
        lf,...
        'BackgroundColor',col1,...
        'String',frames,'UserData',eval(frames));
% 'ForegroundGolor',col3,...
end;

% Help
uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',posinpanel([0.02 .5 cw .45],pcntp),...
    bf,...
    'Callback',@heelp,...
    'tag','?',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','?',...
    'ToolTipString','Show Help');

uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',posinpanel([0.03+cw .5 cw .45],pcntp),...
    bf,...
    'Callback',@editwin,...
    'tag','Ed',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','Ed',...
    'ToolTipString','Edit Selected Files');

uicontrol(fg,...
   'Style','pushbutton',...
   'units','normalized',...
   'Position',posinpanel([0.04+2*cw .5 cw .45],pcntp),...
   bf,...
   'Callback',@select_rec,...
   'tag','Rec',...
   'ForegroundColor',col3,...
   'BackgroundColor',col1,...
   'String','Rec',...
   'ToolTipString','Recursively Select Files with Current Filter');

% Done
dne = uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',posinpanel([0.05+3*cw .5 0.45-3*cw .45],pcntp),...
    bf,...
    'Callback',@(h,e)delete(h),...
    'tag','D',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','Done',...
    'Enable','off',...
    'DeleteFcn',@null);

if numel(already)>=n(1) && numel(already)<=n(2),
    set(dne,'Enable','on');
end;

% Filter Button
uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',posinpanel([0.51 .5 0.1 .45],pcntp),...
    bf,...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'Callback',@clearfilt,...
    'String','Filt');

% Filter
uicontrol(fg,...
    'style','edit',...
    'units','normalized',...
    'Position',posinpanel([0.61 .5 0.37 .45],pcntp),...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    lf,...
    'Callback',@(ob,ev)update(ob),...
    'tag','regexp',...
    'String',filt,...
    'UserData',sfilt);

% Directories
db = uicontrol(fg,...
    'style','listbox',...
    'units','normalized',...
    'Position',posinpanel([0.02 0 0.47 1],pfdp),...
    lf,...
    'Callback',@click_dir_box,...
    'tag','dirs',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'Max',1,...
    'Min',0,...
    'String','',...
    'UserData',wd,...
    'Value',1);

% Files
tmp = uicontrol(fg,...
    'style','listbox',...
    'units','normalized',...
    'Position',posinpanel([0.51 0 0.47 1],pfdp),...
    lf,...
    'Callback',@click_file_box,...
    'tag','files',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'UserData',n,...
    'Max',10240,...
    'Min',0,...
    'String','',...
    'Value',1);
c0 = uicontextmenu('Parent',fg);
set(tmp,'uicontextmenu',c0);
uimenu('Label','Select All', 'Parent',c0,'Callback',@select_all);

% Drives
if strcmpi(computer,'PCWIN') || strcmpi(computer,'PCWIN64'),
    % get fh for lists
    tmp=uicontrol('style','text','string','X',lf,...
        'units','normalized','visible','off');
    fnp = get(tmp,'extent');
    delete(tmp);
    fh = 2*fnp(4); % Heuristics: why do we need 2*
    sz = get(db,'Position');
    sz(4) = sz(4)-fh-2*0.01;
    set(db,'Position',sz);
    uicontrol(fg,...
        'style','text',...
        'units','normalized',...
        'Position',posinpanel([0.02 1-fh-0.01 0.10 fh],pfdp),...
              'HorizontalAlignment','left',...
        lf,...
        'BackgroundColor',get(fg,'Color'),...
        'ForegroundColor',col3,...
        'String','Drive');
    uicontrol(fg,...
        'style','popupmenu',...
        'units','normalized',...
        'Position',posinpanel([0.12 1-fh-0.01 0.37 fh],pfdp),...
        lf,...
        'Callback',@setdrive,...
        'tag','drive',...
        'BackgroundColor',col1,...
        'ForegroundColor',col3,...
        'String',listdrives(false),...
        'Value',1);
end;

[pd,vl] = prevdirs([wd filesep]);

% Previous dirs
uicontrol(fg,...
    'style','popupmenu',...
    'units','normalized',...
    'Position',posinpanel([0.12 .05 0.86 .95*1/3],pdirp),...
    lf,...
    'Callback',@click_dir_list,...
    'tag','previous',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'String',pd,...
    'Value',vl);
uicontrol(fg,...
    'style','text',...
    'units','normalized',...
    'Position',posinpanel([0.02 0 0.10 .95*1/3],pdirp),...
          'HorizontalAlignment','left',...
    lf,...
    'BackgroundColor',get(fg,'Color'),...
    'ForegroundColor',col3,...
    'String','Prev');

% Parent dirs
uicontrol(fg,...
    'style','popupmenu',...
    'units','normalized',...
    'Position',posinpanel([0.12 1/3+.05 0.86 .95*1/3],pdirp),...
    lf,...
    'Callback',@click_dir_list,...
    'tag','pardirs',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'String',pardirs(wd));
uicontrol(fg,...
    'style','text',...
    'units','normalized',...
    'Position',posinpanel([0.02 1/3 0.10 .95*1/3],pdirp),...
          'HorizontalAlignment','left',...
    lf,...
    'BackgroundColor',get(fg,'Color'),...
    'ForegroundColor',col3,...
    'String','Up');

% Directory
uicontrol(fg,...
    'style','edit',...
    'units','normalized',...
    'Position',posinpanel([0.12 2/3 0.86 .95*1/3],pdirp),...
    lf,...
    'Callback',@edit_dir,...
    'tag','edit',...
    'BackgroundColor',col1,...
    'ForegroundColor',col3,...
    'String','');
uicontrol(fg,...
    'style','text',...
    'units','normalized',...
    'Position',posinpanel([0.02 2/3 0.10 .95*1/3],pdirp),...
          'HorizontalAlignment','left',...
    lf,...
    'BackgroundColor',get(fg,'Color'),...
    'ForegroundColor',col3,...
    'String','Dir');

resize_fun(fg);
set(fg, 'ResizeFcn',@resize_fun);
update(sel,wd)
set(fg,'windowstyle', 'modal');

waitfor(dne);
drawnow;
if ishandle(sel),
    t  = get(sel,'String');
    if isempty(t)
        t = {''};
    elseif sfilt.code == -1
        % canonicalise non-empty folder selection
        t = cellfun(@(t1)cpath(t1, pwd), t, 'UniformOutput',false);
    end;
    ok = 1;
end;
if ishandle(fg),  delete(fg); end;
drawnow;
return;
%=======================================================================

%=======================================================================
function apos = posinpanel(rpos,ppos)
% Compute absolute positions based on panel position and relative
% position
apos = [ppos(1:2)+ppos(3:4).*rpos(1:2) ppos(3:4).*rpos(3:4)];
%=======================================================================

%=======================================================================
function [pselp, pcntp, pfdp, pdirp] = panelpositions(fg, sellines)
if nargin == 1
    na = numel(get(findobj(fg,'Tag','selected'),'String'));
    n  = get(findobj(fg,'Tag','files'),'Userdata');
    sellines = min([max([n(2) na]), 4]);
end
lf = cfg_get_defaults('cfg_ui.lfont');
bf = cfg_get_defaults('cfg_ui.bfont');
% Create dummy text to estimate character height
t=uicontrol('style','text','string','Xg','units','normalized','visible','off',lf);
lfh = 1.05*get(t,'extent');
delete(t)
t=uicontrol('style','text','string','Xg','units','normalized','visible','off',bf);
bfh = 1.05*get(t,'extent');
delete(t)
% panel heights
% 3 lines for directory, parent and prev directory list
% variable height for dir/file navigation
% 2 lines for buttons, filter etc
% sellines plus scrollbar for selected files
pselh = sellines*lfh(4) + 1.2*lfh(4);
pselp = [0 0 1 pselh];
pcnth = 2*bfh(4);
pcntp = [0 pselh 1 pcnth];
pdirh = 3*lfh(4);
pdirp = [0 1-pdirh 1 pdirh];
pfdh  = 1-(pselh+pcnth+pdirh);
pfdp  = [0 pselh+pcnth 1 pfdh];

%=======================================================================

%=======================================================================
function null(varargin)
%=======================================================================

%=======================================================================
function omsg = msg(ob,str)
ob = sib(ob,'msg');
omsg = get(ob,'String');
set(ob,'String',str);
if nargin>=3,
    set(ob,'ForegroundColor',[1 0 0],'FontWeight','bold');
else
    set(ob,'ForegroundColor',[0 0 0],'FontWeight','normal');
end;
drawnow;
return;
%=======================================================================

%=======================================================================
function setdrive(ob,varargin)
st = get(ob,'String');
vl = get(ob,'Value');
update(ob,st{vl});
return;
%=======================================================================

%=======================================================================
function resize_fun(fg,varargin)
% do nothing
return;
[pselp pcntp pfdp pdirp] = panelpositions(fg);
set(findobj(fg,'Tag','msg'), 'Position',pselp);
set(findobj(fg,'Tag','pcnt'), 'Position',pcntp);
set(findobj(fg,'Tag','pfd'), 'Position',pfdp);
set(findobj(fg,'Tag','pdir'), 'Position',pdirp);
return;
%=======================================================================

%=======================================================================
function [d,mch] = prevdirs(d)
persistent pd
if ~iscell(pd), pd = {}; end;
if nargin == 0
    d = pd;
else
    if ~iscell(d)
        d = cellstr(d);
    end
    d   = unique(d(:));
    mch = cellfun(@(d1)find(strcmp(d1,pd)), d, 'UniformOutput',false);
    sel = cellfun(@isempty, mch);
    npd = numel(pd);
    pd  = [pd(:);d(sel)];
    mch = [mch{~sel} npd+(1:nnz(sel))];
    d = pd;
end
return;
%=======================================================================

%=======================================================================
function pd = pardirs(wd)
if ispc
    fs = '\\';
else
    fs = filesep;
end
pd1 = textscan(wd,'%s','delimiter',fs,'MultipleDelimsAsOne',1);
if ispc
    pd = cell(size(pd1{1}));
    pd{end} = pd1{1}{1};
    for k = 2:numel(pd1{1})
        pd{end-k+1} = fullfile(pd1{1}{1:k},filesep);
    end
else
    pd = cell(numel(pd1{1})+1,1);
    pd{end} = filesep;
    for k = 1:numel(pd1{1})
        pd{end-k} = fullfile(filesep,pd1{1}{1:k},filesep);
    end
end
%=======================================================================

%=======================================================================
function clearfilt(ob,varargin)
set(sib(ob,'regexp'),'String','.*');
update(ob);
return;
%=======================================================================

%=======================================================================
function click_dir_list(ob,varargin)
vl = get(ob,'Value');
ls = get(ob,'String');
update(ob,deblank(ls{vl}));
return;
%=======================================================================

%=======================================================================
function edit_dir(ob,varargin)
update(ob,get(ob,'String'));
return;
%=======================================================================

%=======================================================================
function c = get_current_char(lb)
fg = sib(lb, mfilename);
c = get(fg, 'CurrentCharacter');
if ~isempty(c)
    % reset CurrentCharacter
    set(fg, 'CurrentCharacter', char(13));
end
%=======================================================================

%=======================================================================
function click_dir_box(lb,varargin)
c = get_current_char(lb);
if isempty(c) || isequal(c,char(13))
    vl  = get(lb,'Value');
    str = get(lb,'String');
    pd  = get(sib(lb,'edit'),'String');
    while ~isempty(pd) && strcmp(pd(end),filesep)
        pd=pd(1:end-1);      % Remove any trailing fileseps
    end
    sel = str{vl};
    if strcmp(sel,'..'),     % Parent directory
        [dr odr] = fileparts(pd);
    elseif strcmp(sel,'.'),  % Current directory
        dr = pd;
        odr = '';
    else
        dr = fullfile(pd,sel);
        odr = '';
    end;
    update(lb,dr);
    if ~isempty(odr)
        % If moving up one level, try to set focus on previously visited
        % directory
        cdrs = get(lb, 'String');
        dind = find(strcmp(odr, cdrs));
        if ~isempty(dind)
            set(lb, 'Value',dind(1));
        end
    end
end
return;
%=======================================================================

%=======================================================================
function re = getfilt(ob)
ob  = sib(ob,'regexp');
ud  = get(ob,'UserData');
re  = struct('code',ud.code,...
             'frames',get(sib(ob,'frame'),'UserData'),...
             'ext',{ud.ext},...
             'filt',{{get(sib(ob,'regexp'),'String')}});
return;
%=======================================================================

%=======================================================================
function update(lb,dr)
lb = sib(lb,'dirs');
if nargin<2 || isempty(dr) || ~ischar(dr)
    dr = get(lb,'UserData');
end;
if ~(strcmpi(computer,'PCWIN') || strcmpi(computer,'PCWIN64'))
    dr    = [filesep dr filesep];
else
    dr    = [dr filesep];
end;
dr(strfind(dr,[filesep filesep])) = [];
[f,d] = listfiles(dr,getfilt(lb));
if isempty(d),
    dr    = get(lb,'UserData');
    [f,d] = listfiles(dr,getfilt(lb));
else
    set(lb,'UserData',dr);
end;
set(lb,'Value',1,'String',d);
set(sib(lb,'files'),'Value',1,'String',f);
set(sib(lb,'pardirs'),'String',pardirs(dr),'Value',1);
[ls,mch] = prevdirs(dr);
set(sib(lb,'previous'),'String',ls,'Value',mch);
set(sib(lb,'edit'),'String',dr);

if numel(dr)>1 && dr(2)==':',
    str = char(get(sib(lb,'drive'),'String'));
    mch = find(lower(str(:,1))==lower(dr(1)));
    if ~isempty(mch),
        set(sib(lb,'drive'),'Value',mch);
    end;
end;
return;
%=======================================================================

%=======================================================================
function update_frames(lb,varargin)
str = get(lb,'String');
%r   = get(lb,'UserData');
try
    r = eval(['[',str,']']);
catch
    msg(lb,['Failed to evaluate "' str '".'],'r');
    beep;
    return;
end;
if ~isnumeric(r),
    msg(lb,['Expression non-numeric "' str '".'],'r');
    beep;
else
    set(lb,'UserData',r);
    msg(lb,'');
    update(lb);
end;
%=======================================================================

%=======================================================================
function select_all(ob,varargin)
lb = sib(ob,'files');
set(lb,'Value',1:numel(get(lb,'String')));
drawnow;
click_file_box(lb);
return;
%=======================================================================

%=======================================================================
function click_file_box(lb,varargin)
c = get_current_char(lb);
if isempty(c) || isequal(c, char(13))
    vlo  = get(lb,'Value');
    if isempty(vlo),
        msg(lb,'Nothing selected');
        return;
    end;
    lim  = get(lb,'UserData');
    ob   = sib(lb,'selected');
    str3 = get(ob,'String');
    str  = get(lb,'String');
    lim1  = min([max([lim(2)-numel(str3),0]),numel(vlo)]);
    if lim1==0,
        msg(lb,['Selected ' num2str(size(str3,1)) '/' num2str(lim(2)) ' already.']);
        beep;
        set(sib(lb,'D'),'Enable','on');
        return;
    end;

    vl   = vlo(1:lim1);
    msk  = false(size(str,1),1);
    if vl>0, msk(vl) = true; else msk = []; end;
    str1 = str( msk);
    str2 = str(~msk);
    dr   = get(sib(lb,'edit'), 'String');
    str1 = strcat(dr, str1);

    set(lb,'Value',min([vl(1),numel(str2)]),'String',str2);
    r    = (1:numel(str1))+numel(str3);
    str3 = [str3(:);str1(:)];
    set(ob,'String',str3,'Value',r);
    if numel(vlo)>lim1,
        msg(lb,['Retained ' num2str(lim1) '/' num2str(numel(vlo))...
            ' of selection.']);
        beep;
    elseif isfinite(lim(2))
        if lim(1)==lim(2),
            msg(lb,['Selected ' num2str(numel(str3)) '/' num2str(lim(2)) ' files.']);
        else
            msg(lb,['Selected ' num2str(numel(str3)) '/' num2str(lim(1)) '-' num2str(lim(2)) ' files.']);
        end;
    else
        if size(str3,1) == 1, ss = ''; else ss = 's'; end;
        msg(lb,['Selected ' num2str(numel(str3)) ' file' ss '.']);
    end;
    if ~isfinite(lim(1)) || numel(str3)>=lim(1),
        set(sib(lb,'D'),'Enable','on');
    end;
end
return;
%=======================================================================

%=======================================================================
function obj = sib(ob,tag)
persistent fg;
if isempty(fg) || ~ishandle(fg)
    fg = findobj(0,'Tag',mfilename);
end
obj = findobj(fg,'Tag',tag);
return;
%if isempty(obj),
%    cfg_message('matlabbatch:usage',['Can''t find object with tag "' tag '".']);
%elseif length(obj)>1,
%    cfg_message('matlabbatch:usage',['Found ' num2str(length(obj)) ' objects with tag "' tag '".']);
%end;
%return;
%=======================================================================

%=======================================================================
function unselect(lb,varargin)
vl      = get(lb,'Value');
if isempty(vl), return; end;
str     = get(lb,'String');
msk     = true(numel(str),1);
if vl~=0, msk(vl) = false; end;
str2    = str(msk);
set(lb,'Value',min(vl(1),numel(str2)),'String',str2);
lim = get(sib(lb,'files'),'UserData');
if numel(str2)>= lim(1) && numel(str2)<= lim(2),
    set(sib(lb,'D'),'Enable','on');
else
    set(sib(lb,'D'),'Enable','off');
end;

if numel(str2) == 1, ss1 = ''; else ss1 = 's'; end;
%msg(lb,[num2str(size(str2,1)) ' file' ss ' remaining.']);
if numel(vl) == 1, ss = ''; else ss = 's'; end;
msg(lb,['Unselected ' num2str(numel(vl)) ' file' ss '. ' ...
        num2str(numel(str2)) ' file' ss1 ' remaining.']);
return;
%=======================================================================

%=======================================================================
function unselect_all(ob,varargin)
lb = sib(ob,'selected');
set(lb,'Value',[],'String',{},'ListBoxTop',1);
msg(lb,'Unselected all files.');
lim = get(sib(lb,'files'),'UserData');
if lim(1)>0, set(sib(lb,'D'),'Enable','off'); end;
return;
%=======================================================================

%=======================================================================
function [f,d] = listfiles(dr,filt)
try
    ob = sib(gco,'msg');
    domsg = ~isempty(ob);
catch
    domsg = false;
end
if domsg
    omsg = msg(ob,'Listing directory...');
end
if nargin<2, filt = '';  end;
if nargin<1, dr   = '.'; end;
de      = dir(dr);
if ~isempty(de),
    d     = {de([de.isdir]).name};
    if ~any(strcmp(d, '.'))
        d = [{'.'}, d(:)'];
    end;
    if filt.code~=-1,
        f = {de(~[de.isdir]).name};
    else
        % f = d(3:end);
        f = d;
    end;
else
    d = {'.','..'};
    f = {};
end;

if domsg
    msg(ob,['Filtering ' num2str(numel(f)) ' files...']);
end
f  = do_filter(f,filt.ext);
f  = do_filter(f,filt.filt);
ii = cell(1,numel(f));
if filt.code==1 && (numel(filt.frames)~=1 || filt.frames(1)~=1),
    if domsg
        msg(ob,['Reading headers of ' num2str(numel(f)) ' images...']);
    end
    for i=1:numel(f),
        try
            ni = nifti(fullfile(dr,f{i}));
            dm = [ni.dat.dim 1 1 1 1 1];
            d4 = (1:dm(4))';
        catch
            d4 = 1;
        end;
        if all(isfinite(filt.frames))
            msk = false(size(filt.frames));
            for j=1:numel(msk), msk(j) = any(d4==filt.frames(j)); end;
            ii{i} = filt.frames(msk);
        else
            ii{i} = d4;
        end;
    end
elseif filt.code==1 && (numel(filt.frames)==1 && filt.frames(1)==1),
    for i=1:numel(f),
        ii{i} = 1;
    end;
end;

if domsg
    msg(ob,['Listing ' num2str(numel(f)) ' files...']);
end

[f,ind] = sortrows(f(:));
ii      = ii(ind);
msk     = true(1,numel(f));
for i=2:numel(f),
    if strcmp(f{i-1},f{i}),
        if filt.code==1,
            tmp      = sort([ii{i}(:) ; ii{i-1}(:)]);
            tmp(~diff(tmp,1)) = [];
            ii{i}    = tmp;
        end;
        msk(i-1) = false;
    end;
end;
f        = f(msk);
if filt.code==1,
    % Combine filename and frame number(s)
    ii       = ii(msk);
    nii      = cellfun(@numel, ii);
    c        = cell(sum(nii),1);
    fi       = cell(numel(f),1);
    for k = 1:numel(fi)
        fi{k} = k*ones(1,nii(k));
    end
    ii = [ii{:}];
    fi = [fi{:}];
    for i=1:numel(c),
        c{i} = sprintf('%s,%d', f{fi(i)}, ii(i));
    end;
    f        = c;
elseif filt.code==-1,
    fs = filesep;
    for i=1:numel(f),
        f{i} = [f{i} fs];
    end;
end;
f        = f(:);
d        = unique(d(:));
if domsg
    msg(ob,omsg);
end
return;
%=======================================================================

%=======================================================================
function [f,ind] = do_filter(f,filt)
t2 = false(numel(f),1);
filt_or = sprintf('(%s)|',filt{:});
t1 = regexp(f,filt_or(1:end-1));
if numel(f)==1 && ~iscell(t1), t1 = {t1}; end;
for i=1:numel(t1),
    t2(i) = ~isempty(t1{i});
end;
ind = find(t2);
f   = f(t2);
return;
%=======================================================================

%=======================================================================
function heelp(ob,varargin)
[col1,col2,col3,fn] = colours;
fg = sib(ob,mfilename);
t  = uicontrol(fg,...
    'style','listbox',...
    'units','normalized',...
    'Position',[0.01 0.01 0.98 0.98],...
    fn,...
    'BackgroundColor',col2,...
    'ForegroundColor',col3,...
    'Max',0,...
    'Min',0,...
    'tag','HelpWin',...
    'String','                   ');
c0 = uicontextmenu('Parent',fg);
set(t,'uicontextmenu',c0);
uimenu('Label','Done', 'Parent',c0,'Callback',@helpclear);

str  = cfg_justify(t, {[...
'File Selection help. You can return to selecting files via the right mouse button (the "Done" option). '],...
'',[...
'The panel at the bottom shows files that are already selected. ',...
'Clicking a selected file will un-select it. To un-select several, you can ',...
'drag the cursor over the files, and they will be gone on release. ',...
'You can use the right mouse button to un-select everything.'],...
'',[...
'Directories are navigated by editing the name of the current directory (where it says "Dir"), ',...
'by going to one of the previously entered directories ("Prev"), ',...
'by going to one of the parent directories ("Up") or by navigating around ',...
'the parent or subdirectories listed in the left side panel.'],...
'',[...
'Files matching the filter ("Filt") are shown in the panel on the right. ',...
'These can be selected by clicking or dragging.  Use the right mouse button if ',...
'you would like to select all files.  Note that when selected, the files disappear ',...
'from this panel.  They can be made to reappear by re-specifying the directory ',...
'or the filter. '],...
'',[...
'Both directory and file lists can also be browsed by typing the leading ',...
'character(s) of a directory or file name.'],...
'',[...
'Note that the syntax of the filter differs from that used by most other file selectors. ',...
'The filter works using so called ''regular expressions''. Details can be found in the ',...
'MATLAB help text on ''regexp''. ',...
'The following is a list of symbols with special meaning for filtering the filenames:'],...
'    ^     start of string',...
'    $     end of string',...
'    .     any character',...
'    \     quote next character',...
'    *     match zero or more',...
'    +     match one or more',...
'    ?     match zero or one, or match minimally',...
'    {}    match a range of occurrances',...
'    []    set of characters',...
'    [^]   exclude a set of characters',...
'    ()    group subexpression',...
'    \w    match word [a-z_A-Z0-9]',...
'    \W    not a word [^a-z_A-Z0-9]',...
'    \d    match digit [0-9]',...
'    \D    not a digit [^0-9]',...
'    \s    match white space [ \t\r\n\f]',...
'    \S    not a white space [^ \t\r\n\f]',...
'    \<WORD\>    exact word match',...
'',[...
'Individual time frames of image files can also be selected.  The frame filter ',...
'allows specified frames to be shown, which is useful for image files that ',...
'contain multiple time points.  If your images are only single time point, then ',...
'reading all the image headers can be avoided by specifying a frame filter of "1". ',...
'The filter should contain a list of integers indicating the frames to be used. ',...
'This can be generated by e.g. "1:100", or "1:2:100".'],...
'',[...
'The recursive selection button (Rec) allows files matching the regular expression to ',...
'be recursively selected.  If there are many directories to search, then this can take ',...
'a while to run.'],...
'',[...
'There is also an edit button (Ed), which allows you to edit your selection of files. ',...
'When you are done, then use the menu-button of your mouse to either cancel or accept your changes.'],''});
set(t,'String',str);
return;
%=======================================================================

%=======================================================================
function helpclear(ob,varargin)
ob = get(ob,'Parent');
ob = get(ob,'Parent');
ob = findobj(ob,'Tag','HelpWin');
delete(ob);
%=======================================================================

%=======================================================================
function t = cpath(t,d)
switch filesep,
case '/',
    mch = '^/';
    fs  = '/';
    fs1 = '/';
case '\',
    mch = '^.:\\';
    fs  = '\';
    fs1 = '\\';
otherwise;
    cfg_message('matlabbatch:usage','What is this filesystem?');
end

if isempty(regexp(t,mch,'once')),
    if (nargin<2)||isempty(d), d = pwd; end;
    t = [d fs t];
end;

% Replace occurences of '/./' by '/' (problems with e.g. /././././././')
re = [fs1 '\.' fs1];
while ~isempty(regexp(t,re, 'once' )),
    t  = regexprep(t,re,fs);
end;
t  = regexprep(t,[fs1 '\.' '$'], fs);

% Replace occurences of '/abc/../' by '/'
re = [fs1 '[^' fs1 ']+' fs1 '\.\.' fs1];
while ~isempty(regexp(t,re, 'once' )),
    t  = regexprep(t,re,fs,'once');
end;
t  = regexprep(t,[fs1 '[^' fs1 ']+' fs1 '\.\.' '$'],fs,'once');

% Replace '//'
t  = regexprep(t,[fs1 '+'], fs);
%=======================================================================

%=======================================================================
function editwin(ob,varargin)
[col1,col2,col3,lf,bf] = colours;
fg   = gcbf;
lb   = sib(ob,'selected');
str  = get(lb,'String');
ac   = allchild(fg);
acv  = get(ac,'Visible');
h    = uicontrol(fg,...
    'Style','Edit',...
    'units','normalized',...
    'String',str,...
    lf,...
    'Max',2,...
    'Tag','EditWindow',...
    'HorizontalAlignment','Left',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'Position',[0.01 0.08 0.98 0.9],...
    'Userdata',struct('ac',{ac},'acv',{acv}));
ea   = uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',[.01 .01,.32,.07],...
    bf,...
    'Callback',@editdone,...
    'tag','EditWindowAccept',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','Accept');
ee   = uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',[.34 .01,.32,.07],...
    bf,...
    'Callback',@editeval,...
    'tag','EditWindowEval',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','Eval');
ec   = uicontrol(fg,...
    'Style','pushbutton',...
    'units','normalized',...
    'Position',[.67 .01,.32,.07],...
    bf,...
    'Callback',@editclear,...
    'tag','EditWindowCancel',...
    'ForegroundColor',col3,...
    'BackgroundColor',col1,...
    'String','Cancel');
set(ac,'visible','off');
%=======================================================================

%=======================================================================
function editeval(ob,varargin)
ob  = get(ob, 'Parent');
ob  = sib(ob, 'EditWindow');
str = get(ob, 'String');
if ~isempty(str)
    [out,sts] = cfg_eval_valedit(char(str));
    if sts && (iscellstr(out) || ischar(out))
        set(ob, 'String', cellstr(out));
    else
        fgc = get(ob, 'ForegroundColor');
        set(ob, 'ForegroundColor', 'red');
        pause(1);
        set(ob, 'ForegroundColor', fgc);
    end
end
%=======================================================================

%=======================================================================
function editdone(ob,varargin)
ob  = get(ob,'Parent');
ob  = sib(ob,'EditWindow');
str = cellstr(get(ob,'String'));
if isempty(str) || isempty(str{1})
    str = {};
else
    dstr = deblank(str);
    if ~isequal(str, dstr)
        c = questdlg(['Some of the filenames contain trailing blanks. This may ' ...
            'be due to copy/paste of strings between MATLAB and the ' ...
            'edit window. Do you want to remove any trailing blanks?'], ...
            'Trailing Blanks in Filenames', ...
            'Remove', 'Keep', 'Remove');
        switch lower(c)
            case 'remove'
                str = dstr;
        end
    end
    filt = getfilt(ob);
    if filt.code >= 0 % filter files, but not dirs
        [p,n,e] = cellfun(@fileparts, str, 'uniformoutput',false);
        fstr = strcat(n, e);
        [fstr1,fsel] = do_filter(fstr, filt.ext);
        str = str(fsel);
    end
end
lim = get(sib(ob,'files'),'UserData');
if numel(str)>lim(2),
    msg(ob,['Retained ' num2str(lim(2)) ' of the ' num2str(numel(str)) ' files.']);
    beep;
    str = str(1:lim(2));
elseif isfinite(lim(2)),
    if lim(1)==lim(2),
        msg(ob,['Selected ' num2str(numel(str)) '/' num2str(lim(2)) ' files.']);
    else
        msg(ob,['Selected ' num2str(numel(str)) '/' num2str(lim(1)) '-' num2str(lim(2)) ' files.']);
    end;
else
    if numel(str) == 1, ss = ''; else ss = 's'; end;
    msg(ob,['Specified ' num2str(numel(str)) ' file' ss '.']);
end;
if ~isfinite(lim(1)) || numel(str)>=lim(1),
    set(sib(ob,'D'),'Enable','on');
else
    set(sib(ob,'D'),'Enable','off');
end;
set(sib(ob,'selected'),'String',str,'Value',[]);
acs = get(ob,'Userdata');
fg = gcbf;
delete(findobj(fg,'-regexp','Tag','^EditWindow.*'));
set(acs.ac,{'Visible'},acs.acv);
%=======================================================================

%=======================================================================
function editclear(ob,varargin)
fg = gcbf;
acs = get(findobj(fg,'Tag','EditWindow'),'Userdata');
delete(findobj(fg,'-regexp','Tag','^EditWindow.*'));
set(acs.ac,{'Visible'},acs.acv);
%=======================================================================

%=======================================================================
function [c1,c2,c3,lf,bf] = colours
c1 = [1 1 1];
c2 = [1 1 1];
c3 = [0 0 0];
lf = cfg_get_defaults('cfg_ui.lfont');
bf = cfg_get_defaults('cfg_ui.bfont');
if isempty(lf)
    lf = struct('FontName',get(0,'FixedWidthFontName'), ...
                'FontWeight','normal', ...
                'FontAngle','normal', ...
                'FontSize',14, ...
                'FontUnits','points');
end
if isempty(bf)
    bf = struct('FontName',get(0,'FixedWidthFontName'), ...
                'FontWeight','normal', ...
                'FontAngle','normal', ...
                'FontSize',14, ...
                'FontUnits','points');
end
%=======================================================================

%=======================================================================
function select_rec(ob, varargin)
start = get(sib(ob,'edit'),'String');
filt  = get(sib(ob,'regexp'),'Userdata');
filt.filt = {get(sib(ob,'regexp'), 'String')};
fob = sib(ob,'frame');
if ~isempty(fob)
    filt.frames = get(fob,'Userdata');
else
    filt.frames = [];
end;
ptr    = get(gcbf,'Pointer');
try
    set(gcbf,'Pointer','watch');
    sel    = select_rec1(start,filt);
catch
    set(gcbf,'Pointer',ptr);
    sel = {};
end;
set(gcbf,'Pointer',ptr);
already= get(sib(ob,'selected'),'String');
fb     = sib(ob,'files');
lim    = get(fb,'Userdata');
limsel = min(lim(2)-size(already,1),size(sel,1));
set(sib(ob,'selected'),'String',[already(:);sel(1:limsel)],'Value',[]);
msg(ob,sprintf('Added %d/%d matching files to selection.', limsel, size(sel,1)));
if ~isfinite(lim(1)) || size(sel,1)>=lim(1),
    set(sib(ob,'D'),'Enable','on');
else
    set(sib(ob,'D'),'Enable','off');
end;
%=======================================================================

%=======================================================================
function sel=select_rec1(cdir,filt)
sel={};
[t,d] = listfiles(cdir,filt);
if ~isempty(t)
    sel = strcat([cdir,filesep],t);
end;
for k = 1:numel(d)
    if ~any(strcmp(d{k},{'.','..'}))
        sel1 = select_rec1(fullfile(cdir,d{k}),filt);
        sel  = [sel(:); sel1(:)];
    end;
end;
%=======================================================================

%=======================================================================
function sfilt=mk_filter(typ,filt,frames)
if nargin<3, frames  = 1;     end;
if nargin<2, filt    = '.*';    end;
if nargin<1, typ     = 'any';   end;
switch lower(typ),
case {'any','*'}, code = 0; ext = {'.*'};
case {'image'},   code = 1; ext = {'.*\.nii(,\d+){0,2}$','.*\.img(,\d+){0,2}$','.*\.NII(,\d+){0,2}$','.*\.IMG(,\d+){0,2}$'};
case {'mesh'},    code = 0; ext = {'.*\.gii$','.*\.GII$','.*\.mat$','.*\.MAT$'};
case {'nifti'},   code = 0; ext = {'.*\.nii$','.*\.img$','.*\.NII$','.*\.IMG$'};
case {'gifti'},   code = 0; ext = {'.*\.gii$','.*\.GII$'};
case {'extimage'},   code = 1; ext = {'.*\.nii(,[0-9]*){0,2}$',...
                            '.*\.img(,[0-9]*){0,2}$',...
                            '.*\.NII(,[0-9]*){0,2}$',...
                            '.*\.IMG(,[0-9]*){0,2}$'};
case {'xml'},     code = 0; ext = {'.*\.xml$','.*\.XML$'};
case {'mat'},     code = 0; ext = {'.*\.mat$','.*\.MAT$','.*\.txt','.*\.TXT'};
case {'batch'},   code = 0; ext = {'.*\.mat$','.*\.MAT$','.*\.m$','.*\.M$','.*\.xml$','.*\.XML$'};
case {'dir'},     code =-1; ext = {'.*'};
case {'extdir'},     code =-1; ext = {['.*' filesep '$']};
otherwise,        code = 0; ext = {typ};
end;
sfilt = struct('code',code,'frames',frames,'ext',{ext},...
                             'filt',{{filt}});
%=======================================================================

%=======================================================================
function drivestr = listdrives(reread)
persistent mydrivestr;
if isempty(mydrivestr) || reread
    driveLett = strcat(cellstr(char(('C':'Z')')), ':');
    dsel = false(size(driveLett));
    for i=1:numel(driveLett)
        dsel(i) = exist([driveLett{i} '\'],'dir')~=0;
    end
    mydrivestr = driveLett(dsel);
end;
drivestr = mydrivestr;
%=======================================================================

%=======================================================================


text: any(cellfun(@exist,{'add_one_source_c', 'find_existing_source_transpose', 'find_existing_source_transpose_min', 'find_sources_complement_grid_fast_int_all_active_c', 'find_sources_complement_grid_fast_int_c', 'vectorize_quad_single'})==3)

repo: REST


text:
%DEMO_DERIVATIVEOBS  Regression problem demonstration with derivative
%                    observations
%
%  Description
%    The regression problem consist of a data with one input variable,
%    two output variables with Gaussian noise; observations and
%    derivative observations. The constructed model is full GP with
%    Gaussian likelihood.
%
%    The covariance matrix K includes now also covariances between
%    derivative observations and between derivative and latent
%    observations. With derivative observations, the K matrix is a
%    block matrix with following blocks:
%
%        K = [K_ll K_Dl'; K_Dl K_DD]
%
%    Where D refers to derivative and l to latent observation and
%       K_ll = k(x_i, x_j | th)
%       K_Dl = d k(x_i, x_j | th) / dx_i
%       K_DD = d^2 k(x_i, x_j | th) / dx_i dx_j
%
%    To include derivative observations in the inference:
%
%       - provide partial derivative observations in the
%       observation vector after output observations
%       y=[y;dy_1;...;dy_n]; for ex. if size(x)=[10 2] ->
%       size(y)=[30 1]
%
%       - gp_set(gp, 'derivobs', 'on')
%
%   The demo is organised in two parts:
%     1) data analysis without derivative observations
%     2) data analysis with derivative observations
%
%  See also  DEMO_REGRESSION1
%

% Copyright (c) 2010 Tuomas Nikoskinen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Create the data
tp=9;                                  %number of training points -1
x=[-2:4/tp:2]';
y=sin(x).*cos(x).^2;                   % The underlying process
dy=cos(x).^3 - 2*sin(x).^2.*cos(x);    % Derivative of the process
ns=0.06;                               % noise standard deviation

% Add noise
y=y + ns*randn(size(y));
% derivative observations are also noisy
dy=dy + ns*randn(size(dy));
% observation vector with derivative observations
y2=[y;dy];

% test points
xt=[-3:0.05:3]';
nt=length(xt);

%========================================================
% PART 1 GP model without derivative obs
%========================================================
disp('GP model without derivative obs')

% Covariance function
pl = prior_t();
pm = prior_sqrtt();
gpcf = gpcf_sexp('lengthScale', 0.5, 'magnSigma2', .5, ...
                 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
% Use default Gaussian likelihood
gp = gp_set('cf', gpcf);

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3,'DerivativeCheck','on');
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);
% Do the prediction
[Eft, Varft] = gp_pred(gp, x, y, xt);

% PLOT THE DATA

figure
%m=shadedErrorBar(p,Eft(1:size(xt)),2*sqrt(Varft(1:size(xt))),{'k','lineWidth',2});
subplot(2,1,1)
m=plot(xt,Eft,'k','lineWidth',2);
hold on
plot(xt,Eft+2*sqrt(Varft),'k--')
hold on
m95=plot(xt,Eft-2*sqrt(Varft),'k--');
hold on
hav=plot(x, y(1:length(x)), 'ro','markerSize',7,'MarkerFaceColor','r');
hold on
h=plot(xt,sin(xt).*cos(xt).^2,'b--','lineWidth',2);
%legend([m.mainLine m.patch h hav],'prediction','95%','f(x)','observations');
legend([m m95 h hav],'prediction','95%','f(x)','observations');
title('GP without derivative observations')
xlabel('input x')
ylabel('output y')

%========================================================
% PART 2 GP model with derivative obs
%========================================================
disp('GP model with derivative obs')

% Option derivobs set so that the derivatives are in use
gp = gp_set('cf', gpcf, 'derivobs', 'on');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3,'DerivativeCheck','on');
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y2,'opt',opt);
% Do the prediction
[Eft2, Varft2] = gp_pred(gp, x, y2, xt);
% Use predictions for function values only
Eft2=Eft2(1:nt);Varft2=Varft2(1:nt);

% PLOT THE DATA
% plot lines indicating the derivative

subplot(2,1,2)
m=plot(xt,Eft2,'k','lineWidth',2);
hold on
plot(xt,Eft2+2*sqrt(Varft2),'k--')
hold on
m95=plot(xt,Eft2-2*sqrt(Varft2),'k--');
hold on
hav=plot(x, y(1:length(x)), 'ro','markerSize',7,'MarkerFaceColor','r');
hold on
h=plot(xt,sin(xt).*cos(xt).^2,'b--','lineWidth',2);

xlabel('input x')
ylabel('output y')
title('GP with derivative observations')

i1=0;
a=0.1;
ddx=zeros(2*length(x),1);
ddy=zeros(2*length(x),1);
for i=1:length(x)
  i1=i1+1;
  ddx(i1)=x(i)-a;
  ddy(i1)=y(i)-a*dy(i);
  i1=i1+1;
  ddx(i1)=x(i)+a;
  ddy(i1)=y(i)+a*dy(i);
end

for i=1:2:length(ddx)
  hold on
  dhav=plot(ddx(i:i+1), ddy(i:i+1),'r','lineWidth',2);
end
legend([m m95 h hav dhav],'prediction','95%','f(x)','observations','der. obs.');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% MAIN  --  minimum acceleration trajectory
%
% For a simple pendulum:
%
% x = position
% v = velocity
% u = torque
%
% ddx = f(x,dx,u);     <-- dynamics
%
% cost = integral(  ddx^2  );     <-- cost function
%
% subject to:
%   x(0) = 0;
%   x(1) = pi;
%   dx(0) = 0;
%   dx(1) = pi;
%
% How to pose as a standard trajectory optimization problem?
%
% dx = v1;
% dv1 = f(x,v1,u1)
%
% v2 == v1;   % <-- Key line.
% dv2 = u2;
% cost = integral(  u2^2  );
%
%
% NOTES:
%
%   z = [x;v1;v2];
%   u = [u1;u2];
%
clc; clear;
addpath ../../..

%%%% Specify boundary conditions
t0 = 0;
tF = 5;

maxTorque = 1.0;
z0 = [0;0;0];
zF = [pi;0;0];

param.k = 1.0;  % gravity torque constant for pendulum model
param.b = 0.1;  % viscous damping constant

%%%% Pack up boundary conditions
problem.bounds.initialTime.low = t0;
problem.bounds.initialTime.upp = t0;

problem.bounds.finalTime.low = tF;
problem.bounds.finalTime.upp = tF;

problem.bounds.initialState.low = z0;
problem.bounds.initialState.upp = z0;

problem.bounds.finalState.low = zF;
problem.bounds.finalState.upp = zF;

problem.bounds.control.low = [-maxTorque; -inf];
problem.bounds.control.upp = [maxTorque; inf];

%%%% Initialize trajectory with a straight line
problem.guess.time = [t0,tF];
problem.guess.state = [z0, zF];
problem.guess.control = [zeros(2,1), zeros(2,1)];

%%%% Pack up function handles
problem.func.dynamics = @(t,z,u)(  dynamics(z,u,param)  );
problem.func.pathObj = @(t,z,u)(  pathObjective(u)  );
problem.func.pathCst = @(t,z,u)(  pathConstraint(z)  );

%%%% Choice of solver:
method = 'trapezoid';

switch method
    case 'chebyshev'
        problem.options.method = method;
        problem.options.chebyshev.nColPts = 25;
    case 'trapezoid'
        problem.options.method = method;
    case 'hermiteSimpson'
        problem.options.method = method;
        problem.options.hermiteSimpson.nSegment = 15;
        problem.options.nlpOpt.MaxFunEvals = 5e4;
    case 'gpops'
        problem.options.method = 'gpops';
    otherwise
        error('invalid method')
end

%%%% Solve
soln = optimTraj(problem);

%%%% Unpack the solution

tGrid = soln.grid.time;
xGrid = soln.grid.state(1, :);
v1Grid = soln.grid.state(2, :);
v2Grid = soln.grid.state(3, :);
u1Grid = soln.grid.control(1, :);
dv2Grid = soln.grid.control(2, :);

t = linspace(tGrid(1), tGrid(end), 100);
z = soln.interp.state(t);
u = soln.interp.control(t);
x = z(1,:);
v1 = z(2,:);
v2 = z(3,:);
u1 = u(1,:);
dv2 = u(2,:);

%%%% Plot the trajectory against time
figure(1); clf;

subplot(2,2,1); hold on;
plot(t,x)
plot(tGrid,xGrid,'ko','MarkerSize',8,'LineWidth',2);
title('angle')

subplot(2,2,2); hold on;
plot(t,v1)
plot(t,v2)
plot(tGrid,v1Grid,'ko','MarkerSize',8,'LineWidth',2);
plot(tGrid,v2Grid,'ko','MarkerSize',8,'LineWidth',2);
title('angular rate')
legend('v1','v2')

subplot(2,2,3); hold on;
plot(t([1,end]),[1,1]*maxTorque,'k--','LineWidth',1);
plot(t([1,end]),-[1,1]*maxTorque,'k--','LineWidth',1);
plot(t,u1)
plot(tGrid,u1Grid,'ko','MarkerSize',8,'LineWidth',2);
title('torque')

subplot(2,2,4); hold on;
plot(t,dv2)
plot(tGrid,dv2Grid,'ko','MarkerSize',8,'LineWidth',2);
title('angular acceleration')

author: MatthewPeterKelly

repo: OptimTraj


text:
function rgb = CIELab_to_sRGB(Lab)
% Convert a matrix of CIELAB L* a* b* values to sRGB R G B values.
%
% (c) 2018-2020 Stephen Cobeldick
%
%%% Syntax:
% rgb = CIELab_to_sRGB(Lab)
%
% https://en.wikipedia.org/wiki/Lab_color_space
% https://en.wikipedia.org/wiki/SRGB
%
%% Examples %%
%
% >> CIELab_to_sRGB([68.18,2.14,-43.8])*255
% ans =
%   109.0000  169.0000  245.0000
%
% >> CIELab_to_sRGB([0,0,0;100,0,0])
% ans =
%          0         0         0
%     1.0000    1.0000    1.0000
%
%% Inputs and Outputs
%
%%% Input Argument:
% Lab = Numeric Array, size Nx3 or RxCx3, where the last dimension
%       encodes CIELAB values [L*,a*,b*] in the range 0<=L*<=100.
%
%%% Output Argument:
% rgb = Numeric Array, same size as <Lab>, where the last dimension
%       encodes sRGB values [R,G,B] in the range 0<=RGB<=1.
%
% See also CIELAB_TO_DIN99 SRGB_TO_CIELAB SRGB_TO_CAM02UCS SRGB_TO_OSAUCS
% MAXDISTCOLOR MAXDISTCOLOR_VIEW MAXDISTCOLOR_DEMO

%% Input Wrangling %%
%
isz = size(Lab);
assert(isnumeric(Lab),...
	'SC:CIELab_to_sRGB:Lab:NotNumeric',...
	'1st input <Lab> must be numeric.')
assert(isreal(Lab),...
	'SC:CIELab_to_sRGB:Lab:ComplexValue',...
	'1st input <Lab> cannot be complex.')
assert(isz(end)==3,...
	'SC:CIELab_to_sRGB:Lab:InvalidSize',...
	'1st input <Lab> last dimension must have size 3 (e.g. Nx3 or RxCx3).')
Lab = reshape(Lab,[],3);
assert(all(Lab(:,1)>=0&Lab(:,1)<=100),...
	'SC:CIELab_to_sRGB:Lab:OutOfRange',...
	'1st input <Lab> L values must be within the range 0<=L<=100')
%
if ~isfloat(Lab)
	Lab = double(Lab);
end
%
%% Lab2RGB %%
%
M = [... High-precision sRGB to XYZ matrix:
	0.4124564,0.3575761,0.1804375;...
	0.2126729,0.7151522,0.0721750;...
	0.0193339,0.1191920,0.9503041];
% Source: http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
%
wpt = [0.95047,1,1.08883]; % D65
%
% Approximately equivalent to this function, requires Image Toolbox:
%rgb = applycform(lab,makecform('lab2srgb','AdaptedWhitePoint',wpt))
%
% Lab2XYZ
tmp = bsxfun(@rdivide,Lab(:,[2,1,3]),[500,Inf,-200]);
tmp = bsxfun(@plus,tmp,(Lab(:,1)+16)/116);
idx = tmp>(6/29);
tmp = idx.*(tmp.^3) + ~idx.*(3*(6/29)^2*(tmp-4/29));
XYZ = bsxfun(@times,tmp,wpt);
%
% XYZ2RGB
rgb = max(0,min(1,sGammaCor(XYZ / M.')));
%
rgb = reshape(rgb,isz);
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Lab_to_sRGB
function rgb = sGammaCor(rgb)
% Gamma correction of sRGB data.
idx = rgb <= 0.0031308;
rgb(idx) = 12.92 * rgb(idx);
rgb(~idx) = real(1.055 * rgb(~idx).^(1/2.4) - 0.055);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sGammaCor

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function a = dtiGetAlgoHeaderSizeWordAligned()
% Assume a word is 4 bytes

a = dtiGetAlgoHeaderSize();
a = ceil(a/4)*4;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% This simple example shows taking gaussian derivatives on an image
% using the function gfilter.

% set scale, play with this :)
sigma = 2;

% load image and normalize
im  = imread('peppers.png');
im = double( rgb2gray( im ) )/255;

% create figure
figure(1); clf;

% now just calculate and put in subplots...

tmp = gfilter(im,sigma);
subplot(331);
imshow(tmp);
xlabel('Only smoothed');

tmp = gfilter(im,sigma,[0 1]);
subplot(332);
imshow(tmp,[]);
xlabel('Lx');

tmp = gfilter(im,sigma,[0 2]);
subplot(333);
imshow(tmp,[]);
xlabel('Lxx');

tmp = gfilter(im,sigma,[1 0]);
subplot(334);
imshow(tmp,[]);
xlabel('Ly');

tmp = gfilter(im,sigma,[1 1]);
subplot(335);
imshow(tmp,[]);
xlabel('Lxy');

tmp = gfilter(im,sigma,[1 2]);
subplot(336);
imshow(tmp,[]);
xlabel('Lxxy');

tmp = gfilter(im,sigma,[2 0]);
subplot(337);
imshow(tmp,[]);
xlabel('Lyy');

tmp = gfilter(im,sigma,[2 1]);
subplot(338);
imshow(tmp,[]);
xlabel('Lxyy');

tmp = gfilter(im,sigma,[2 2]);
subplot(339);
imshow(tmp,[]);
xlabel('Lxxyy');

author: Sable

repo: mcbench-benchmarks


text:
function pick = nmsMe(boxes, overlap)
% Non-maximum suppression. (FAST VERSION)
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected
% detection.
% NOTE: This is adapted from Pedro Felzenszwalb's version (nms.m),
% but an inner loop has been eliminated to significantly speed it
% up in the case of a large number of boxes

% modified based on Tomasz Malisiewicz's esvm code

if isempty(boxes)
  pick = [];
  return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);
s = boxes(:,end);

area = (x2-x1+1) .* (y2-y1+1);
[~, I] = sort(s);

pick = s*0;
counter = 1;
while ~isempty(I)

  last = length(I);
  i = I(last);
  pick(counter) = i;
  counter = counter + 1;

  xx1 = max(x1(i), x1(I(1:last-1)));
  yy1 = max(y1(i), y1(I(1:last-1)));
  xx2 = min(x2(i), x2(I(1:last-1)));
  yy2 = min(y2(i), y2(I(1:last-1)));

  w = max(0.0, xx2-xx1+1);
  h = max(0.0, yy2-yy1+1);

  o = w.*h ./ area(I(1:last-1));

  I([last; find(o>overlap)]) = [];
end

pick = pick(1:(counter-1));

author: jianxiongxiao

repo: ProfXkit


text:
classdef InertialControlFrame < AbstractControlFrame
    %InertialControlFrame

    properties(Constant)
        enum = ControlFramesEnum.InertialFrame;
    end

    methods
        function obj = InertialControlFrame()

        end

        function enum = getControlFrameEnum(obj)
            enum = ControlFramesEnum.InertialFrame;
        end

        function dcm = computeDcmToInertialFrame(obj, ut, rVect, vVect, bodyInfo, gammaAng, betaAng, alphaAng, baseFrame)
            arguments
                obj(1,1) InertialControlFrame
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                gammaAng(1,1) double
                betaAng(1,1) double
                alphaAng(1,1) double
                baseFrame(1,1) AbstractReferenceFrame
            end

            frame = bodyInfo.getBodyCenteredInertialFrame();
            ce = CartesianElementSet(ut, rVect, vVect, frame);
            ce = ce.convertToFrame(baseFrame);

%             [~, ~, ~, base_frame_2_inertial] = baseFrame.getOffsetsWrtInertialOrigin(ut, ce, baseFrame.getOriginBody());
            base_frame_2_inertial = baseFrame.getRotMatToInertialAtTime(ut, ce, baseFrame.getOriginBody());

            body_2_base_frame = eul2rotmARH([alphaAng,betaAng,gammaAng],'zyx');
            dcm = base_frame_2_inertial * body_2_base_frame;
        end

        function [gammaAngle, betaAngle, alphaAngle] = getAnglesFromInertialBodyAxes(obj, attState, ut, rVect, vVect, bodyInfo, inFrame)
            arguments
                obj(1,1) InertialControlFrame
                attState(1,1) LaunchVehicleAttitudeState
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                inFrame(1,1) AbstractReferenceFrame
            end

            [gammaAngle, betaAngle, alphaAngle] = attState.getEulerAnglesInBaseFrame(ut, rVect, vVect, bodyInfo, inFrame);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function pass = test_nonlinSys1Breaks_US(pref)
% Test 2x2 system (sin/cos). This is pecewiseificaion of the test
%       test_nonlinearSystem1
%
% Asgeir Birkisson, April 2014.

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

tol = 1e-10;

d = [-pi 0 pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Piecewise (ultraS):
pref.discretization = @ultraS;

A = chebop(@(x,u,v) [u - diff(v,2) + u.^2; diff(u) + sin(v)],d);
A.lbc = @(u,v) u-1;
A.rbc = @(u,v) [v-1/2; diff(v)];

u = mldivide(A, f, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = chebfun(A.rbc(u1,u2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(1) = norm( chebfun(A(x, u1, u2))) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
pass(3) = norm(u1jump) < tol && norm(u2jump) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function model = yalmip2coneprog(interfacedata);

options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
Q       = interfacedata.Q;
lb      = interfacedata.lb;
ub      = interfacedata.ub;

ops = options.coneprog;
switch options.verbose
case 0
    ops.Display = 'off';
case 1
    ops.Display = 'final';
otherwise
    ops.Display = 'iter';
end

Aeq = -F_struc(1:1:K.f,2:end);
beq = F_struc(1:1:K.f,1);
A =-F_struc(K.f+1:K.f+K.l,2:end);
b = F_struc(K.f+1:K.f+K.l,1);

if any(K.q)
    top = K.f + K.l + 1;
    for i = 1:length(K.q)
        base = F_struc(top:top+K.q(i)-1,:);
        top = top + K.q(i);
        socConstraints(i) = secondordercone(base(2:end,2:end),-base(2:end,1),base(1,2:end),-base(1));
    end
else
    socConstraints = [];
end

model.f = c;
model.A = A;
model.b = b;
model.Aeq = Aeq;
model.beq = beq;
model.lb = lb;
model.ub = ub;
model.ops = ops;
model.x0 = interfacedata.x0;
model.socConstraints = socConstraints;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function h = vl_plotframe(frames,varargin)
% VL_PLOTFRAME  Plot a geometric frame
%  VL_PLOTFRAME(FRAME) plots the feature frame FRAME. The frame can be
%  either a 2D point, a circle, an oriented circle, an ellipse, or an
%  oriented ellipse, as follows:
%
%  Point::
%    FRAME has 2 components. FRAME(1:2) are the x,y coordinates of the
%    point.
%
%  Circle::
%    FRAME has 3 components. FRAME(1:2) are the x,y coordinates of the
%    center and FRAME(3) is its radius.
%
%  Oriented circle::
%    FRAME has 4 components. FRAME(1:2) are the x,y coordiantes of the
%    center of the circle, FRAME(3) is the radius, and FRAME(4) is the
%    orientation, expressed as a rotation in radians of the standard
%    oriented frame (see below). Positive rotations appear clockwise
%    since the image coordiante system is left-handed.
%
%  Ellipse::
%    FRAME has 5 components. FRAME(1:2) are the x,y coordiantes of the
%    center and FRAME(3:5) are the elements S11, S12, S22 of a 2x2
%    covariance matrix S (a positive semidefinite matrix) defining the
%    ellipse shape. The ellipse is the set of points {x + T: x' inv(S)
%    x = 1}, where T is the ellipse center.
%
%  Oriented ellipse::
%    FAME has 6 components. FRAME(1:2) are the coordiantes T=[x;y] of
%    the center. FRAME(3:6) is the column-wise stacking of a 2x2
%    matrix A. The oriented ellipse is obtained by applying the affine
%    transformation (A,T) to the standard oriented frame (see below).
%
%  A standard unoriented frame is a circle of unit radius centered at
%  the origin; a standard oriented frame is the same, but marked with
%  a radius pointing towards the positive Y axis (downwards) to
%  represent the frame orientation. All other frames can be obtained
%  as affine transformations of these two. In the case of unoriented
%  frames, this transformation is ambiguous up to a rotation.
%
%  VL_PLOTFRAME(FRAMES), where FRAMES is a D x N matrix, plots N
%  frames, one per column. This is significantly more efficient than
%  looping over frames explicitly.
%
%  H = VL_PLOTFRAME(...) returns the handle H of the graphical object
%  representing the frames.
%
%  VL_PLOTFRAME(FRAMES,...) passes any extra argument to the
%  underlying plotting function. The first optional argument, in
%  particular, can be a line specification string such as the one used
%  by PLOT().
%
%  See also: SIFT,
%  covariant detectors,
%  VL_FRAME2OELL(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% Copyright (C) 2013 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% number of vertices drawn for each frame
np = 40 ;

lineprop = {} ;
if length(varargin) > 0
  lineprop = vl_linespec2prop(varargin{1}) ;
  lineprop = {lineprop{:}, varargin{2:end}} ;
end

% --------------------------------------------------------------------
%                                         Handle various frame classes
% --------------------------------------------------------------------

% if just a vector, make sure it is column
if(min(size(frames))==1)
  frames = frames(:) ;
end

[D,K] = size(frames) ;
zero_dimensional = D==2 ;

% just points?
if zero_dimensional
  h = plot(frames(1,:),frames(2,:),'g.',lineprop{:}) ;
  return ;
end

% reduce all other cases to ellipses/oriented ellipses
frames = vl_frame2oell(frames) ;
do_arrows = (D==4 || D==6) ;

% --------------------------------------------------------------------
%                                                                 Draw
% --------------------------------------------------------------------

K   = size(frames,2) ;
thr = linspace(0,2*pi,np) ;

% allx and ally are nan separated lists of the vertices describing the
% boundary of the frames
allx = nan*ones(1, np*K+(K-1)) ;
ally = nan*ones(1, np*K+(K-1)) ;

if do_arrows
  % allxf and allyf are nan separated lists of the vertices of the
  allxf = nan*ones(1, 3*K) ;
  allyf = nan*ones(1, 3*K) ;
end

% vertices around a unit circle
Xp = [cos(thr) ; sin(thr) ;] ;

for k=1:K
  % frame center
  xc = frames(1,k) ;
  yc = frames(2,k) ;

  % frame matrix
  A = reshape(frames(3:6,k),2,2) ;

  % vertices along the boundary
  X = A * Xp ;
  X(1,:) = X(1,:) + xc ;
  X(2,:) = X(2,:) + yc ;

  % store
  allx((k-1)*(np+1) + (1:np)) = X(1,:) ;
  ally((k-1)*(np+1) + (1:np)) = X(2,:) ;

  if do_arrows
    allxf((k-1)*3 + (1:2)) = xc + [0 A(1,2)] ;
    allyf((k-1)*3 + (1:2)) = yc + [0 A(2,2)] ;
  end
end

if do_arrows
  h = line([allx nan allxf], ...
           [ally nan allyf], ...
           'Color','g','LineWidth',3, ...
           lineprop{:}) ;
else
  h = line(allx, ally, ...
           'Color','g','LineWidth',3, ...
           lineprop{:}) ;
end

author: thomaskuestner


text:
function [gmu, gsigmavar, factors] = gpPosteriorGradMeanCovar(model, X);

% GPPOSTERIORGRADMEANCOVAR Gadient of the mean and variances of the posterior at points given by X.
% FORMAT
% DESC computes the gradient of the mean and covariances of the
% posterior distribution of a Gaussian process with respect to the
% input locations.
% ARG model : the model for which gradients are to be computed.
% ARG X : the input locations where gradients are to be computed.
% RETURN gmu : the gradient of the posterior mean with respect to
% the input locations.
% RETURN gCovar : the gradients of the posterior covariance with
% respect to the input locations. By raw, we mean that the gradient
% has not yet been multiplied by any output scale in that direction
% (as is done for gpPosteriorGradMeanCovar). The gradients are
% stored in a cell array of dimension MODEL.q x MODEL.d.
%
% DESC computes the gradient of the mean and covariances of the
% posterior distribution of a Gaussian process with respect to the
% input locations. Returns a compact representation for the
% covariances which separates the factors associated with the
% different dimensions from the covariance gradients.
% ARG model : the model for which gradients are to be computed.
% ARG X : the input locations where gradients are to be computed.
% RETURN gmu : the gradient of the posterior mean with respect to
% the input locations.
% RETURN grCovar : the 'raw' gradient of the posterior covariance with
% respect to the input locations. By raw, we mean that the gradient
% has not yet been multiplied by any output scale in that direction
% (as is done for gpPosteriorGradMeanCovar). The gradients are
% stored in a cell array of length MODEL.q.
% RETURN factors : the factors for multiplying the 'raw' gradients
% of the covariances by.
%
% SEEALSO : gpCreate, gpPosteriorMeanVar
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2009

% GP

if ~isfield(model, 'alpha')
  model = gpComputeAlpha(model);
end

switch model.approx
 case 'ftc'
  gX = kernGradX(model.kern, X, model.X);
  kX_star = kernCompute(model.kern, X, model.X)';
 case {'dtc', 'dtcvar', 'fitc', 'pitc'}
  gX = kernGradX(model.kern, X, model.X_u);
  kX_star = kernCompute(model.kern, X, model.X_u)';
 otherwise
  error('Unrecognised approximation type');

end
K = kernGradX(model.kern, X);

if ~model.isMissingData
  for i = 1:model.q
    switch model.approx
     case 'ftc'
      KinvgK = model.invK_uu*squeeze(gX(:, i, :));
     case {'dtc', 'dtcvar', 'fitc', 'pitc'}
      KinvgK = (model.invK_uu - (1/model.beta)*model.Ainv)*squeeze(gX(:, i, :));
     otherwise
      error('Unrecognised approximation type');
    end
    kXTKinvgK = kX_star'*KinvgK;
    gCovar{i} = squeeze(K(:, i, :))-kXTKinvgK - diag(diag(kXTKinvgK));
    gmu{i} = squeeze(gX(:, i, :))'*model.alpha.*repmat(model.scale, ...
                                                      size(X, 1), 1);
  end

  % Deal with scaling.
  if nargout < 3
    for i = 1:model.q
      for j = 1:model.d
        gsigmavar{i, j} = gCovar{i}*model.scale(j)*model.scale(j);
      end
    end
  else
    factors = model.scale.*model.scale;
    gsigmavar = gCovar;
  end
else
  error('Not yet implemented for models trained on missing data.');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function atlasView = atlasInit(view,atlasName,nAtlases,newAtlasNum,atlasDim)
%
%   atlasView = atlasInit(view,atlasName,nAtlases,newAtlasNum,atlasDim)
%
% Author:
% Purpose:
%    We should use initHiddenFlat for this routine. But that function insists
% on recomputing the coords, and the paths to the flat mat are hard-coded
% and not machine independent, so it often breaks.

atlasView.name = atlasName;
atlasView.viewType = 'Flat';
atlasView.subdir = view.subdir;
atlasView.map = [];
atlasView.mapName = '';
atlasView.curDataType = newAtlasNum;
atlasView.coords = view.coords;
atlasView.grayCoords = view.grayCoords;
atlasView.leftPath = view.leftPath;
atlasView.rightPath = view.rightPath;
atlasView.co = cell(1, nAtlases);
atlasView.amp = cell(1, nAtlases);
atlasView.ph = cell(1, nAtlases);
atlasView.ui.imSize = view.ui.imSize;

% Fill up the co amp and ph with NaNs until the atlas is created.
for(ii=1:nAtlases)
    atlasView.co{ii} = repmat(NaN, atlasDim);
    atlasView.amp{ii} = repmat(NaN, atlasDim);
    atlasView.ph{ii} = repmat(NaN, atlasDim);
end

% In this model, there could be an existing atlas for the other hemisphere.
atlasView = loadCorAnal(atlasView);

return;


text:
figure;

%% Cartesian
% subplot(1,2,1);
% hold on;
% truth_100 = truth.*100;
% %h1 = plot(ax(1),measurement(1,k),measurement(2,k),'k*','MarkerSize', 10);
% meas = 100*obs.heval_inv([measurement(1,:);measurement(2,:)]);
% h1 = plot(meas(obs.Mapping(1),:),meas(obs.Mapping(2),:),'r+','MarkerSize', 15);
% h2 = plot(truth_100(1,1:end),truth_100(2,1:end),'k--','LineWidth',2);
% h3 = plot(truth_100(1,1),truth_100(2,1),'ko','MarkerSize', 20, 'MarkerFaceColor','Green');
% h4 = plot(truth_100(1,end),truth_100(2,end),'ko','MarkerSize', 20, 'MarkerFaceColor','Red');
% h5 = plot(0,0,'rd','MarkerSize', 20, 'MarkerFaceColor','Blue');
% legend([h1,h2,h3,h4,h5],'Measurements','Trajectory', 'Start', 'End','Radar','Location','southeast')
% str = sprintf('Vessel trajectory');
% title(str)
% xlabel('X (m)')
% ylabel('Y (m)')
% axis([0,2500,0,1500])
% box on

%% Polar
% subplot(1,2,2);
% truth2meas = obs.heval(truth_100(1:2,:));
% % [a,b]=obs.heval(truth_100(1:2,:));
% %h1 = polarplot(a(x_start),b(y_start),'ko','LineWidth',2,'MarkerSize',15,'MarkerFaceColor','w');hold on;
% polarplot(truth2meas(1,:),truth2meas(2,:),'k-','LineWidth',2,'MarkerSize',10,'MarkerFaceColor','w');hold on;
% thetalim([0 90])
% %h2 = polarplot(a(x_end),b(y_end),'k^','LineWidth',2,'MarkerSize',15,'MarkerFaceColor','w');hold on;
% h3 = polarplot(measurement(1,:),100*measurement(2,:),'+r','MarkerSize',15);
% pax = gca;
% pax.ThetaAxisUnits = 'radians';
% %rlabel('Range (m)');
% %thetalabel('Bearing (rad)');

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [f,ier]=nufft3d3(x,y,z,c,isign,eps,s,t,u)
%NUFFT3D3: Nonuniform FFT in R^3 - Type 3.
%
%  [FK,IER] = NUFFT3D3(NJ,XJ,YJ,ZJ,CJ,IFLAG,NK,SK,TK,UK);
%
%                 1  nj
%     fk(k)    = -- SUM cj(j) exp(+/-i (s(k),t(k),u(k))*(xj(j),yj(j),zj(j)))
%                nj j=1
%
%     If (isign .ge.0) the + sign is used in the exponential.
%     If (isign .lt.0) the - sign is used in the exponential.
%
%  Input parameters:
%
%     nj     number of sources   (integer)
%     xj,yj,zj  location of sources (real *8)
%
%            on interval [-pi,pi].
%
%     cj     strengths of sources (complex *16)
%     isign  determines sign of FFT (see above)
%     eps    precision request  (between 1.0e-15 and 1.0e-1)
%     nk     number of (noninteger) Fourier modes computed
%     sk,tk,uk  k-values (locations) of desired Fourier modes
%
%  Output parameters:
%
%     fk     Fourier transform values (complex *16)
%     ier    error return code
%            ier = 0  => normal execution.
%            ier = 1  => precision eps requested is out of range.
%
%

nj=numel(x);
nk=numel(s);
if numel(y)~=nj, error('y must have the same number of elements as x'); end
if numel(z)~=nj, error('z must have the same number of elements as x'); end
if numel(c)~=nj, error('c must have the same number of elements as x'); end
if numel(t)~=nk, error('t must have the same number of elements as s'); end
if numel(u)~=nk, error('u must have the same number of elements as s'); end

f=zeros(nk,1)+1i*zeros(nk,1);
ier=0;

mex_id_ = 'nufft3d3f90(i int[x], i double[], i double[], i double[], i dcomplex[], i int[x], i double[x], i int[x], i double[], i double[], i double[], io dcomplex[], io int[x])';
[f, ier] = nufft3d(mex_id_, nj, x, y, z, c, isign, eps, nk, s, t, u, f, ier, 1, 1, 1, 1, 1);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [bound_coeff,vertex_coeff] =  eq_diam_coeff(dim,N)
%EQ_DIAM_COEFF Coefficients of diameter bound and vertex diameter of EQ partition
%
%Syntax
% [bound_coeff,vertex_coeff] = eq_diam_coeff(dim,N);
%
%Description
% [BOUND_COEFF,VERTEX_COEFF] = EQ_DIAM_COEFF(dim,N) does the following:
% 1) uses the recursive zonal equal area sphere partitioning algorithm to
%    partition the unit sphere S^dim into N regions,
% 2) finds the maximum of the per-region diameter bound over all the regions
%    of the partition,
% 3) sets BOUND_COEFF to be the diameter bound coefficient, defined as the
%    solution to
%
%    max_diam_bound == BOUND_COEFF N^(-1/dim),
%
% 4) optionally finds the maximum vertex diameter over all the regions of the
%    partition, and
% 5) optionally sets VERTEX_COEFF to be the vertex diameter coefficient,
%    defined as the solution to
%
%    max_vertex_diam == VERTEX_COEFF N^(-1/dim).
%
% The argument dim must be a positive integer.
% The argument N must be a positive integer or an array of positive integers.
% The result BOUND_COEFF and the optional result VERTEX_COEFF will be arrays of
% the same size as N.
%
%Examples
% > bound_coeff=eq_diam_coeff(2,10)
%  bound_coeff =
%      5.2915
%
% > [bound_coeff,vertex_coeff]=eq_diam_coeff(3,1:6)
%  bound_coeff =
%      2.0000    2.5198    2.8845    3.1748    3.4200    3.6342
%
%  vertex_coeff =
%      2.0000    2.5198    2.8845    3.1748    3.4200    3.6342
%
%See also
% EQ_DIAM_BOUND, EQ_VERTEX_DIAM, EQ_REGIONS, EQ_VERTEX_DIAM_COEFF

% Copyright 2004-2005 Paul Leopardi for the University of NSW.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

%
% Check number of arguments
%
error(nargchk(2,2,nargin));
error(nargoutchk(0,2,nargout));

if nargout < 2
    bound_coeff = eq_diam_bound(dim,N) .* N.^(1/dim);
else
    %
    % Flatten N into a row vector.
    %
    shape = size(N);
    n_partitions = prod(shape);
    N = reshape(N,1,n_partitions);

    bound_coeff =  zeros(size(N));
    vertex_coeff = zeros(size(N));
    for partition_n = 1:n_partitions
        n = N(partition_n);
        regions = eq_regions(dim,n);
        scale = n^(1/dim);
        bound_coeff(partition_n) =  max_diam_bound_of_regions(regions)  * scale;
        vertex_coeff(partition_n) = max_vertex_diam_of_regions(regions) * scale;
    end
    %
    % Reshape output to same array size as original N.
    %
    bound_coeff =  reshape(bound_coeff,shape);
    vertex_coeff = reshape(vertex_coeff,shape);
end
%
%end function

repo: gspbox


text:
classdef ShiftingFunctionComputer < handle

    properties (Access = private)
        LHS
        RHS
    end

    properties (Access = private)
        meshDisc
        mesh
        corrector
        interpolator
    end

    methods (Access = public)

        function obj = ShiftingFunctionComputer(cParams)
            obj.init(cParams);
        end

        function sF = compute(obj)
            obj.computeLHS();
            obj.computeRHS();
            uC = obj.solveSystem();
            In = obj.interpolator;
            u  = In*uC;
            u = reshape(u,obj.mesh.nnodeElem,[]);
            s.mesh = obj.mesh;
            s.fValues(1,:,:) = u;
            sF = P1DiscontinuousFunction(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh     = cParams.mesh;
            obj.corrector    = cParams.corrector;
            obj.interpolator = cParams.interpolator;
            obj.meshDisc     = obj.mesh.createDiscontinuousMesh();
        end

        function computeLHS(obj)
            K = obj.computeStiffnessMatrix();
            In = obj.interpolator;
            K = In'*K*In;
            obj.LHS = K;
        end

        function K = computeStiffnessMatrix(obj)
            s.mesh = obj.mesh;
            s.type = 'StiffnessMatrix';
            s.fun  = P1DiscontinuousFunction.create(obj.mesh, 1);
            lhs = LHSintegrator.create(s);
            K = lhs.compute();
        end

        function computeRHS(obj)
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('QUADRATIC');
            cG = obj.corrector.computeGradient(q);

            s.mesh = obj.meshDisc;
            s.type = 'ShapeDerivative';
            s.quadratureOrder = 'QUADRATIC';
            rhs  = RHSintegrator.create(s);
            rhsF = rhs.compute(cG);
            In = obj.interpolator;
            rhsV = In'*rhsF.fValues;
            obj.RHS = rhsV;
        end

        function u = solveSystem(obj)
            a.type = 'DIRECT';
            s = Solver.create(a);
            u = s.solve(obj.LHS,obj.RHS);
            u = u(1:end);
        end

    end

end


text:
% put a vector in selected positions of a matrix.
%
function [output] = F_copyVec2Mat(input_layer, curr_layer)
input = input_layer.a;

[D,T,N] = size(input);
D1 = curr_layer.targetDims(1);
D2 = curr_layer.targetDims(2);
index2copy = curr_layer.index2copy;

if IsInGPU(input(1))
    output = gpuArray.zeros(D1, D2, T, N);
    tmp = gpuArray.zeros(D1,D2);
else
    output = zeros(D1, D2, T, N);
    tmp = zeros(D1,D2);
end

for t=1:T
    for n = 1:N
        tmp(index2copy) = input(:,t,n);
        output(:,:,t,n) = tmp;
    end
end

end

author: singaxiong

repo: SignalGraph


text:
clear all; close all;
m=256; n=256;
a=50;
b=180;
I=a+(b-a)*rand(m,n);
figure;
subplot(121);  imshow(uint8(I));
subplot(122);  imhist(uint8(I));

author: HuangCongQing


text:
%% Function Serial to Parallel I/P msg and No. of Parallel Channels

function p_data = s2p(s_data,N)
    l = length(s_data);
    mode = mod(l,N);
    if mode ~= 0
        z_add = zeros(1,N-mode);
        data = [s_data z_add];
    else
        data = s_data;
    end

    M = length(data)/N;
    p_data = reshape(data,N,M);
end


text:
function y = feval(varargin)
%FEVAL   Evaluate a DISKFUN at one or more points.
%   Y = FEVAL( F, X, Y) evaluates a diskfun F at a point (X,Y) in Cartesian
%   cooridnates, where X and Y are doubles.
%
%   Y = FEVAL( F, THETA, R, 'polar') evaluates a diskfun F in polar
%   coordinates (THETA,R).  Here THETA and R are doubles representing the
%   central angle (in radians) and radius in polar coordinates and must be
%   points in the unit disk.
%
%   Y = FEVAL(F, c), where c is a complex-valued chebfun representing a
%   contour, evaluates F along the contour.
%
% See also SUBSREF.

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

% Figure out if Cartesian or polar coordinates should be used.
iscart = 1;
% Search for user-supplied 'polar' flag in arguments:
isPolar = find(strcmp(varargin, 'polar'));
if ( any( isPolar ) )
    iscart = 0;
end

% Now evaluate
f = varargin{1};
if nargin < 3 % Eval on a contour parametrized as complex chebfun
    c1 = varargin{2};
    y = chebfun(@(t) feval(f, real(c1(t)), imag(c1(t))), c1.domain, 'vectorize' );
else
    c1 = varargin{2};
    c2 = varargin{3};
    if ( isnumeric(c1) && isnumeric(c2) ) % Eval at a point
        tns =0;
        if ( ndims(c1) >= 3 && isequal(size(c1), size(c2)) )
            % x and y are tensors.
            sizec1 = size(c1);
            c1 = c1(:);
            c2 = c2(:);
            tns =1;
        end
        if iscart
            [theta,r] = cart2pol(c1,c2); % Convert to polar
            if ((any(r > 1+1e-8) )) % Check for points off disk
                error('CHEBFUN:DISKFUN:FEVAL:pointsNotOnDisk',...
                    ['The specified points to evaluate the function do not '...
                    'lie sufficiently close to the unit disk.']);
            end
            y = feval@separableApprox(f, theta, r);
        else
            theta = c1;
            r = c2;
            if ( any(r > 1+1e-8) ) % Check for points off disk
                error('CHEBFUN:DISKFUN:FEVAL:pointsNotOnDisk',...
                    ['The specified points to evaluate the function do not '...
                    'lie sufficiently close to the unit disk.']);
            end
            y = feval@separableApprox(f, theta, r);
        end

        if tns==1
            y = reshape(y, sizec1);
        end
    elseif ( strcmp(c1, ':') && strcmp(c2, ':') ) % Return the diskfun
        y = f;
    elseif ( strcmp(c1, ':') && isnumeric(c2) ) % Angular slice
        y = chebfun(@(t) feval(f, bsxfun(@times, c2, cos(t) ), ...
            bsxfun(@times, c2, sin(t) ) , 'cart'), [-pi, pi], 'trig');
    elseif (isnumeric(c1) && strcmp(c2, ':')) % Radial slice
        y = chebfun(@(t) feval(f, bsxfun(@times, cos(c1), t ),...
            bsxfun(@times, sin(c1), t ), 'cart') );
    else
        error('CHEBFUN:DISKFUN:feval:argin',['Unknown input '...
            'feval(%s,%s)',c1,c2]);
    end
end

end


text:
  function [img sino] = fbp_helix_gh(cg, ig, proj, varargin)
%|function [img sino] = fbp_helix_gh(cg, ig, proj, varargin)
%|
%| A single slice rebinning method for cone-beam tomography data
%| collected with a helical source trajectory
%|
%| in
%|	cg			ct_geom()
%|	ig			image_geom()
%|	proj	[ns nt na]	cone-beam projection views (line integrals)
%|
%| option
%|	'short'	1|0		1 for short-scan fan beam (default); 0 for 360
%|	'chat'	1|0		verbosity
%|
%| out
%|	img	[nx ny nz]	reconstructed image
%|	sino	[ns na_rebin nz] fan-beam sinogram
%|
%| See helix_example, for information on how to call
%|
%| Equations used are taken from
%| Noo F, Defrise M, Clackdoyle R and Kudo H; Phys. Med. Biol. 44:561-70, 1999
%| "Single-slice rebinning for helical cone-beam CT"
%| @u doi 10.1088/0031-9155/44/2/019
%|
%| Copyright 2010-07-21, Gregory Handy and Jeff Fessler, University of Michigan

if nargin == 1 && streq(cg, 'test')
	run_mfile_local helix_example
return
end

if nargin < 3, ir_usage, end

warn 'this file is obsolete; use rebin_helix and fbp_helix_stat'

arg.short = true;
arg.chat = false;
arg = vararg_pair(arg, varargin);

[img sino] = fbp_helix_do(proj, ig.mask_or, ...
		ig.z, ig.dz, ig.dx, ig.nx, ig.ny, ig.nz, ...
		cg.na, cg.ns, cg.ds, cg.s, cg.t, cg.nt, cg.dt, ...
		cg.dsd, cg.dso, cg.dod, cg.dfs, cg.offset_s, cg.wt, ...
		cg.orbit, cg.orbit_start, cg.pitch, cg.rmax, ...
		cg.source_zs, arg.short, arg.chat);

end % fbp_helix_gh()

%
% fbp_helix_do()
%
function [img sino] = fbp_helix_do(proj, mask2, ...
		zslice, dz, dx, nx, ny, nz, ...
		na, ns, ds, spoints, tpoints, nt, dt, ...
		dsd, dso, dod, dfs, offset_s, wt, ...
		orbit, orbit_start, pitch, rmax, source_zs, short, chat)

% step 1: z-sampling
na1 = na/orbit*360; % # of views in one turn - should be integer!
if abs(round(na1) - na1) > 1e-4
	keyboard
	fail 'bad na/orbit'
end
betas = mod(deg2rad(orbit_start + 360 * [0:(na1-1)]'/na1), 2*pi);

% check to see if the data is used "fully"
if dz > (dso/dsd * dt) * 1.1
	warn('Full use of CB data is not achieved')
end

myPitch = pitch * nt * dso / dsd * dt;
delta = asin(rmax/dso); % fan angle (one sided)

% dist = d in the paper, and allows for either a short scan,
% or for a full 360 scan for each z-slice
if short == 1
	dist = 0.5 * myPitch * (pi + 2 * delta)/(2*pi);
else
	dist = 0.5 * myPitch;

	max_t = (max(abs(spoints)) + dsd^2) / (dso*dsd) * dist;
	if max_t >= max(abs(tpoints))
		error('CT geometry does not allow for 360 degree rebinning')
	end
end

% step 2: rebin the cone-beam data into fan-beam projections
fan_beam_proj = zeros(ns, na1, nz);

% loop over the different view angles
if chat
	printm('Rebinning step beginning')
end
for ia=0:na-1
	ticker(mfilename, ia, na)
	currentZ = source_zs(ia+1);

	upper_limit = currentZ + dist;
	lower_limit = currentZ - dist;

	% acceptable range of z-slices
	iz_list = find((lower_limit <= zslice) & (zslice <= upper_limit));
	if isempty(iz_list), continue, end

	% loop over the acceptable z-slices for the current view angle
	for iz=iz_list'
		% Calculate the values of t to be used for each s
		deltaZ = zslice(iz) - currentZ;
		tpoints = (((spoints).^2 + dsd^2) / (dso*dsd)) * deltaZ; % [ns]

		t_index = tpoints / dt + wt; % [ns]
		it0 = floor(t_index); % [ns]
		it1 = 1 + it0;
		alpha = t_index - it0; % [ns] for linear interpolation
		it0 = it0 + 1; % matlab indexing
		it1 = it1 + 1; % matlab indexing

		if any(it0 < 1) || any(it0 > nt), fail 'bug', end
		if any(it1 < 1) || any(it1 > nt), fail 'bug', end

		% scaling factor due to different ray lengths
		scale = sqrt(spoints.^2 + dsd^2) ...
			./ sqrt(spoints.^2 + tpoints.^2 + dsd^2); % [ns]

		tmp = sub2ind([ns nt], 1:ns, it0') + ia * ns * nt;
		y0 = proj(tmp');
		tmp = sub2ind([ns nt], 1:ns, it1') + ia * ns * nt;
		y1 = proj(tmp');
		fan_beam_proj(:, mod(ia, na1) + 1, iz) = ...
			scale .* (alpha.*y1 + (1-alpha).*y0);
	end
end

if short
	if chat, printm('Preparing sinogram for FBP2'), end

	fb_orbit_start = zeros(nz,1);
	% beginning orbit for each z-slice
	new_orbit = deg2rad(orbit_start);
	zloc = source_zs(1);
	change_betas = betas(2)-betas(1);
	for iz = 1:nz
		orbitFound = 0;
		while orbitFound == 0
			if zslice(iz) < zloc + dist
				fb_orbit_start(iz) = new_orbit;
				orbitFound = 1;
			else
				new_orbit = new_orbit + change_betas;
				zloc = zloc + (source_zs(2) - source_zs(1));
			end
		end
	end

	% # of angles for each z-slice sinogram
%	newNA = floor(2 * dist / myPitch * na1); % GH
	newNA = ceil(2 * dist / myPitch * na1); % JF
	new_fan_beam_proj = zeros(ns, newNA, nz);

	% delete the empty rows of the sinograms
	for iz = 1:nz
		angle = round((fb_orbit_start(iz)-deg2rad(orbit_start))...
				/ change_betas);
		for ia = 0:newNA-1
			nextAngle = angle + ia;
			new_fan_beam_proj(:, ia+1, iz) = ...
			fan_beam_proj(:, mod(nextAngle, na1)+1, iz);
		end
	end

	% same image geometry as before, except for the z direction
	ig = image_geom('nx', nx, 'ny', ny, 'dx', dx, 'mask', mask2);

	img = nans(nx, ny, nz);

	if chat, printm('Performing the filter back projection.'), end
	% filter back project for each z-slice
	for iz = 1:nz
		sg = sino_geom('fan', 'ns', ns, 'na', newNA, ...
			'ds', ds, 'offset_s', offset_s, ...
			'dsd', dsd, 'dod', dod, 'dfs', dfs, ...
			'orbit', 'short', ...
			'orbit_start', rad2deg(fb_orbit_start(iz)));

		if iz == 1 % trick: Parker weights do not depend on orbit_start
			[parker_wt scale180] = fbp_fan_short_wt(sg);
		end

		% apply Parker weighting and scaling
		tmp = scale180 * parker_wt .* new_fan_beam_proj(:,:,iz);
		geomsino = fbp2(sg, ig);
		img(:,:,iz) = fbp2(tmp, geomsino);
	end

	if nargout >= 2
		sino = new_fan_beam_proj;
	end

else % 360
	ig = image_geom('nx', nx, 'ny', ny, 'dx', dx);
	mask2 = true([nx ny]);
	mask2(end) = 0; % trick: test it
	ig.mask = repmat(mask2, [1 1]);
	clear mask2

	img = zeros(nx,ny,nz);

	if chat, printm('Performing the filter back projection.'), end
	% filter back project for each z-slice
	for iz = 1:nz
		sg = sino_geom('fan', 'ns', ns, 'na', na1, ...
				'ds', ds, 'orbit', 360, ...
				'orbit_start', orbit_start, ...
				'offset_s', offset_s, ...
				'dsd', dsd, 'dod', dod, 'dfs', dfs);
		geomsino = fbp2(sg, ig);
		img(:,:,iz) = fbp2(fan_beam_proj(:,:,iz), geomsino);
	end

	if nargout >= 2
		sino = fan_beam_proj;
	end
end

end % fbp_helix_do()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function enviwrite(filename,envi_data,datatype,resolu,UL,interleave,zone)
%  write envi format file of uint8 and int16
%  for example: enviwrite('testzz',data,'int16',[30,30],jiUL,'bsq',ZC,'this is for test only');
%  i,j => x, y => rows and cols
%  Input  1)  'filename' is the name of the file
%  Input  2)  'envi_data' is the data matrix
%  Input  3)  'dataype' is the datatype to write 'uint8' or 'int16'
%  Input  4)  'resolu' is the resolution of the pixel
%  Input  5)  'UL' is the UpperLeftPoint X Y of the UL pixel (not center of UL pixel)
%  Input  6)  'interleave' is the bsq bip bil format
%  Input  7)  'zone' is the UTM zone
% History: Projection of Landsat ARD is fixed. 18 Aug., 2020. Shi

if strcmp(datatype,'uint8')
    envi_data=uint8(envi_data);
    dt=1;
elseif strcmp(datatype,'int16')
    envi_data=int16(envi_data);
    dt=2;
elseif strcmp(datatype,'int32')
    envi_data=int32(envi_data);
    dt=3;
elseif strcmp(datatype,'uint16')
    envi_data=uint16(envi_data);
    dt=12;
elseif strcmp(datatype,'single')
    envi_data=single(envi_data);
    dt=4;
else
    fprintf('Invalid write data type!\n');
    return;
end

n_dims=size(envi_data);

nrows=n_dims(1);
ncols=n_dims(2);
bands=1;

if length(n_dims)>=3
    bands=n_dims(3);
end

multibandwrite(envi_data,filename,interleave);

filename_hdr=[filename,'.hdr'];

fid_out=fopen(filename_hdr,'wt');

fprintf(fid_out,'ENVI\n');
fprintf(fid_out,'description = {Landsat Scientific Data}\n');
fprintf(fid_out,'samples = %d\n',ncols); % samples is for j
fprintf(fid_out,'lines   = %d\n',nrows); % lines is for i
fprintf(fid_out,'bands   = %d\n',bands);
fprintf(fid_out,'header offset = 0\n');
fprintf(fid_out,'file type = ENVI Standard\n');
fprintf(fid_out,'data type = %d\n',dt);
fprintf(fid_out,'interleave = %s\n',interleave);
fprintf(fid_out,'sensor type = Landsat\n');
fprintf(fid_out,'byte order = 0\n');

if isempty(zone)
    %label as Landsat ARD
    fprintf(fid_out, 'map info = {Albers Conical Equal Area, 1, 1, %d, %d, %d, %d, WGS-84}\n',UL(1),UL(2),resolu(1),resolu(2));
    fprintf(fid_out, 'projection info = {9, 6378140, 6356755.288157528, 23, -96, 0, 0, 29.5, 45.5,WGS-84, Albers Conical Equal Area}\n');
    fprintf(fid_out, 'coordinate system string = {PROJCS["WGS_1984_Albers",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378140.0,298.257]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-96.0],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_origin",23.0],UNIT["Meter",1.0]]}');
else
    if (zone > 0)
        fprintf(fid_out, 'map info = {UTM, 1.000, 1.000, %d, %d, %d, %d, %d, North, WGS-84, units=Meters}',UL(1),UL(2),resolu(1),resolu(2),zone);
    else
        fprintf(fid_out, 'map info = {UTM, 1.000, 1.000, %d, %d, %d, %d, %d, South, WGS-84, units=Meters}',UL(1),UL(2),resolu(1),resolu(2),-zone);
    end
end

fclose(fid_out);
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC