text:
function succ = install_matpower(modify, save_it, verbose, rm_oldpaths)
%INSTALL_MATPOWER  Assist user in setting path to install MATPOWER.
%   INSTALL_MATPOWER
%   INSTALL_MATPOWER(MODIFY)
%   INSTALL_MATPOWER(MODIFY, SAVE_IT)
%   INSTALL_MATPOWER(MODIFY, SAVE_IT, VERBOSE)
%   INSTALL_MATPOWER(MODIFY, SAVE_IT, VERBOSE, RM_OLDPATHS)
%   SUCCESS = INSTALL_MATPOWER(...)
%
%   Assists the user in setting up the proper MATLAB/Octave path to
%   be able to use MATPOWER and run its tests. With no input arguments
%   it prompts interactively to determine how to handle the paths.
%
%   Inputs (all are optional):
%       MODIFY : select how to set path
%           0 (default) - generate relevant ADDPATH() commands, but
%               don't execute them
%           1 - modify the path directly executing the relevant
%               ADDPATH() commands
%       SAVE_IT : indicates whether or not to save the results
%           0 or [] (default) - don't save any results
%           if MODIFY is 0
%               SAVE_IT = <string> : the relevant ADDPATH() commands
%                   are saved to a file whose name is provided in SAVE_IT
%               SAVE_IT = <other true value> : the relevant ADDPATH() commands
%                   are saved to a file named 'startup.m' in the current
%                   directory
%               otherwise : the commands are displayed, but not saved
%           if MODIFY is 1
%               SAVE_IT = <any true value> : the path will be modified
%                   and saved with SAVEPATH
%               otherwise : the path will be modified but not saved
%       VERBOSE : prints the relevant ADDPATH commands if true (default),
%               silent otherwise
%       RM_OLDPATHS : remove existing installation
%           0 (default) - do NOT remove existing MATPOWER from path
%           1 - remove existing MATPOWER paths first
%
%   Outputs (all are optional):
%       SUCCESS : 1 if all commands succeeded, 0 otherwise
%
%   Examples:
%       install_matpower;           %% print the required ADDPATH() commands
%       install_matpower(0, 1);     %% save the commands to startup.m
%       install_matpower(1, 1);     %% modify my path and save
%       install_matpower(1, 0, 0);  %% modify my path temporarily and silently
%       install_matpower(0, 'matpower6');   %% save the commands to matpower6.m
%       install_matpower(0, 0, 1, 1);   %% uninstall MATPOWER from path
%                                   %% (must call SAVEPATH separately
%                                   %%  to make permanent)
%
%   See also ADDPATH, SAVEPATH.

%   MATPOWER
%   Copyright (c) 2017-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% installation data for each component
min_ver.Octave = '4.0.0';
min_ver.MATLAB = '7.5.0';
install = struct( ...
    'name', { ...       %% (not currently used)
        'matpower', ...
        'most', ...
        'mp-opt-model', ...
        'mips', ...
        'mptest', ...
        'maxloadlim', ...
        'misc', ...
        'reduction', ...
        'sdp_pf', ...
        'se', ...
        'simulink_matpower', ...
        'smartmarket', ...
        'state_estimator', ...
        'syngrid' }, ...
    'dirs', { ...       %% directories to add to path
        {{'lib'}, {'lib', 't'}, {'data'}}, ...
        {{'most', 'lib'}, {'most', 'lib', 't'}}, ...
        {{'mp-opt-model', 'lib'}, {'mp-opt-model', 'lib', 't'}}, ...
        {{'mips', 'lib'}, {'mips', 'lib', 't'}}, ...
        {{'mptest', 'lib'}, {'mptest', 'lib', 't'}}, ...
        {{'extras', 'maxloadlim'}, {'extras', 'maxloadlim', 'tests'}, ...
            {'extras', 'maxloadlim', 'examples'}}, ...
        {{'extras', 'misc'}}, ...
        {{'extras', 'reduction'}}, ...
        {{'extras', 'sdp_pf'}}, ...
        {{'extras', 'se'}}, ...
        {{'extras', 'simulink_matpower'}, ...
            {'extras', 'simulink_matpower', 'helper_functions'}, ...
            {'extras', 'simulink_matpower', 'helper_functions', 'getter_setter_functions'}}, ...
        {{'extras', 'smartmarket'}}, ...
        {{'extras', 'state_estimator'}}, ...
        {{'extras', 'syngrid', 'lib'}, {'extras', 'syngrid', 'lib', 't'}} }, ...
    'fcns', { ...       %% functions to search for to find old paths to remove
        {'mpver', 'test_matpower', 'case9'}, ...
        {'mostver', 'test_most'}, ...
        {'mpomver', 'test_mp_opt_model'}, ...
        {'mipsver', 'test_mips'}, ...
        {'t_begin', 'test_mptest'}, ...
        {'maxloadlim', 'test_mll_main', 'example_ieee9.m'}, ...
        {'checklimits'}, ...
        {'MPReduction'}, ...
        {'sdp_pf_ver'}, ...
        {'run_se'}, ...
        {'SimulinkMATPOWERbase.slx', 'ac_testbed_setup', 'get_losses_from_matpower'}, ...
        {'runmarket'}, ...
        {'runse'}, ...
        {'syngrid', 'test_syngrid'} } ...
);
ni = length(install);   %% number of components

%% default arguments
interactive = 0;
if nargin < 4
    rm_oldpaths = 0;
    if nargin < 3
        verbose = 1;
        if nargin < 2
            save_it = 0;
            if nargin < 1
                modify = [];
                interactive = 1;
            end
        end
    end
end

%% get path to new MATPOWER installation root
[root, n, e] = fileparts(which('install_matpower'));

%% MATLAB or Octave
if exist('OCTAVE_VERSION', 'builtin') == 5
    sw = 'Octave';
else
    sw = 'MATLAB';
end

%% check for required version of MATLAB or Octave
vstr = '';
switch sw
    case 'Octave'
        v = ver('octave');
    case 'MATLAB'
        v = ver('matlab');
        if length(v) > 1
            warning('The built-in VER command is behaving strangely, probably as a result of installing a 3rd party toolbox in a directory named ''matlab'' on your path. Check each element of the output of ver(''matlab'') to find the offending toolbox, then move the toolbox to a more appropriately named directory.');
            v = v(1);
        end
end
if ~isempty(v) && isfield(v, 'Version') && ~isempty(v.Version)
    vstr = v.Version;
    if vstr2num(vstr) < vstr2num(min_ver.(sw))
        error('\n\n%s\n  MATPOWER requires %s %s or later.\n      You are using %s %s.\n%s\n\n', repmat('!',1,45), sw, min_ver.(sw), sw, vstr, repmat('!',1,45));
    end
else
    warning('\n\n%s\n  Unable to determine your %s version. This indicates\n  a likely problem with your %s installation.\n%s\n', repmat('!',1,60), sw, sw, repmat('!',1,60));
end

%% get installation options interactively, if necessary
div_line = sprintf('\n-------------------------------------------------------------------\n\n');
if interactive
    fprintf(div_line);
    fprintf('MATPOWER Installation Options:\n\n');
    fprintf('   1. Do NOT modify the %s path\n', sw);
    fprintf('          (just generate the required ADDPATH commands)\n');
    fprintf('   2. DO modify the %s path, but only temporarily\n', sw);
    fprintf('          (you will have to do it again next time you run %s)\n', sw);
    fprintf('   3. DO modify the %s path, and SAVE the updated path\n', sw);
    fprintf('          (so you will not have to do it again next time you run %s)\n\n', sw);
    s = 0;
    while ~isempty(s) && (s < 1 || s > 3 || rem(s, 1))
        s = str2num(input('Please enter your selection [1, 2, 3] (default = 1) : ', 's'));
    end
    if isempty(s)
        s = 1;              %% default to option 1
    end
    switch s
        case 1              %% don't modify
            modify = 0;
        case 2              %% modify, don't save
            modify = 1;
            save_it = 0;
        case 3              %% modify and save
            modify = 1;
            save_it = 1;
        otherwise
    end
end

%% find available new components
available = zeros(ni, 1);
for i = 1:ni
    if exist(fullfile(root, install(i).dirs{1}{:}), 'dir')
        available(i) = 1;
    end
    % fprintf('%d %s\n', available(i), fullfile(root, install(i).dirs{:}));
end

%% find paths for currently installed components
oldpaths = {};
for i = 1:ni
    if available(i)     %% only if we have a new replacement
        for k = 1:length(install(i).fcns)
            [p, n, e] = fileparts(which(install(i).fcns{k}));
            if ~isempty(p) && ~ismember(p, oldpaths)
                oldpaths{end+1} = p;
            end
        end
    end
end

if ~isempty(oldpaths)
    if exist('mpver', 'file')
        v = mpver;
    else
        v = '';
    end
    str = sprintf('It appears you already have MATPOWER %s or some of its\n', v);
    str = sprintf('%scomponents installed in the following directories:\n\n', str);
    rm_path_str = '';
    for k = 1:length(oldpaths)
        rm_path_str = sprintf('%s    %s\n', rm_path_str, oldpaths{k});
    end
    if interactive
        fprintf(div_line);
        fprintf('%s%s\n', str, rm_path_str);
        fprintf('Unless you are sure you know what you are doing, we recommend\n');
        fprintf('removing these directories from your %s path before proceeding.\n\n', sw);
        s = '?';
        while ~isempty(s) && (s(1) ~= 'Y' && s(1) ~= 'N')
            s = upper(input('Would you like me to do it for you? [Y or N] (default = N) : ', 's'));
        end
        if isempty(s)
            s = 'N';                %% default to N
        end
        if s(1) == 'Y'
            rm_oldpaths = 1;
        end
    elseif ~rm_oldpaths
        error('install_matpower: %s%s\nPlease remove the old installation first, or re-run\nINSTALL_MATPOWER in interactive mode (no arguments).', str, rm_path_str);
    end

    %% remove old paths
    if rm_oldpaths
        rmpath(oldpaths{:});
        if verbose
            fprintf(div_line);
            fprintf('The following directories were removed from your %s path:\n\n%s\n', sw, rm_path_str);
            if ~modify
                fprintf('You will need to manually use SAVEPATH to make the changes permanent.\n');
                if interactive
                    s = input('Hit any key to continue ...', 's');
                end
            end
        end
    end
end

%% generate paths to add
newpaths = {};
for i = 1:ni
    if available(i)     %% only available components
        for k = 1:length(install(i).dirs)
            p = fullfile(root, install(i).dirs{k}{:});
            if ~isempty(p) && ~ismember(p, newpaths)
                newpaths{end+1} = p;
            end
        end
    end
end

%% build addpath string
cmd = sprintf('addpath( ...\n');
for k = 1:length(newpaths)
    cmd = sprintf('%s    ''%s'', ...\n', cmd, newpaths{k});
end
cmd = sprintf('%s    ''%s'' );\n', cmd, '-end');

%% print addpath
if verbose
    fprintf(div_line);
    if modify
        fprintf('Your %s path will be updated using the following command:\n\n%s', sw, cmd);
        if interactive
            s = input(sprintf('\nHit any key to continue ...'), 's');
        end
    else
        fprintf('Use the following command to add MATPOWER to your %s path:\n\n%s\n', sw, cmd);
    end
end

%% add the new paths
if modify
    addpath(newpaths{:}, '-end');
    if verbose
        fprintf('Your %s path has been updated.\n', sw);
    end
end

%% prompt about saving the command to a file
if interactive && ~modify
    s = '?';
    while ~isempty(s) && (s(1) ~= 'Y' && s(1) ~= 'N')
        s = upper(input('Would you like to save this command to a file? [Y or N] (default = N) : ', 's'));
    end
    if isempty(s)
        s = 'N';                %% default to N
    end
    if s(1) == 'Y'
        %% prompt for file name
        s = '?';
        while ~isempty(s) && s(1) == '?'
            s = input('Please enter the name of file? (default = ''startup.m'') : ', 's');
        end
        if isempty(s)
            s = 'startup.m';    %% default to 'startup.m'
        end
        save_it = s;
    else
        save_it = 0;
    end
end

%% save path
success = 1;
if save_it
    if modify                   %% modify the path directly
        savepath;                   %% save the updated path
        if verbose
            fprintf(div_line);
            fprintf('Your updated %s path has been saved using SAVEPATH.\n', sw);
        end
    else                        %% don't modify the path directly
        %% set up file name
        if ~ischar(save_it)
            save_it = 'startup.m';
        end
        [p, n, e] = fileparts(save_it);
        if isempty(e)
            e = '.m';
        end
        fn = fullfile(p, [n e]);

        %% write file
        if exist(fn, 'file')        %% file exists, print a warning
            if verbose
                fprintf(div_line);
                fprintf('The file ''%s'' was not written. A file with that name already exists.', fn);
            else
                error('install_matpower: file ''%s'' not written, a file with that name already exists', fn);
            end
        else                        %% create the file
            fid = fopen(fn, 'a');
            if fid == -1
                success = 0;
            else
                fprintf(fid, '%%%s\n\n', upper(n));
                fprintf(fid, '%%%% add MATPOWER paths\n');
                fprintf(fid, '%s', cmd);
                fclose(fid);
            end
            if verbose
                fprintf(div_line);
                fprintf('The file ''%s'' containing the commands to\nadd MATPOWER to your %s path has been created.\n', fn, sw);
            end
        end
    end
end

if verbose
    fprintf(div_line);
    if modify
        fprintf('Now that you have added the required directories to your %s path\n', sw);
        fprintf('MATPOWER is installed and ready to use.\n\n');
    else
        fprintf('Once you have added the required directories to your %s path\n', sw);
        fprintf('MATPOWER will be installed and ready to use.\n\n');
    end
    fprintf('You may want to begin by typing: mpver\n');
    fprintf('to see the list of installed MATPOWER related software versions.\n\n');
    fprintf('Or to run the MATPOWER test suite to ensure everything is\n');
    fprintf('working correctly, type: test_matpower\n\n');

    if interactive && modify
        s = '?';
        while ~isempty(s) && (s(1) ~= 'Y' && s(1) ~= 'N')
            s = upper(input('Would you like to run the MATPOWER tests now? [Y or N] (default = N) : ', 's'));
        end
        if isempty(s)
            s = 'N';                %% default to N
        end
        if s(1) == 'Y'
            test_matpower;
        end
    end
end

if nargout
    succ = success;
end

function num = vstr2num(vstr)
% Converts version string to numerical value suitable for < or > comparisons
% E.g. '3.11.4' -->  3.011004
pat = '\.?(\d+)';
[s,e,tE,m,t] = regexp(vstr, pat);
b = 1;
num = 0;
for k = 1:length(t)
    num = num + b * str2num(t{k}{1});
    b = b / 1000;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [f, h] = plotBrainGrid(brainGridData, ax, brain_figure, black_brain)
% function plotBrainGrid([brainGridData], [ax])
%
% To plot the wire mesh data loaded from brainGridData.npy.

if nargin<1 || isempty(brainGridData)
    mf = mfilename('fullpath');
    brainGridData = readNPY(fullfile(fileparts(mf), 'brainGridData.npy'));
end

bp = double(brainGridData);
bp(sum(bp,2)==0,:) = NaN; % when saved to uint16, NaN's become zeros. There aren't any real vertices at (0,0,0) and it shouldn't look much different if there were

if nargin<2||isempty(ax)
    if nargin<3||isempty(brain_figure)
        brain_figure = figure('Name','Brain View');
    end

    ax = axes('Parent', brain_figure);
end

if nargin<4||~black_brain
    h = plot3(ax, bp(:,1), bp(:,2), bp(:,3), 'Color', [0 0 0 0.3]);
else
    h = plot3(ax, bp(:,1), bp(:,2), bp(:,3), 'Color', [.7 .7 .7 0.3]);
    set(get(ax, 'Parent'),'color','k')
end

set(ax, 'ZDir', 'reverse')
axis(ax, 'equal');
axis(ax, 'vis3d');
axis(ax, 'off');
f = get(ax, 'Parent');

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


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

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

% Last Modified by GUIDE v2.5 07-Sep-2006 18:11:23

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

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

% --- Executes just before cluster_tool_tools is made visible.
function cluster_tool_tools_OpeningFcn(hObject, eventdata, handles, varargin)
movegui(hObject,'center');
set(handles.title,'string',[cluster_tool_getver]);
handles.pwd=pwd;handles.fdir=pwd;
for k=1:size(varargin(:),1)
    if strcmp(varargin{k},'cl'),handles.cl=varargin{k+1};
    elseif strcmp(varargin{k},'fdir'),handles.fdir=varargin{k+1};
    end
end

% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to cluster_tool_tools (see VARARGIN)

% Choose default command line output for cluster_tool_tools
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes cluster_tool_tools wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in quit.
function quit_Callback(hObject, eventdata, handles)

delete(handles.figure1);

% --- Executes on button press in ret.
function ret_Callback(hObject, eventdata, handles)
if isfield(handles,'cl'),cluster_tool('cl',handles.cl,'fdir',handles.fdir);
else cluster_tool('fdir',handles.fdir);end
delete(handles.figure1);

% --- Executes on button press in save.
function save_Callback(hObject, eventdata, handles)
cluster_tool_savefunc

% --- Executes on button press in mask.
function mask_Callback(hObject, eventdata, handles)
cl_mask=cluster_tool_callgui('cluster_tool_masktypeq');
if cl_mask
    handles.pwd=pwd;cd(handles.fdir);
    [FileName,PathName]=uigetfile('*.mat','Select .mat file containing cl structure');
    cd(handles.pwd);
    if ~PathName,return,end
    handles.fdir=PathName;
    load([PathName FileName]);
    try mask=cl;catch,beep,disp('Error: No cl structure in .mat file'),return,end
else
    handles.pwd=pwd;cd(handles.fdir);
    [FileName,PathName]=uigetfile('*.img','Select mask .img');
    cd(handles.pwd);
    if ~PathName,return,end
    handles.fdir=PathName;
    mask=iimg_read_img([PathName FileName],1);
end
if ~isfield(mask,'M')
    if ~isfield(mask,'mat'),beep,disp('Error: no .M or .mat field in source space.'),return,end
    mask.M=mask.mat;
end
if ~isequal(handles.cl(1).voxSize(:),diag(mask(1).M(1:3,1:3)))
    warndlg('Warning: The mask you have selected is not in the same image space as the cl structure. The cl structure will be transferred into the space of the mask. Note that this may lead to undesirable behavior if the mask space has a higher spatial resolution than the space of the cl structure.')
end
if ~cl_mask
    handles.cl=cl_img_intersect(handles.cl,mask);
else
    handles.cl=cl_cl_intersect(mask,handles.cl,'switch');
end
guidata(hObject,handles);

% --- Executes on button press in write.
function write_Callback(hObject, eventdata, handles)
handles.pwd=pwd;cd(handles.fdir);
[FileName,PathName]=uiputfile('*.img','Save .img file');
cd(handles.pwd);
if ~PathName,return,end
handles.fdir=PathName;
use_values=cluster_tool_callgui('cluster_tool_writemasktypeq');
if use_values
    Z=[];
end
XYZ=[];
for k=1:length(handles.cl)
    if size(handles.cl(k).XYZ,1)~=3
        handles.cl(k).XYZ=handles.cl(k).XYZ';
    end
    XYZ(:,end+1:end+handles.cl(k).numVox)=handles.cl(k).XYZ;
    if use_values
        Z(end+1:end+handles.cl(k).numVox)=handles.cl(k).Z;
    end
end
handles.cl(1).nvox=prod(handles.cl(1).dim(1:3));
ind=sub2ind(handles.cl(1).dim(1:3),XYZ(1,:),XYZ(2,:),XYZ(3,:));
dat=zeros(handles.cl(1).nvox,1);
if use_values
    dat(ind)=Z;
else
    dat(ind)=1;
end
if ~isfield(handles.cl,'mat')
    if ~isfield(handles.cl,'M')
        error('No affine matrix in cl structure. Please use the ''Check structure integrity'' tool (cl_check.m).')
    else
        handles.cl(1).mat=handles.cl(1).M;
    end
end
iimg_write_images(dat,handles.cl(1),[PathName FileName]);
% guidata(hObject,handles);

% --- Executes on button press in subdivide.
function subdivide_Callback(hObject, eventdata, handles)
cluster_tool_subdivide('cl',handles.cl,'fdir',handles.fdir);
close(handles.figure1);
% --- Executes on button press in recluster.
function recluster_Callback(hObject, eventdata, handles)

% --- Executes on button press in table.
function table_Callback(hObject, eventdata, handles)
cluster_tool_table('cl',handles.cl,'fdir',handles.fdir);
close(handles.figure1);

% --- Executes on button press in edit.
function edit_Callback(hObject, eventdata, handles)
cluster_tool_edit('cl',handles.cl,'fdir',handles.fdir);
close(handles.figure1);

% --- Executes on button press in check.
function check_Callback(hObject, eventdata, handles)
warndlg('Please watch the matlab command window and answer any questions asked there')
handles.cl=check_cl(handles.cl);
guidata(hObject,handles);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
delete(hObject);

% --- Executes on button press in get_betas.
function get_betas_Callback(hObject, eventdata, handles)
cluster_tool_getbetas('cl',handles.cl,'fdir',handles.fdir);
close(handles.figure1);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [rs, archetypes] = arch(s, na, mode)

%tstoolbox/@signal/arch
%   Syntax:
%     * [rs, archetypes]=arch(s, na, mode='normalized')
%
%   Input arguments:
%     * na - number of generated archetypes
%     * mode - mode can be one of the following : 'normalized' , 'mean',
%       'raw' (optional)
%
%   Archetypal analysis of column orientated data set:
%     * each row of data is one 'observation', e.g. the sample values of
%       all channels in a multichannel measurement at one point in time
%     * in mode 'normalized' each column of data is centered by removing
%       its mean and then normalized by dividing through its standard
%       deviation before the covariance matrix is calculated
%     * in mode 'mean' only the mean of every column of data is removed
%     * in mode 'raw' no preprocessing is applied to data
%
%   Default value for mode is 'normalized'.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,3);

if nargin < 3
	mode = 'normalized';
end

if ndim(s)~=2
	error('arch needs a signal with two dimensions as input')
end

dat = data(s);
n = dlens(s,1);
m = dlens(s,2);

htext = {''};
htext{end+1} = 'Applied Archetypal Analysis (arch)';

[Z,A,B] = arch2(data(s)', na, mode)

a = achse(unit, 1,1);
a = setname(a, 'Archetyp');

rs = signal(core(A'), s);
rs = addhistory(rs, htext);
rs = addcommandlines(rs, 's = arch(s', na, mode);
rs = setaxis(rs, 2, a);
rs = setplothint(rs, 'multigraph');

archetypes = signal(core(Z),s);
archetypes = addhistory(archetypes, htext);
archetypes = addcommandlines(archetypes, '[dummy, s] = arch(s', na, mode);
archetypes = setaxis(archetypes, 2, a);
archetypes = setaxis(archetypes, 1, getaxis(s,2));

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [u, ind, baseDataTYPE] = getScanGroups(vw, baseDT, confirm)
% Group the scans within a dataTYPE into subgroups with identical annotations.
% The number of subgroups equals the number of scans with unique
% annotations.
%
% [u, ind, baseDataTYPE] = getScanGroups(vw, baseDT)
%
% Purpose:
%   Return the list of unique scan annotatons (u) and the corresponding scan
%   numbers (ind). This is useful if you want to do calculations on all
%   scans with the same annotation (e.g., you might want to make averages
%   of these scans, etc.).
%
% Example
%   [u, ind, baseDataTYPE] = getScanGroups;

% var check
mrGlobals;
%----------------------------------------------------------------------
% variable check
if notDefined('vw'),      vw = getCurView;                       end
if notDefined('baseDT'),  baseDT = viewGet(vw, 'curdatatype');   end
if notDefined('confirm'), confirm = false;                       end
%----------------------------------------------------------------------

% set view to base dataTYPE
vw = viewSet(vw, 'currentDataTYPE', baseDT);

% get the dataTYPE name
baseDataTYPE = dtGet(dataTYPES(baseDT), 'name');

% count the scans
nScans = viewGet(vw, 'nScans');

% get the annotation for each scan
annotation = cell(1,nScans);
for scan = 1:nScans;
    annotation{scan} =  dtGet(dataTYPES(baseDT), 'annotation', scan);
end

% get a list of the unique scan annotations
u = unique(annotation);

% count them
nGroups = length(u);

% get the scan numbers that correspond to each unique scan
ind = cell(1,nGroups);
for scan =1:nGroups;
    ind{scan} = find(strcmp(u{scan},annotation));
end

% confirm the groupings
if confirm
    for group = 1:nGroups
        q{group}  = [u{group} sprintf('\n') ...
            sprintf('%s scans ', baseDataTYPE)...
            num2str(ind{group})...
            sprintf(' -> Group %d\n', group)];
    end
    theanswer = questdlg(q, mfilename);
    if ~isequal(theanswer, 'Yes')
        fprintf('[%s]: Aborting....\n', mfilename);
        u = []; ind = []; baseDataTYPE =[];
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function c = tapas_hgf_ar1_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) for AR(1)
% processes and continuous inputs.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% This file refers to CONTINUOUS inputs (Eqs 48ff in Mathys et al., (2011));
% for binary inputs, refer to tapas_hgf_ar1_binary_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% The phis are estimated in 'logit space' because they are confined to the interval from 0 to 1.
% 'Logit-space' is a logistic sigmoid transformation of native space with a variable upper bound
% a>0:
%
% tapas_logit(x) = ln(x/(a-x)); x = a/(1+exp(-tapas_logit(x)))
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_ar1_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.phi        row vector of phis
%         est.p_prc.m          row vector of ms
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%         est.p_prc.al         alpha
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.dau         input prediction error
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_ar1_config', 'tapas_bayes_optimal_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_ar1';

% Number of levels (minimum: 2)
c.n_levels = 2;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% PLACEHOLDER VALUES
% It is often convenient to set some priors to values
% derived from the inputs. This can be achieved by
% using placeholder values. The available placeholders
% are:
%
% 99991   Value of the first input
%         Usually a good choice for mu_0mu(1)
% 99992   Variance of the first 20 inputs
%         Usually a good choice for mu_0sa(1)
% 99993   Log-variance of the first 20 inputs
%         Usually a good choice for logsa_0mu(1)
%         and logalmu
% 99994   Log-variance of the first 20 inputs minus two
%         Usually a good choice for ommu(1)

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first level, this is usually best
% kept fixed to 1 (determines origin on x_i-scale).
c.mu_0mu = [99991, 1];
c.mu_0sa = [99992, 0];

c.logsa_0mu = [99993, log(0.1)];
c.logsa_0sa = [    1,        1];

% Phis
% Format: row vector of length n_levels.
% Phi is estimated in logit-space because it is
% bounded between 0 and 1
% Fix this to zero (leading to a Gaussian random walk) by
% setting logitphimu = -Inf; logitphisa = 0;
c.logitphimu = [tapas_logit(0.1,1), -Inf];
c.logitphisa = [        10^2,    0];

% ms
% Format: row vector of length n_levels.
% This should be fixed for all levels where the omega of
% the next lowest level is not fixed because that offers
% an alternative parametrization of the same model.
c.mmu = [99991, c.mu_0mu(2)];
c.msa = [99992,           0];

% Kappas
% Format: row vector of length n_levels-1.
% This should be fixed (preferably to 1) if the observation model
% does not use mu_i+1 (kappa then determines the scaling of x_i+1).
c.logkamu = [log(1)];
c.logkasa = [     0];

% Omegas
% Format: row vector of length n_levels
c.ommu = [99994,   -6];
c.omsa = [ 10^2, 10^2];

% Alpha
% Format: scalar
% Fix this to zero (no percpeptual uncertainty) by setting
% logalmu = -Inf; logalsa = 0;
c.logalmu = 99993;
c.logalsa = 2^2;

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.logitphimu,...
    c.mmu,...
    c.logkamu,...
    c.ommu,...
    c.logalmu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.logitphisa,...
    c.msa,...
    c.logkasa,...
    c.omsa,...
    c.logalsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 4*c.n_levels+2*(c.n_levels-1)+2;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_ar1;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_ar1_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function [q,Q_la,Q_lb] = lines2Epoint(La,Lb)

% LINES2EPOINT Intersection point of 2 Plucker lines. Result in Euclidean.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

na = La(1:3);
nb = Lb(1:3);
vb = Lb(4:6);

if nargout == 1

    q = cross(na,nb)/dot(na,vb);

else % jac

    % work in homogeneous for easiear Jacobians

    [pe,PE_na,PE_nb] = crossJ(na,nb);
    [ph,PH_na,PH_vb] = dotJ(na,vb);

    p = [pe;ph];
    Z33 = zeros(3);
    Z13 = zeros(1,3);
    P_la = [PE_na Z33;PH_na Z13];
    P_lb = [PE_nb Z33;Z13 PH_vb];

    [q,Q_p] = hm2eu(p);

    Q_la = Q_p*P_la;

    Q_lb = Q_p*P_lb;

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [rotImg3M,rotMask3M] = perturbImageRotation(img3M,mask3M,angl)
% function [rotImg3M,rotMask3M] = perturbImageRotation(img3M,mask3M,angl)
% angl is in degrees.
%
% Rotates img3M and mask3M by angle angl.
%
% APA, 2/25/2019

rotMask3M = zeros(size(mask3M),'like',mask3M);
rotImg3M = zeros(size(img3M),'like',img3M);
sizV = size(mask3M);
iCtr = round(sizV(1)/2);
jCtr = round(sizV(2)/2);
for slc = 1:size(mask3M,3)
    %     fullSiz = size(mask3M(:,:,slc));
    %     [iV,jV] = find3d(mask3M(:,:,slc));
    %     maskM = mask3M(min(iV):max(iV),min(jV):max(jV),slc);
    maskM = mask3M(:,:,slc);
    imgM = img3M(:,:,slc);
    %     di = max(iV) - min(iV);
    %     dj = max(jV) - min(jV);
    %     iCtr = ceil((max(iV) + min(iV))/2);
    %     jCtr = ceil((max(jV) + min(jV))/2);
    %     minRow = -di*pctJitter/100;
    %     maxRow = di*pctJitter/100;
    %     minCol = -dj*pctJitter/100;
    %     maxCol = dj*pctJitter/100;
    %     numRows = minRow + (maxRow-minRow)*rand(1);
    %     numCols = minCol + (maxCol-minCol)*rand(1);
    %     minAng = -pctJitter/100*180;
    %     maxAng = pctJitter/100*180;
    %     angl = minAng + (maxAng-minAng)*rand(1);
    %     scl = (100-pctJitter)/100 + 2*pctJitter/100*rand(1);
    %     maskM = imtranslate(maskM,[numRows,numCols],'nearest','FillValues',0);
    maskM = imrotate(maskM,angl,'nearest');
    imgM = imrotate(imgM,angl,'nearest');
    %     maskM = imresize(maskM, scl, 'nearest');
    newSiz = size(maskM);
    iStart = 1;
    jStart = 1;
    iEnd = 0;
    jEnd = 0;
    iMin = iCtr - ceil(newSiz(1)/2);
    if iMin < 0
        iStart = 1-iMin;
        iMin = 1;
    else
        iMin = iMin + 1;
    end
    iMax = iCtr + floor(newSiz(1)/2);
    if iMax > sizV(1)
        iEnd = iMax - sizV(1);
        iMax = sizV(1);
    end
    jMin = jCtr - ceil(newSiz(2)/2);
    if jMin < 0
        jStart = 1-jMin;
        jMin = 1;
    else
        jMin = jMin + 1;
    end
    jMax = jCtr + floor(newSiz(2)/2);
    if jMax > sizV(2)
        jEnd = jMax - sizV(2);
        jMax = sizV(2);
    end
    %mask3M(:,:,slc) = 0;
    rotMask3M(iMin:iMax,jMin:jMax,slc) = maskM(iStart:end-iEnd,jStart:end-jEnd);
    rotImg3M(iMin:iMax,jMin:jMax,slc) = imgM(iStart:end-iEnd,jStart:end-jEnd);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function p = he_polynomial_value ( m, n, x )

%*****************************************************************************80
%
%% HE_POLYNOMIAL_VALUE evaluates He(i,x).
%
%  Discussion:
%
%    He(i,x) represents the probabilist's Hermite polynomial.
%
%  Differential equation:
%
%    ( exp ( - 0.5 * x^2 ) * He(n,x)' )' + n * exp ( - 0.5 * x^2 ) * He(n,x) = 0
%
%  First terms:
%
%   1
%   X
%   X^2  -  1
%   X^3  -  3 X
%   X^4  -  6 X^2 +   3
%   X^5  - 10 X^3 +  15 X
%   X^6  - 15 X^4 +  45 X^2 -   15
%   X^7  - 21 X^5 + 105 X^3 -  105 X
%   X^8  - 28 X^6 + 210 X^4 -  420 X^2 +  105
%   X^9  - 36 X^7 + 378 X^5 - 1260 X^3 +  945 X
%   X^10 - 45 X^8 + 630 X^6 - 3150 X^4 + 4725 X^2 - 945
%
%  Recursion:
%
%    He(0,X) = 1,
%    He(1,X) = X,
%    He(N,X) = X * He(N-1,X) - (N-1) * He(N-2,X)
%
%  Orthogonality:
%
%    Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * He(M,X) He(N,X) dX
%    = sqrt ( 2 * pi ) * N! * delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the probabilist's Hermite polynomials
%    of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = x(1:m,1) .* p(1:m,j) - ( j - 1 ) * p(1:m,j-1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function event = read_spmeeg_event(filename, varargin)

% read_spmeeg_event() - import evtns from SPM5 and SPM8 meeg datasets
%
% Usage:
%   >> header = read_spmeeg_event(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   event   - FILEIO toolbox event structure
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Vladimir Litvak

if nargin < 1
    help read_spmeeg_event;
    return;
end

header = ft_getopt(varargin, 'header');

if isempty(header)
    header = read_spmeeg_header(filename);
end

D = header.orig;

event = [];

if isfield(D, 'Radc') % SPM5
    for i = 1:D.Nevents
        if isfield(D, 'events') && isfield(D.events, 'code') && length(D.events.code) == D.Nevents
            value = D.events.code(i);
        else
            value = [];
        end

        event = [event struct('type', 'trial', 'sample', (i-1)*header.nSamples + 1,...
            'value', value, 'offset', -header.nSamplesPre, 'duration', header.nSamples)];
    end

    if D.Nevents == 1 && isfield(D, 'events') && isfield(D.events, 'time')
        event = [event struct('type', 'spm5_event', 'sample', num2cell(D.events.time),...
            'value', num2cell(D.events.code), 'offset', -header.nSamplesPre, 'duration', header.nSamples)];
    end

elseif all(isfield(D, {'type', 'Nsamples', 'Fsample', 'timeOnset'})) % SPM8

    for i = 1:numel(D.trials)
        if isfield(D.trials, 'events')
            cevent = D.trials(i).events;
            if ~isempty(cevent)
                for j = 1:numel(cevent)
                    if ~strcmp(cevent(j).type, 'trial')
                        csample = round((cevent(j).time - D.trials(i).onset)*header.Fs + 1);
                        if csample > 0 && csample <= header.nSamples
                            tmp = rmfield(cevent(j), 'time');
                            if strcmp(D.type, 'continuous') && (D.timeOnset ~= 0)
                                tmp.sample  = (cevent(j).time - D.timeOnset)* header.Fs + 1;
                            else
                                tmp.sample = (i-1)*header.nSamples + 1 +(cevent(j).time - D.trials(i).onset)*header.Fs;
                            end
                            tmp.sample = round(tmp.sample);
                            tmp.offset = 0;
                            event = [event tmp];
                        end
                    end
                end
            end
        end
    end

else
    ft_error('Cannot recognize an SPM EEG header format');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% MESH_BOOLEAN Compute boolean csg operations on "solid", consistently oriented
% meshes.
%
% [W,H] = mesh_boolean(V,F,U,G,operation)
% [W,H] = mesh_boolean(V,F,U,G,operation,'ParameterName',paramter_value, ...)
% [W,H] = mesh_boolean({V1,V2,..,Vn},{F1,F2,...,Fn}, ...
%    operation,'ParameterName',paramter_value, ...)
%
% Inputs:
%   V  #V by 3 list of vertex positions of first mesh
%   F  #F by 3 list of triangle indices into V
%   U  #U by 3 list of vertex positions of second mesh
%   G  #G by 3 list of triangle indices into U
%   operation  followed by operation to perform as a string, one of: 'union',
%     'intersect', 'minus', 'xor', or 'resolve'
%     Optional:
%       'BooleanLib' followed by boolean library back-end to use, one of:
%         {'libigl'}  uses CGAL's exact arithmetic kernel and is believed to be
%                     correct.
%         'cork'  is faster but may give incorrect results.
%         'libigl-try-cork-resolve'  libigl boolean extraction but tries to use
%                                    cork's fast resolve, if intersections
%                                    persist, then resolves remaining with
%                                    libigl's resolve. This adds a "layer of
%                                    robustness" on top of cork, but since it's
%                                    not understood _how_ cork is failing, it
%                                    is unknown whether this will lead to
%                                    correct results.
% Outputs:
%   W  #W by 3 list of vertex positions of boolean result mesh
%   H  #H by 3 list of triangle indices into W
%   J  #H list of indices into [FA;FB] of facet birth parents
%
% See also: self_intersect
%

author: alecjacobson

repo: gptoolbox


text:
classdef nnmnist < nntest
  properties (TestParameter)
    networkType = {'dagnn', 'simplenn'}
  end

  methods (TestClassSetup)
    function init(test)
      addpath(fullfile(vl_rootnn, 'examples', 'mnist'));
    end
  end

  methods (Test)
    function valErrorRate(test, networkType)
      clear mex ; % will reset GPU, remove MCN to avoid crashing
                  % MATLAB on exit (BLAS issues?)
      if strcmp(test.dataType, 'double'), return ; end
      switch test.currentDevice
        case 'cpu'
          gpus = [];
        case 'gpu'
          gpus = 1;
      end
      trainOpts = struct('numEpochs', 1, 'continue', false, 'gpus', gpus, ...
        'plotStatistics', false);
      if strcmp(networkType, 'simplenn')
        trainOpts.errorLabels = {'error', 'top5err'} ;
      end
      [~, info] = cnn_mnist('train', trainOpts, 'networkType', networkType);
      test.verifyLessThan(info.train.error, 0.08);
      test.verifyLessThan(info.val.error, 0.025);
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function point_num = sparse_grid_cc_size_old ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_CC_SIZE_OLD sizes a sparse grid of Clenshaw Curtis points.
%
%  Discussion:
%
%    This function has been replaced by a much faster version.
%
%    It is retained for historical interest only.
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      0 <= LEVEL <= LEVEL_MAX.
%
%    This routine works on an abstract set of nested grids.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique points in the grids.
%

%
%  Special case.
%
  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  The outer loop generates LEVELs from 0 to LEVEL_MAX.
%
  point_num = 0;

  for level = 0 : level_max
%
%  The middle loop generates the next partition that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_closed ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  The inner (hidden) loop generates all points corresponding to given grid.
%
      grid_index = multigrid_index0 ( dim_num, order_1d, order_nd );
%
%  Adjust these grid indices to reflect LEVEL_MAX.
%
      grid_index = multigrid_scale_closed ( dim_num, order_nd, level_max, ...
        level_1d, grid_index );
%
%  Determine the first level of appearance of each of the points.
%
      grid_level = abscissa_level_closed_nd ( level_max, dim_num, order_nd, ...
        grid_index );
%
%  Only keep those points which first appear on this level.
%
      for point = 1 : order_nd

        if ( grid_level(point) == level )

          point_num = point_num + 1;

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
%SF_SIMP_P1 Linear Lagrange shape function for simplices (P1).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_SIMP_P1( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming linear P1 Lagrange shape functions on simplices with
%   values defined in the nodes. XI is Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1-3            Number of space dimensions
%       n_vert      scalar: 2-4            Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim+1]             Local coordinates of evaluation point
%       aInvJac     [n,n_sdim+1*n_sdim]    Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SFLAG1

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function cheby_u_test ( )

%*****************************************************************************80
%
%% CHEBY_U_TEST tests the use of the MEX file CHEBY_U.CPP
%
%  Discussion:
%
%    The file cheby_u.cpp is a C++ function which computes the
%    Chebyshev U polynomials.
%
%    This M file "compiles" cheby_u.cpp, and then shows how it can be called.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_U_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the MEX compiler,\n' );
  fprintf ( 1, '  which allows MATLAB to call C++ functions.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.\n' );
  fprintf ( 1, '  The file CHEBY_U.CPP should be there.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compile the file CHEBY_U.CPP.\n' );

  mex cheby_u.cpp

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.  A new file should show up,\n' );
  fprintf ( 1, '  containing the compiled information.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now use CHEBY_U as though it were a MATLAB M-file\n' );
  fprintf ( 1, '  function.\n' );

  n = 12;
  x = 0.2;

  cx_check(1:n+1) = [ ...
     1.0000000000, ...
     0.4000000000, ...
    -0.8400000000, ...
    -0.7360000000, ...
     0.5456000000, ...
     0.9542400000, ...
    -0.1639040000, ...
    -1.0198016000, ...
    -0.2440166400, ...
     0.9221949440, ...
     0.6128946176, ...
    -0.6770370970, ...
    -0.8837094564 ];

  cx = cheby_u ( n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Evaluation point X = %f\n', x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N           U(N,X)          U(N,X)\n' );
  fprintf ( 1, '           (computed)     (tabulated)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n + 1

    fprintf ( 1, '  %2d  %14f  %14f\n', i-1, cx(i), cx_check(i) );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_U_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function x_stdft = stdft(x, N, K, N_fft)

frames      = 1:K:(length(x)-N);
x_stdft     = zeros(length(frames), N_fft);

w           = hanning(N);
x           = x(:);

for i = 1:length(frames)
    ii              = frames(i):(frames(i)+N-1);
	x_stdft(i, :) 	= fft(x(ii).*w, N_fft);
end

author: mpariente

repo: pystoi

save_path: github-repos/MATLAB/mpariente-pystoi


text:
function RT_moco()
% Display and save motion information at real time. It also shows images and the
% progress of scanning, and allows to check motion information for previous
% series/patients.
%
% To make this work, you will need:
%  1. Set up shared folder at the computer running RT_moco.
%     The folder default to ../incoming_DICOM, but better set to your own by
%     setpref('dicm2nii_gui_para', 'incomingDcm', '/mypath/myIncomingDicom');
%     The result log will be saved into incoming_DICOM/RTMM_log/ folder.
%  2. Set up real time image transfer at Siemens console.

% 200207 xiangrui.li at gmail.com first working version inspired by FIRMM

% Create/re-use GUI and start timer
fh = findall(0, 'Type', 'figure', 'Tag', 'RT_moco');
if ~isempty(fh), figure(fh); return; end
res = get(0, 'ScreenSize');
fh = figure('mc'*[256 1]'); clf(fh);
set(fh, 'MenuBar', 'none', 'ToolBar', 'none', 'NumberTitle', 'off', ...
	'DockControls', 'off', 'CloseRequestFcn', @closeFig, 'Color', [1 1 1]*0.94, ...
    'Name', 'Real Time Image Monitor ', 'Tag', 'RT_moco', 'Visible', 'off');
try fh.WindowState = 'maximized'; catch, fh.Position = [1 60 res(3) res(4)-80]; end
hs.fig = fh;
hs.rootDir = getpref('dicm2nii_gui_para', 'incomingDcm', '../incoming_DICOM/');
hs.backupDir = getpref('dicm2nii_gui_para', 'backupDir', '');
fullName = dicm2nii('', 'fullName', 'func_handle');
hs.rootDir = fullName(hs.rootDir);
if isfolder(hs.backupDir), hs.backupDir = fullName(hs.backupDir); end
hs.logDir = [hs.rootDir 'RTMM_log/'];
if ~isfolder(hs.logDir), mkdir(hs.logDir); end % folder to save subj.mat

h = uimenu(fh, 'Label', '&Patient');
hs.menu(1) = uimenu(h, 'Label', 'Load Patient', 'Callback', @loadSubj);
hs.menu(2) = uimenu(h, 'Label', 'Redo Patient', 'Callback', @redoSubj);
hs.menu(3) = uimenu(h, 'Label', 'Close Patient', 'Callback', @closeSubj);

h = uimenu(fh, 'Label', '&Series');
uimenu(h, 'Label', 'View Selected Series in 3D', 'Callback', @view_3D);
uimenu(h, 'Label', 'Overlay Selected Series onto Anatomy', 'Callback', @overlay);
hs.derived = uimenu(h, 'Label', 'Skip DERIVED Series', 'Checked', 'on', ...
    'Callback', @toggleChecked, 'Separator', 'on');
hs.SBRef = uimenu(h, 'Label', 'Skip *_SBRef Series', 'Callback', @toggleChecked, 'Checked', 'on');

h = uimenu(fh, 'Label', '&View');
uimenu(h, 'Label', 'Reset Brightness', 'Callback', @setCLim);
uimenu(h, 'Label', 'Increase Brightness', 'Callback', @setCLim);
uimenu(h, 'Label', 'Decrease Brightness', 'Callback', @setCLim);
hDV = uimenu(h, 'Label', '&DVARS Threshold', 'Separator', 'on');
for i = [0.01 0.02 0.05 0.1 0.12 0.15 0.2 0.4]
    uimenu(hDV, 'Label', num2str(i), 'Callback', @DV_yLim);
end
uimenu(h, 'Label', 'Show FD plot', 'Callback', @toggleFD, 'Separator', 'on');
hFD = uimenu(h, 'Label', '&FD Axis Range');
for i = [0.18 0.3:0.3:1.5 2.4 3 6]
    uimenu(hFD, 'Label', num2str(i), 'Callback', @FD_range)
end

panel = @(pos)uipanel(fh, 'Position', pos, 'BorderType', 'none');
if res(3) < res(4) % Portrait
    pa1 = panel([0 0.62 1 0.38]); % img and label
    pa2 = panel([0 0 1 0.62]); % table and plot
    axPos = [0.05 0 0.65 1]; % img axis
    lbPos = [0.65 0 0.34 1]; % label axis
    subjPos = [1 0.95]; seriesPos = [1 0.1]; msPos = [1 0.7]; ha = 'right';
else
    pa1 = panel([0 0 0.38 1]);
    pa2 = panel([0.38 0 0.62 1]);
    axPos = [0.05 0.31 0.9 0.67];
    lbPos = [0.05 0.01 0.9 0.3];
    subjPos = [0 0.98]; seriesPos = [1 0.1]; msPos = [0 0.6]; ha = 'left';
end

dy = 0.12 * (0:3);
hs.ax = axes(pa2, 'Position', [0.07 0.5 0.86 0.38], ...
    'NextPlot', 'add', 'XLim', [0.5 300.5], 'UserData', dy, ...
    'TickDir', 'out', 'TickLength', 0.002*[1 1], 'ColorOrder', [0 0 1; 1 0 1]);
xlabel(hs.ax, 'Volume Number');
hs.slider = uicontrol(pa2, 'Units', 'normalized', 'Position', [0.05 0.96 0.9 0.03], ...
    'Style', 'slider', 'Value', 1, 'Min', 1, 'Max', 300, 'Callback', @sliderCB, ...
    'BackgroundColor', 0.5*[1 1 1], 'SliderStep', [1 1]./300);

yyaxis left; ylabel(hs.ax, 'DVARS');
set(hs.ax, 'YTick', dy, 'YLim', dy([1 4]));
c3 = [0 0.8 0;  0.8 0.8 0;  0.5 0 0];
for i = 3:-1:1
    rectangle(hs.ax, 'Position', [0.5 dy(i) 2000 dy(i+1)-dy(i)], ...
        'FaceColor', c3(i,:), 'EdgeColor', c3(i,:), 'LineWidth', 0.01);
end
hs.dv = plot(hs.ax, nan, '.:');

yyaxis right; ylabel(hs.ax, 'Framewise Displacement (mm)');
set(hs.ax, 'YTick', 0:0.4:1.2, 'YLim', [0 1.2]);

txt = @(a)text(hs.ax, 'Units', 'normalized', 'Position', a, 'FontSize', 12, ...
    'BackgroundColor', [1 1 1]*0.94, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
hs.pct(1) = txt([0.995 0.32]); hs.pct(2) = txt([0.995 0.65]);

hs.fd = plot(hs.ax, 0, '.:', 'Visible', 'off');
hs.ax.YAxis(2).Visible = 'off';

ax = axes(pa2, 'Position', [0.07 0.88 0.86 0.04], 'Color', fh.Color);
hs.resp = plot(ax, nan, ones(1,3), 'o', 'Color', 'none', 'MarkerSize', 8);
set(hs.resp(1), 'MarkerFaceColor', [0 0 0]);
set(hs.resp(2), 'MarkerFaceColor', [1 0 0]);
set(hs.resp(3), 'MarkerFaceColor', [0 0.8 0]);
title(ax, ' ', 'FontSize', 14, 'interpreter', 'tex');
set(ax, 'XLim', [0.5 300.5], 'YLim', [0.5 1.5], 'Visible', 'off');

vars = {'Description' 'Series' 'Instances' 'Green' ...
    'Yellow' 'MeanFD'};
w2 = [90 100 90 90 100]; w2 = num2cell([res(3)*pa2.Position(3)*0.94-sum(w2)-24 w2]);
hs.table = uitable(pa2, 'Units', 'normalized', 'Position', [0.02 0.01 0.96 0.42], ...
    'FontSize', 14, 'RowName', [], 'CellSelectionCallback', @tableCB, ...
    'ColumnName', strcat('
', vars, '
'), 'ColumnWidth', w2);

ax = axes(pa1, 'Position', axPos, 'YDir', 'reverse', 'Visible', 'off', 'CLim', [0 1]);
hs.img = image(ax, 'CData', ones(2)*0.94, 'CDataMapping', 'scaled');
axis equal; colormap(ax, 'gray');
hs.instnc = text(ax, 'Units', 'normalized', 'Position', [0.99 0.01], 'Color', 'y', ...
    'FontSize', 14, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom');

ax = axes(pa1, 'Position', lbPos, 'Visible', 'off');
hs.subj = text(ax, 'Position', subjPos, 'FontSize', 24, 'FontWeight', 'bold', ...
    'HorizontalAlignment', ha, 'VerticalAlignment', 'top', 'Interpreter', 'none');
hs.series = text(ax, 'Position', seriesPos, 'FontSize', 18, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom', ...
    'BackgroundColor', fh.Color, 'Interpreter', 'none', 'UserData', [1 0]);
hs.MMSS = text(ax, 'Position', msPos, 'FontSize', 18, 'FontWeight', 'bold', ...
    'HorizontalAlignment', ha, 'VerticalAlignment', 'top', 'Interpreter', 'none', ...
    'BackgroundColor', fh.Color, 'String', {'' ''}, 'Color', 'b');

set(fh, 'HandleVisibility', 'callback', 'Visible', 'on'); % fh.Resize = 'off';

% set up serial port
delete(instrfind('Tag', 'RTMM'));
if ispc, port = 'COM1'; else, port = '/dev/ttyUSB0'; end % change this to yours
hs.serial = serial(port, 'BaudRate', 115200, 'Terminator', '', 'Tag', 'RTMM', ...
    'Timeout', 0.3, 'UserData', struct('fig', fh, 'send', false), ...
    'BytesAvailableFcnCount', 1, 'BytesAvailableFcnMode', 'byte', ...
    'BytesAvailableFcn', @serialRead); %#ok
try fopen(hs.serial); catch, end

hs.timer = timer('StartDelay', 5, 'ObjectVisibility', 'off', 'UserData', fh, ...
    'StopFcn', @saveResult, 'TimerFcn', @doSeries, 'ErrorFcn', @errorLog);

hs.countDown = timer('ExecutionMode', 'fixedRate', 'ObjectVisibility', 'off', ...
    'TimerFcn', {@countDown hs}, 'UserData', 0);

guidata(fh, hs); closeSubj(fh);
start(hs.timer);

%% Log error for debugging
function errorLog(obj, evt) %#ok
hs = guidata(obj.UserData);
vnam = "err_"+hs.subj.String+"_"+hs.series.String{1};
eval(vnam+" = evt.Data;");
fnam = hs.logDir+"errorLog.mat";
if isfile(fnam), save(fnam, vnam, '-append'); else, save(fnam, vnam); end

%% TimerFunc: do a series if avail, then call stopFunc to save result
function doSeries(obj, ~)
hs = guidata(obj.UserData);
hs.fig.Name(end) = 78 - hs.fig.Name(end); % dot/space switch: timer indicator
if hs.timer.StartDelay > 1, return; end % no new series
set(hs.menu, 'Enable', 'off'); set([hs.table hs.slider], 'Enable', 'inactive');

dict = dicm_dict('', {'Rows' 'Columns' 'BitsAllocated' 'InstanceNumber'});
try hs.fig.WindowState = 'maximized'; catch, end
L0 = java.awt.MouseInfo.getPointerInfo().getLocation();
java.awt.Robot().mouseMove(L0.getX+9, L0.getY+9); % wake up screen
pause(0.1); java.awt.Robot().mouseMove(L0.getX, L0.getY);

iRun = hs.series.UserData; % updated in new_series()
f = sprintf('%s/%03u_%06u_', hs.subj.UserData, iRun);
s = dicm_hdr_wait([f '000001.dcm']);
if isempty(s), return; end % non-image dicom, skip series
if all(iRun == 1) % first series: reset GUI
    closeSubj(hs.fig);
    hs.subj.String = regexprep(s.PatientName, '[_\s]', '');
    close(findall(0, 'Type', 'figure', 'Tag', 'nii_viewer')); % last subj if any
    if now-getfield(dir(s.Filename), 'datenum')<5/1440
        save([hs.rootDir '/hdr_' hs.subj.String], 's'); % as new subj flag
    end
end

if hs.derived.Checked=="on" && contains(s.ImageType, 'DERIVED'), return; end
if hs.SBRef.Checked=="on" && endsWith(s.SeriesDescription, '_SBRef'), return; end
if startsWith(s.SequenceName, 'ABCD3d1'), return; end

isDTI = contains(s.ImageType, '\DIFFUSION');
isMOS = contains(s.ImageType, '\MOSAIC');
isMoCo = contains(s.ImageType, '\MOCO');
isEnh = isfield(s, 'NumberOfFrames') && s.NumberOfFrames>1;

nTE = asc_header(s, 'lContrasts', 1);
% fieldmap phase diff series: nTE=2, EchoNumber=2
if (contains(s.ImageType, '\P\') && contains(s.SequenceName, 'fm2d')) || ...
   (contains(s.ImageType, '\MEAN') && s.NumberOfAverages>1) % vNAV RMS
    nTE = 1;
end
if isDTI
    nV = [];
    for i = 1:999
        a = asc_header(s, sprintf('sDiffusion.alAverages[%i]', i-1));
        if isempty(a), break; else, nV(i) = a; end %#ok
    end
    nV = nV(1) + sum(nV(2:end))*asc_header(s, 'sDiffusion.lDiffDirections');
else
    nV = nTE * (asc_header(s, 'lRepetitions', 0)+1);
end
if ~isMoCo, hs.series.String = seriesInfo(s); end
if isMOS, nSL = s.CSAImageHeaderInfo.NumberOfImagesInMosaic; end

if isDTI && (isMOS || isEnh)
    nIN = nV;
elseif contains(s.SequenceName, {'epfid2d' 'ASL'}) && (isMOS || isEnh) % EPI/ASL
    nIN = asc_header(s, 'lRepetitions', 0) + 1;
else % T1, T2, fieldmap etc: show info/img only
    if s.MRAcquisitionType == "3D", nIN = asc_header(s, 'sKSpace.lImagesPerSlab');
    else, nIN = asc_header(s, 'sSliceArray.lSize'); % 2D
    end
    if isEnh, nIN = nIN / s.NumberOfFrames; end
    init_series(hs, s, nIN*nV, nV);
    nam = sprintf('%s%06u.dcm', f, nIN*nTE);
    tEnd = now + 9/86400; % wait several seconds for files to arrive
    while now<tEnd && ~isfile(nam), pause(0.2); end
    if ~isfile(nam)
        hs.instnc.String = '1';
        img = dicm_img(s);
    else
        if isEnh
            hs.instnc.String = '1';
            img = permute(dicm_img(s), [1 2 4 3]);
        else
            hs.instnc.String = '';
            nii = dicm2nii([s.Filename(1:end-10) '*'], ' ', 'no_save');
            iSL = bitand(nii.hdr.dim_info, 48)/16;
            img = nii.img(:,:,:,1);
            if     iSL == 1, img = permute(img, [3 2 1]);
            elseif iSL == 2, img = permute(img, [3 1 2]);
            elseif iSL == 3, img = permute(img, [2 1 3]);
            end
            img = flip(img, 1);
        end
        sz = getpixelposition(hs.img.Parent);
        dim = [size(img) 1 1];
        nMos = round(min(sz(3:4)*2./dim(1:2))).^2;
        img(:,1,:) = 0; img(1,:,:) = 0; % visual divider
        img = vol2mos(img(:,:,unique(round(linspace(1, dim(3), nMos)))));
    end
    hs.slider.Value = 1;
    set_img(hs.img, img);
    return;
end

if isempty(nIN) || endsWith(s.SeriesDescription, '_SBRef'), nIN = 1; end
if nIN>5 && hs.countDown.Running == "off" % work for mosaic without tricks
    a = dir(s.Filename);
    hs.countDown.UserData = (nIN-1) * asc_header(s,'alTR[0]')/864e8 + a.datenum;
    start(hs.countDown);
end
try stopAt = str2double(regexp(s.ImageComments, '(?<=stopAt:)\d*', 'match', 'once'));
catch, stopAt = inf;
end

mos = dicm_img(s);
if isMOS, img = mos2vol(mos, nSL);
else, img = permute(mos, [1 2 4 3]); mos = vol2mos(img);
end
if isEnh
    multiFrameFields = dicm2nii('', 'multiFrameFields', 'func_handle');
    s = multiFrameFields(s);
end
img = brainMask(img, s);
try thk = s.SpacingBetweenSlices; catch, thk = s.SliceThickness; end
p = refVol(img, [s.PixelSpacing' thk]);
if isDTI
    p.mean = p.mean*5; % make it similar to EPI
end
img0 = double(img);
set_img(hs.img, mos);
init_series(hs, s, nIN, nIN);
hs.img.UserData(1) = mean(img(:));
viewer = findall(0, 'Type', 'figure', 'Tag', 'nii_viewer');
if nIN<6 && isempty(viewer), overlay(hs.fig); end

R1 = inv(p.R0);
m6 = zeros(2,6);

nextSeries = sprintf('%s/%03u_%06u_000001.dcm', hs.subj.UserData, iRun+[0 1]);
for i = 2:nIN
    nam = sprintf('%s%06u.dcm', f, (i-1)*nTE+1);
    tEnd = now + 1/1440; % 1 minute no mosaic coming, treat as stopped series
    while ~isfile(nam)
        if isfile(nextSeries) || now>tEnd, return; end
        if hs.countDown.Running=="on", countDown(hs.countDown, 0, hs); end
        if hs.serial.Status=="open", serialRead(hs.serial); end
        pause(0.2);
    end
    s = dicm_hdr_wait(nam, dict); iN = s.InstanceNumber;
    if nTE>1, iN = i; end % fake InstanceNumber for multi-echo EPI
    mos = dicm_img(s);
    if isMOS, img = mos2vol(mos, nSL);
    else, img = permute(mos, [1 2 4 3]); mos = vol2mos(img);
    end
    hs.img.UserData(iN) = mean(img(:));
    set_img(hs.img, mos); hs.instnc.String = num2str(iN);
    hs.slider.Value = iN; % show progress
    img = brainMask(img);

    if isMoCo % FD from dicom hdr, DV uses MoCo img for now
        if isEnh
            s1 = dicm_hdr(nam, dicm_dict('Siemens', {'RBMoCoTrans' 'RBMoCoRot'}));
            m6(2,:) = [s1.RBMoCoTrans; s1.RBMoCoRot];
        else
            s1 = dicm_hdr(nam, dicm_dict('Siemens', 'CSAImageHeaderInfo'));
            m6(2,:) = [s1.CSAImageHeaderInfo.RBMoCoTrans; ...
                s1.CSAImageHeaderInfo.RBMoCoRot];
        end
    else
        p.F.Values = smooth_mc(img, p.sz);
        [m6(2,:), R1] = moco_estim(p, R1);
    end
    a = abs(m6(2,:) - m6(1,:)); m6(1,:) = m6(2,:);
    hs.fd.YData(iN) = sum([a(1:3) a(4:6)*50]); % 50mm: head radius

    img = double(img);
    a = img(:) - img0(:);
    hs.dv.YData(iN) = sqrt(a'*a / numel(a)) / p.mean;
    img0 = img;
    if hs.serial.UserData.send
        fwrite(hs.serial, uint8(hs.dv.YData(iN) / hs.ax.YAxis(1).Limits(2)*255));
    end

    a = hs.dv.YData(1:iN); a = [0 a(~isnan(a))];
    dy = hs.ax.UserData; fd = hs.fd.YData(1:iN);
    N = {numel(a) sum(a<dy(2)) sum(a<dy(3))};
    hs.table.Data(1,3:6) = [N mean(fd, 'omitnan')];
    for j = 1:2, hs.pct(j).String = sprintf('%.3g%%', N{j+1}/N{1}*100); end
    if i==2, set(hs.pct, 'Visible', 'on'); end
    if N{2}>=stopAt, [~, ~] = system(['touch "' hs.rootDir 'StopScan"']); end
    if iN>=nIN, return; end % ISSS alike
    % drawnow; % update instance for offline test
end

%% Reshape mosaic into volume, remove padded zeros
function vol = mos2vol(mos, nSL)
nMos = ceil(sqrt(nSL)); % nMos x nMos tiles for Siemens
[nr, nc] = size(mos); % number of row & col in mosaic
nr = nr / nMos; nc = nc / nMos; % number of row and col in slice
vol = zeros([nr nc nSL], class(mos));
for i = 1:nSL
    r = floor((i-1)/nMos) * nr + (1:nr); % 2nd slice is tile(1,2)
    c =    mod(i-1, nMos) * nc + (1:nc);
    vol(:, :, i) = mos(r, c);
end

%% Initialize GUI for a new series
function init_series(hs, s, nIN, nV)
if size(hs.table.Data,1)>0 && hs.table.Data{1,2}==s.SeriesNumber, return; end
fid = fopen([hs.rootDir 'currentSeries.txt'], 'w');
fprintf(fid, '%s_%s_%s', s.PatientName, asc_header(s, 'tProtocolName'), dicmTime(s));
fclose(fid);

set(hs.slider, 'Max', nIN, 'Value', 1, 'UserData', s.Filename(1:end-10));
if nIN==1, hs.slider.Visible = 'off';
else, set(hs.slider, 'SliderStep', [1 1]./(nIN-1)); hs.slider.Visible = 'on';
end
hs.dv.YData = nan(nV,1); hs.fd.YData = nan(nV,1); hs.img.UserData = nan(nV,1);
hs.fig.UserData.nIN{end+1} = nIN;
if nV>1 && nIN==nV
    set([hs.ax hs.pct], 'Visible', 'on');
    hs.ax.XLim(2) = nV + 0.5;
else
    set([hs.ax hs.pct], 'Visible', 'off');
end
hs.resp(1).Parent.XLim(2) = hs.ax.XLim(2);
set(hs.resp, 'XData', nan, 'YData', 1); update_resp(hs);
set([hs.instnc hs.pct], 'String', '');
figure(hs.fig); drawnow; % bring GUI front if needed
if contains(s.ImageType, '\MOCO'), return; end
hs.table.Data = [{s.SeriesDescription s.SeriesNumber nIN [] [] []}; hs.table.Data];
hs.fig.UserData.hdr{end+1} = s; % 1st instance

pat = asc_header(s,'sPat.lAccelFactPE', 1);
if pat==1, thr = 0.12; else, thr = 0.15; end % arbitrary
h = findobj(hs.fig, 'Type', 'uimenu', 'Label', '&DVARS Threshold');
thrs = str2double(get(h.Children, 'Label'));
[~, i] = min(abs(thrs-thr));
DV_yLim(h.Children(i));

%% Set img and img axis
function set_img(h, img)
d = size(img) + 0.5;
set(h.Parent, 'CLim', [0 imgClim(img)], 'XLim', [0.5 d(2)], 'YLim', [0.5 d(1)]);
h.CData = img;

%% get some series information
function c = seriesInfo(s)
c{1} = s.SeriesDescription;
if numel(c{1})>24, c{1} = [c{1}(1:16) '...' c{1}(end-3:end)]; end
c{2} = sprintf('Series %g', s.SeriesNumber);
if s.StudyID~="1", c{2} = ['Study ' s.StudyID ', ' c{2}]; end
dcmT = dicmTime(s);
c{3} = datestr(datenum(dcmT, 'yymmddHHMMSS.fff'), 'HH:MM:SS AM');
c{4} = datestr(datenum(dcmT, 'yymmdd'), 'ddd, mmm dd, yyyy');
try c{5} = sprintf('TR = %g', asc_header(s, 'alTR[0]')/1000); catch, end

%% toggle FD display on/off
function toggleFD(h, ~)
hs = guidata(h);
if hs.fd.Visible == "on"
    set([hs.fd hs.ax.YAxis(2)], 'Visible', 'off');
    h.Label = 'Show FD plot';
else
    set([hs.fd hs.ax.YAxis(2)], 'Visible', 'on');
    h.Label = 'Hide FD plot';
end

%% Set FD plot y-axis limit
function FD_range(h, ~)
hs = guidata(h);
dy = str2double(h.Label) * (0:3) / 3;
yyaxis(hs.ax, 'right'); set(hs.ax, 'YTick', dy, 'YLim', dy([1 4]));

%% Set DVARS plot y-axis limit, and update table
function DV_yLim(h, ~)
hs = guidata(h);
dy = str2double(h.Label) * (0:3);
hs.ax.UserData = dy;
for i = 1:numel(hs.fig.UserData.DV)
    a = hs.fig.UserData.DV{i}; a = [0 a(~isnan(a))];
    if numel(a)<2, continue; end
    hs.table.Data(end-i+1, 4:5) = {sum(a<dy(2)) sum(a<dy(3))};
end
yyaxis(hs.ax, 'left'); set(hs.ax, 'YTick', dy, 'YLim', dy([1 4]));
rect = findobj(hs.ax, 'type', 'Rectangle');
for i = 1:3, rect(i).Position([2 4]) = dy([i 2]); end
DV = hs.dv.YData; a = [0 DV(~isnan(DV))];
for i = 1:2, hs.pct(i).String = sprintf('%.3g%%', sum(a<dy(i+1))/numel(a)*100); end

%% Table-click callback: show moco/series info and image if avail
function tableCB(h, evt)
if isempty(evt.Indices) || evt.Indices(1,2)>2, return; end
hs = guidata(h);
C = h.Data;
iT = evt.Indices(1,1);
iR = size(C,1) - iT + 1;
dV = hs.fig.UserData.DV{iR};
hs.fd.YData = hs.fig.UserData.FD{iR};
hs.dv.YData = dV;
if numel(dV)<2 || isnan(dV(2))
    set([hs.ax hs.pct], 'Visible', 'off');
else
    set([hs.ax hs.pct], 'Visible', 'on');
    for j = 1:2, hs.pct(j).String = sprintf('%.3g%%', C{iT,j+3}/C{iT,3}*100); end
    hs.resp(1).Parent.XLim(2) = hs.ax.XLim(2);
    try
        for i = 1:3
            a = hs.fig.UserData.resp{iR}{i};
            set(hs.resp(i), 'XData', a, 'YData', ones(size(a)));
        end
    catch, set(hs.resp, 'XData', nan, 'YData', 1);
    end
    update_resp(hs);
end
hs.instnc.String = '';
hs.series.String = C{iT,1}; % in case hdr not saved
try s = hs.fig.UserData.hdr{iR}; catch, set_img(hs.img, inf(2)); return; end
if ~isfile(s.Filename), s.Filename = strrep(s.Filename, hs.rootDir, hs.backupDir); end

nIN = C{iT,3};
iIN = ceil(nIN/2); % start with middle Instance if avail
nam = sprintf('%s%06g.dcm', s.Filename(1:end-10), iIN);
if ~isfile(nam), iIN = 1; nam = s; end
hs.instnc.String = num2str(iIN);
set(hs.slider, 'Max', nIN, 'Value', iIN, 'UserData', s.Filename(1:end-10));
if nIN == 1, hs.slider.Visible = 'off';
else, set(hs.slider, 'SliderStep', [1 1]./(nIN-1)); hs.slider.Visible = 'on';
end
hs.ax.XLim(2) = numel(hs.dv.YData) + 0.5;
hs.series.String = seriesInfo(s);
try img = dicm_img(nam); catch, img = ones(2)*0.94; end
set_img(hs.img, img);

%% Load subj data to review
function loadSubj(h, ~)
hs = guidata(h);
[fname, pName] = uigetfile([hs.logDir '*.mat'], 'Select MAT file for a Patient');
if isnumeric(fname), return; end
load([pName '/' fname], 'T3');
hs.fig.UserData = T3.Properties.UserData;
DV = hs.fig.UserData.DV;
N = size(T3, 1);
C = flip(table2cell(T3), 1); C(:,6) = C(:,3);
dy = hs.ax.UserData;
for i = 1:N
    a = DV{N-i+1}; a = [0 a(~isnan(a))];
    if numel(a)<2, C(i,3:5) = {numel(a) [] []};
    else, C(i,3:5) = {numel(a) sum(a<dy(2)) sum(a<dy(3))};
    end
    try C{i,3} = hs.fig.UserData.nIN{N-i+1}; end
end
hs.table.Data = C;
s = hs.fig.UserData.hdr{end};
hs.subj.String = regexprep(s.PatientName, '[_\s]', '');
[hs.subj.UserData, nam] = fileparts(s.Filename);
hs.series.UserData = sscanf(nam, '%u_%u', 1:2);
tableCB(hs.table, struct('Indices', [1 1])); % show top series

%% close subj
function closeSubj(h, ~)
hs = guidata(h);
hs.table.Data = {};
hs.img.CData = ones(2)*0.94;
hs.fd.YData = nan; hs.dv.YData = nan;
c = {{}};
hs.fig.UserData = struct('FD', c, 'DV', c, 'hdr', c, 'resp', c, 'GM', c, 'nIN', c);
set([hs.subj hs.series hs.instnc], 'String', '');
set(hs.pct, 'Visible', 'off');

%% Re-do current subj: useful in case of error during a session
function redoSubj(h, ~)
hs = guidata(h);
subj = hs.subj.String;
if isempty(subj), return; end
if ~isfolder(hs.subj.UserData)
    fprintf(2, 'Image for %s deleted?\n', subj);
    return;
end
try delete([hs.logDir subj '*.mat']); catch, end
closeSubj(hs.fig)

%% Get reference vol info. Adapted from nii_moco.m
function p = refVol(img, pixdim)
d = size(img);
p.R0 = diag([pixdim 1]); % no need for real xform_mat here
p.R0(1:3, 4) = -pixdim .* (d/2); % make center voxel [0 0 0]

sz = pixdim;
if all(abs(diff(sz)/sz(1))<0.05) && sz(1)>2 && sz(1)<4 % 6~12mm
    sz = 3; % iso-voxel, 2~4mm res, simple fast smooth
else
    sz = 9 ./ sz'; % 9 mm seems good
end

% resample ref vol to isovoxel (often lower-res)
d0 = d-1;
dd = 4 ./ pixdim; % use 4 mm grid for alignmen
[i, j, k] = ndgrid(0:dd(1):d0(1)-0.5, 0:dd(2):d0(2)-0.5, 0:dd(3):d0(3)-0.5);
I = [i(:) j(:) k(:)]';
a = rng('default'); I = I + rand(size(I))*0.5; rng(a); % used by spm
V = smooth_mc(img, sz);
F = griddedInterpolant({0:d0(1), 0:d0(2), 0:d0(3)}, V, 'linear', 'none');
V0 = F(I(1,:), I(2,:), I(3,:)); % ref: 1 by nVox
I(4,:) = 1; % 0-based ijk: 4 by nVox
I = p.R0 * I; % xyz of ref voxels

% compute derivative to each motion parameter in ref vol
dG = zeros(6, numel(V0));
dd = 1e-6; % delta of motion parameter, value won't affect dG much
R0i = inv(p.R0); % speed up a little
for i = 1:6
    p6 = zeros(6,1); p6(i) = dd; % change only 1 of 6
    J = R0i * rigid_mat(p6) * I; %#ok<*MINV>
    dG(i,:) = F(J(1,:), J(2,:), J(3,:)) - V0; % diff now
end
dG = dG / dd; % derivative

% choose voxels with larger derivative for alignment: much faster
a = sum(dG.^2); % 6 derivatives has similar range
ind = a > std(a(~isnan(a)))/10; % arbituray threshold. Also exclude NaN
p.mean = mean(V0);
p.dG = dG(:, ind);
p.V0 = V0(ind);
p.mm = I(:, ind);
F.GridVectors = {0:d(1)-1, 0:d(2)-1, 0:d(3)-1};
p.F = F;
p.sz = sz;

%% motion correction to ref-vol. From nii_moco.m
function [m6, rst] = moco_estim(p, R)
mss0 = inf;
rst = R;
for iter = 1:64
    J = R * p.mm; % R_rst*J -> R0*ijk
    V = p.F(J(1,:), J(2,:), J(3,:));
    ind = ~isnan(V); % NaN means out of range
    dV = p.V0(ind) - V(ind);
    mss = dV*dV' / numel(dV); % mean(dV.^2)
    if mss > mss0, break; end % give up and use previous R
    rst = R; % accecpt only if improving
    if 1-mss/mss0 < 1e-6, break; end % little effect, stop

    a = p.dG(:, ind);
    p6 = (a * a') \ (a * dV'); % dG(:,ind)'\dV' estimate p6 from current R
    R = R * rigid_mat(p6); % inv(inv(rigid_mat(p6)) * inv(R_rst))
    mss0 = mss;
end

R = p.R0 * rst; % inv(R_rst / Rref)
m6 = -[R(1:3, 4)' atan2(R(2,3), R(3,3)) asin(R(1,3)) atan2(R(1,2), R(1,1))];

%% Translation (mm) and rotation (deg) to 4x4 R. Order: ZYXT
function R = rigid_mat(p6)
ca = cosd(p6(4:6)); sa = sind(p6(4:6));
rx = [1 0 0; 0 ca(1) -sa(1); 0 sa(1) ca(1)]; % 3D rotation
ry = [ca(2) 0 sa(2); 0 1 0; -sa(2) 0 ca(2)];
rz = [ca(3) -sa(3) 0; sa(3) ca(3) 0; 0 0 1];
R = rx * ry * rz;
R = [R p6(1:3); 0 0 0 1];

%% Simple gaussian smooth for motion correction, sz in unit of voxels
function out = smooth_mc(in, sz)
out = double(in);
if numel(unique(in))<5, return; end
if all(abs(diff(sz)/sz(1))<0.05) && abs(sz(1)-round(sz(1)))<0.05 ...
        && mod(round(sz(1)),2)==1
    out = smooth3(out, 'gaussian', round(sz)); % sz odd integer
    return; % save time for special case
end

d = size(in);
I = {1:d(1) 1:d(2) 1:d(3)};
n = sz/3;
if numel(n)==1, n = n*[1 1 1]; end
J = {1:n(1):d(1) 1:n(2):d(2) 1:n(3):d(3)};
intp = 'linear';
F = griddedInterpolant(I, out, intp);
out = smooth3(F(J), 'gaussian'); % sz=3
F = griddedInterpolant(J, out, intp);
out = F(I);

%% mask brain for better motion estimate
function img = brainMask(img, s)
persistent msk isDTI;
if nargin>1 % init msk from Volume 1
    isDTI = contains(s.ImageType, '\DIFFUSION');
    nii = nii_tool('init', logical(img)); % xposed, img=dicm_img(s)
    xform_mat = dicm2nii('', 'xform_mat', 'func_handle');
    [~, R, nii.hdr.pixdim(2:4)] = xform_mat(s); R = R(:,[2 1 3 4]); % LPS
    nii.hdr.sform_code = 1;
    nii.hdr.srow_x = -R(1,:); % RAS
    nii.hdr.srow_y = -R(2,:);
    nii.hdr.srow_z = R(3,:);
    if ~isDTI || isempty(msk)
        msk = nii;
        msk.img = EPI_mask(img);
        msk.img = smooth3(double(msk.img), 'guassian', 5) > 1e-6; % dilate
    end
    if isDTI % possibly EPI mask
        msk = nii_xform(msk, nii, [], 'nearest', false);
    end
end
img(~msk.img) = 0;
if ~isDTI, return; end % for EPI: only mask out
img = EPI_mask(img);

%% Quick brain mask for EPI or alike with low contrast
function M = EPI_mask(img)
B = smooth3(double(img), 'box', 5);
mn = median(B(B>mean(B(:)/8)));
mn = min(max(B(:)), mn*20);
[N, edges] = histcounts(img(:), linspace(0,mn,100));
N = movmean(N, 9);
mx = mean(N);
for i = 1:9
    p = islocalmin(N, 'MinProminence', mx);
    if any(p), break; end
    mx = mx / 2;
end
if ~any(p), M = true(size(B)); return; end % give up
B = B > edges(find(p,1))*0.7;
M = false(size(B));
c = round(nii_viewer('LocalFunc', 'img_cog', img));
M(c(1), c(2), c(3)) = true; n1 = 1;
while 1
    for i = 1:3
        M = (M<circshift(M, 1,i) | M<circshift(M,-1,i)) & B | M;
    end
    n2 = sum(M(:));
    if n2>n1, n1 = n2; else, break; end
end

%% new series or new subj: result saved as incoming_DICOM/RTMM_log/subj.mat
% The subj folders (yyyymmdd.PatientName.PatientID) default to ../incoming_DICOM/
% The dcm file names from Siemens push are always in format of
% 001_000001_000001.dcm. All three numbers always start at 1, and are continuous.
% First is study, second is series and third is instance.
function new = new_series(hs)
try setCountDown(hs); catch, end
QC_report(hs.rootDir);
f = hs.subj.UserData;
if ~isempty(f) % check new run for current subj
    iR = hs.series.UserData;
    if isfile(sprintf('%s/%03u_%06u_000001.dcm', f, iR+[0 1]))
        hs.series.UserData(2) = iR(2) + 1; new = true; return;
    elseif isfile(sprintf('%s/%03u_000001_000001.dcm', f, iR(1)+1))
        hs.series.UserData = [iR(1)+1 1];  new = true; return;
    end
end
dirs = dir([hs.rootDir '20*']); % check new subj
dirs(~isfile(strcat(hs.rootDir, {dirs.name}, '/001_000001_000001.dcm'))) = [];
new = false;
for i = numel(dirs):-1:1
    subj = regexp(dirs(i).name, '(?<=\d{8}\.).*?(?=\.)', 'match', 'once');
    subj = regexprep(subj, '[_\s]', '');
    if isfile([hs.logDir subj '.mat']), continue; end
    hs.subj.UserData = [hs.rootDir dirs(i).name];
    hs.series.UserData = [1 1];
    new = true; return;
end

% Move/Delete old subj folder right after mid-night
if ~isfile([hs.rootDir 'host.txt']) || mod(now,1) > 10/86400; return; end
dirs(now-[dirs.datenum]<2) = []; % keep for 2 days
for i = 1:numel(dirs)
    src = [hs.rootDir dirs(i).name];
    if ~isfolder(hs.backupDir), rmdir(src, 's');
    else, system(['cp -p -r ' src ' ' hs.backupDir ';rm -rf ' src]);
    end
end

%% Subfunction: get a parameter in CSA series ASC header: MrPhoenixProtocol
function val = asc_header(s, key, dft)
if nargin>2, val = dft; else, val = []; end
csa = 'CSASeriesHeaderInfo';
if ~isfield(s, csa) % in case of multiframe
    try s.(csa) = s.SharedFunctionalGroupsSequence.Item_1.(csa).Item_1; end
end
if isfield(s, 'Private_0029_1020') && isa(s.Private_0029_1020, 'uint8')
    str = char(s.Private_0029_1020(:)');
    str = regexp(str, 'ASCCONV BEGIN(.*)ASCCONV END', 'tokens', 'once');
    if isempty(str), return; end
    str = str{1};
elseif isfield(s, 'MrPhoenixProtocol') % X20A
    str = s.MrPhoenixProtocol;
elseif ~isfield(s, csa), return; % non-siemens
elseif isfield(s.(csa), 'MrPhoenixProtocol') % most Siemens dicom
    str = s.(csa).MrPhoenixProtocol;
elseif isfield(s.(csa), 'MrProtocol') % older version dicom
    str = s.(csa).MrProtocol;
else, return;
end

% tSequenceFileName  = ""%SiemensSeq%\gre_field_mapping""
expr = ['\n' regexptranslate('escape', key) '\s*=\s*(.*?)\n'];
str = regexp(str, expr, 'tokens', 'once');
if isempty(str), return; end
str = strtrim(str{1});

if strncmp(str, '""', 2) % str parameter
    val = str(3:end-2);
elseif strncmp(str, '"', 1) % str parameter for version like 2004A
    val = str(2:end-1);
elseif strncmp(str, '0x', 2) % hex parameter, convert to decimal
    val = sscanf(str(3:end), '%x', 1);
else % decimal
    val = sscanf(str, '%g', 1);
end

%% Wait till a file copy is complete
function s = dicm_hdr_wait(varargin)
tEnd = now + 1/86400; % wait up to 1 second
while 1
    s = dicm_hdr(varargin{:});
    try %#ok<*TRYNC> maybe too strict to be equal? Test indicates always equal
        if s.PixelData.Start+s.PixelData.Bytes == s.FileSize, return; end
    end
    if now>tEnd, s = []; return; end % give up
    pause(0.1);
end

%% User closing GUI: stop and delete timer
function closeFig(fh, ~)
hs = guidata(fh);
delete(fh);
try delete(hs.serial); catch, end
try hs.timer.StopFcn = ''; catch, end
try tObj = timerfindall; stop(tObj); delete(tObj); catch, end

%% menu callback for both DERIVED and _SBRef
function toggleChecked(h, ~)
if h.Checked == "on", h.Checked = 'off'; else, h.Checked = 'on'; end

%% Increase/Decrease image CLim
function setCLim(h, ~)
hs = guidata(h);
ax = hs.img.Parent;
if startsWith(h.Label, 'Increase'), ax.CLim(2) = ax.CLim(2)*0.8;
elseif startsWith(h.Label, 'Decrease'), ax.CLim(2) = ax.CLim(2)*1.2;
else, ax.CLim(2) = imgClim(hs.img.CData);
end

%% show series in nii_viewer
function view_3D(h, ~)
hs = guidata(h);
nams = dir([hs.slider.UserData '*.dcm']);
if isempty(nams), return; end
nams = strcat(nams(1).folder, '/', {nams.name});
nii = dicm2nii(nams, ' ', 'no_save');
nii_viewer(nii);

%% overlay series onto T1w or Scout if avail
function overlay(h, ~)
hs = guidata(h);
hdrs = hs.fig.UserData.hdr;
if isempty(hdrs), return; end
is3D = @(c)c.MRAcquisitionType=="3D" && ~contains(c.ImageType, 'DERIVED');
is3D = cellfun(is3D, hdrs);
if ~any(is3D), is3D = cellfun(@(c)c.MRAcquisitionType=="3D", hdrs); end
if sum(is3D)>1
    isT1 = cellfun(@(c)contains(c.SequenceName, 'fl3d1'), hdrs);
    isT1 = isT1 & is3D;
    if any(isT1), is3D = isT1; end
end
if ~any(is3D), view_3D(h); return; end % no T1, just show in nii_viewer
is3D = find(is3D, 1, 'last');
a = hdrs{is3D}.Filename(1:end-10);
nams = dir([a '*.dcm']);
if isempty(nams), nams = dir([strrep(a, hs.rootDir, hs.backupDir) '*.dcm']); end
nams = strcat(nams(1).folder, '/', {nams.name});
T1w = dicm2nii(nams, ' ', 'no_save');
nams = dir([hs.slider.UserData '*.dcm']);
nams = strcat(nams(1).folder, '/', {nams.name});
epi = dicm2nii(nams, ' ', 'no_save');
fh = nii_viewer(T1w, epi);
nii_viewer('LocalFunc', 'nii_viewer_cb', [], [], 'center', fh);

%% slider callback: show img if avail
function sliderCB(h, ~)
hs = guidata(h);
if isempty(h.UserData), return; end
h.Value = round(h.Value);
nam = sprintf('%s%06u.dcm', h.UserData, h.Value);
if ~isfile(nam), nam = strrep(nam, hs.rootDir, hs.backupDir); end
if ~isfile(nam), return; end
set_img(hs.img, dicm_img(nam));
hs.instnc.String = num2str(h.Value);

%% Timer StopFunc: Save result, start timer with delay, even after error
function saveResult(obj, ~)
hs = guidata(obj.UserData);
set([hs.menu hs.table hs.slider], 'Enable', 'on');
if size(hs.table.Data,1) > numel(hs.fig.UserData.FD) % new series to save?
    hs.fig.UserData.FD{end+1} = hs.fd.YData;
    hs.fig.UserData.DV{end+1} = hs.dv.YData;
    hs.fig.UserData.GM{end+1} = hs.img.UserData;
    hs.fig.UserData.resp{end+1} = {hs.resp.XData};
    T3 = cell2table(flip(hs.table.Data(:,[1 2 6]), 1), ...
        'VariableNames', {'Description' 'SeriesNumber' 'MeanFD'});
    T3.Properties.UserData = hs.fig.UserData;
    save([hs.logDir hs.subj.String], 'T3');
    hs.serial.UserData.send = false; % stop until asked again
end
if new_series(hs), obj.StartDelay = 0.1; else, obj.StartDelay = 5; end
start(obj);

%% Serial BytesAvail callback: update response: 1=missed, 2=incorrect, 3=correct
function serialRead(s, ~)
if s.BytesAvailable<1, return; end
b = fread(s, 1);
hs = guidata(s.UserData.fig);
if     b == '?', fwrite(s, uint8('RTMM')); return; % identity
elseif b == 'P', fwrite(s, uint8(hs.subj.String)); return; % PatientName
elseif b == 'T', fwrite(s, uint8(strjoin(hs.MMSS.String))); return; % count down
elseif b == 'M', s.UserData.send = true; return; % start to send motion info
elseif b == 'Q' % stim computer asks to stop scan
    [~, ~] = system(['touch "' hs.rootDir 'StopScan"']); return;
elseif b<1 || b>3, return; % ignore for now
end
if hs.timer.StartDelay>1, return; end % not during a series
x = find(~isnan(hs.dv.YData), 1, 'last');
if isempty(x), x = 0; end
hs.resp(b).XData(end+1) = x + 1; hs.resp(b).YData(end+1) = 1;
update_resp(hs);

%% update response text
function update_resp(hs)
n = cellfun(@numel, {hs.resp.XData}) - 1;
h = hs.resp(1).Parent.Title;
if ~any(n>0), h.Visible = 'off'; return; end
h.Visible = 'on';
h.String = "Missed " + n(1) + ", \color{red}Incorrect " + n(2) + ...
    ", \color[rgb]{0 0.8 0}Correct " + n(3) + ", \color{blue}Total " + sum(n);

%% start countdown
function setCountDown(hs)
nam = [hs.rootDir 'SyngoMeas'];
if ~isfile(nam), return; end
c0 = fileread(nam); pause(0.2); c = fileread(nam);
if ~isequal(c0, c), pause(1); c = fileread(nam); end
delete(nam);
% From scanner: "RunStartTime" "ProtocolName" TotalScanTimeSec
c = regexp(c, '"(.*?)" "(.*?)" (\d+)', 'tokens', 'once');
if isempty(c), return; end % MeasFinished?
tStart = datenum(c{1}, 'yyyy-mm-dd HH:MM:SS,fff');
dClock = str2double(fileread([hs.rootDir 'dClock']))/86400;
tFnsh = tStart - dClock + str2double(c{3})/86400;
if tFnsh-now < 2/86400, return; end
hs.countDown.UserData = tFnsh;
if numel(c{2})>24, c{2} = [c{2}(1:16) '...' c{2}(end-3:end)]; end
hs.MMSS.String = {c{2} ''};
if hs.countDown.Running=="off", start(hs.countDown); end

%% timer func to show scanning time
function countDown(tObj, ~, hs)
t = tObj.UserData - now;
if t<1/86400, stop(tObj); hs.MMSS.String = {'' ''}; return; end
hs.MMSS.String{2} = ['Scanning ' datestr(t, 'MM:SS')];
nam = hs.rootDir + "SyngoMeas";
if ~isfile(nam), return; end
if ~startsWith(fileread(nam), "Finished"), return; end
delete(nam);
stop(tObj);
hs.MMSS.String{2} = ['Finished ' datestr(t, 'MM:SS')];

%% get CLim for dicom img
function mx = imgClim(img)
im = double(img(:));
im = im(im>max(im)/10);
mx = mean(im) + 2*std(im);

%% Create QC report
function QC_report(rootDir)
nam = dir([rootDir 'closed_*']);
if isempty(nam) || isfile([rootDir 'EyelinkRecording.mat']), return; end
nam = [rootDir nam(1).name];
done = onCleanup(@()movefile(nam, strrep(nam, 'closed_', 'done_')));
rmQC = onCleanup(@()delete('./tmp_QC_*.pdf'));
subj = regexp(nam, '(?<=closed_)\d{4}\w{2}$', 'match', 'once');
if isempty(subj), return; end
load([rootDir 'RTMM_log/' subj '.mat'], 'T3');
uDat = T3.Properties.UserData;

close all; delete('./tmp_QC_*.pdf');
fig = figure('Position', [10 30 [8.5 11]*96], 'Units', 'normalized');
set(fig, 'Color', 'w', 'PaperUnits', 'normalized', 'PaperPosition', [0 0 1 1]);
ax = axes(fig, 'Position', [0.01 0.955 0.98 0.045], 'Visible', 'off');
try imshow('./logo.png', 'Parent', ax); ax.HandleVisibility = 'off'; catch, end

layout = getpref('nii_viewer_para', 'layout');
if layout ~= 1
    setpref('nii_viewer_para', 'layout', 1);
    cln = onCleanup(@()setpref('nii_viewer_para', 'layout', layout));
end

ax = axes(fig, 'Position', [0.1 0.92 0.8 0.03], 'Visible', 'off');
text(ax, 0.5, 1, subj, 'FontSize', 18, 'HorizontalAlignment', 'center');
s = uDat.hdr{1};
dat = datestr(datenum(dicmTime(s,1), 'yymmdd'), 'dddd mmm dd, yyyy');
text(ax, 0.5, 0, dat, 'FontSize', 12, 'HorizontalAlignment', 'center');
tbl = cell(1, 5);
dict = dicm_dict('', {'AcquisitionDateTime' 'AcquisitionDate' ...
    'AcquisitionTime' 'SeriesNumber' 'SeriesDescription'});
for i = 1:99
    nams = dir([s.Filename(1:end-17) sprintf('%06i_',i) '*.dcm']);
    if isempty(nams), break; end
    s = dicm_hdr([nams(1).folder '/' nams(1).name], dict);
    d = dicmTime(s); d = d(7:8)+":"+d(9:10)+":"+d(11:12);
    try a = T3.MeanFD{T3.SeriesNumber == s.SeriesNumber}; catch, a = []; end
    tbl(i,:) = {s.SeriesNumber d numel(nams) s.SeriesDescription a};
end
tbl = cellfun(@num2str, tbl, 'UniformOutput', false); % for left-align
vName = {'SeriesNumber' 'Time' 'TotalInstances' 'Description' 'meanFD'};
h = uitable(fig, 'Units', 'normalized', 'Position', [0.04 0.1 0.9 0.8], ...
    'FontSize', 12, 'RowName', [], 'ColumnName', vName, 'Data', tbl, ...
    'ColumnWidth', {96 96 108 342 82});
if i>40, h.FontSize = max(8, fix(450/i)); end
newPage(fig); y = 0.95;

el = {}; clear EL;
nam = dir([rootDir 'RTMM_log/' subj '_*.edf']);
while ~isempty(nam)
    tm = []; tg = []; pa = [];
    for i = numel(nam):-1:1 % in case of multiple files
        [~, a] = evalc("edfmex('"+nam(i).folder+"/"+nam(i).name+"')");
        tm = [a.FSAMPLE.time tm]; %#ok
        tg = [a.FSAMPLE.buttons tg]; %#ok
        pa = [a.FSAMPLE.pa pa]; %#ok
    end
    if numel(unique(pa(1,:)))>9, pa = pa(1, :);
    elseif numel(unique(pa(2,:)))>9, pa = pa(2, :);
    else, break;
    end
    i = find(arrayfun(@(c)startsWith(char(c.message), 'DICM_'), a.FEVENT), 1, 'last');
    a = double(tm)/1000 - double(a.FEVENT(i).sttime)/1000 + ...
        sscanf(a.FEVENT(i).message, 'DICM_secs=%g');
    EL.Hz = 100; % resample pa to 100 Hz
    EL.t0 = a(1); % recording start time in Syngo secs of the day
    EL.tg = a(diff([0 bitget(tg, 5)])>0); % trigger time
    EL.pa = interp1(a, pa/55^2, a(1):1/EL.Hz:a(end)); % convert roughly to degree
    break;
end

for i = 1:numel(uDat.hdr)
    s = uDat.hdr{i};
    series = sprintf('%s (Series %g)', s.ProtocolName, s.SeriesNumber);
    if contains(s.SequenceName, {'epfid2d' 'mbPCASL' 'ep_b0'}) % EPI/ASL/Diff
        fd = uDat.FD{i};
        N = find(~isnan(fd), 1, 'last');
        if isempty(N) || N<11 || all(fd==0), continue; end % skip slice check etc
        if y<0.35, newPage(fig); y = 0.95; end
        if exist('EL', 'var') && contains(s.SequenceName, 'epfid2d')
            t0 = mod(datenum(dicmTime(s), 'yymmddHHMMSS.fff'), 1) * 86400;
            [mi, j] = min(abs(EL.tg - t0));
            if mi<5 % recording started late if mi too large
                TR = asc_header(s, 'alTR[0]')*1e-6;
                while j>1 && abs(diff(EL.tg([j-1 j]))-TR)<0.1, j = j-1; end
                % fprintf('%4.2f\n', EL.tg(j)-t0); % normally <1
                j = floor((EL.tg(j)-EL.t0) * EL.Hz);
                nP = round(N * TR * EL.Hz);
                el{end+1} = {EL.pa(j+(1:nP)) series}; %#ok save for users
                y = y - 0.4;
                ax0 = axes(fig, 'Position', [0.1 y+0.27 0.8 0.08]);
                plot(ax0, el{end}{1}, 'b');
                set(ax0, 'XLim', [1 nP], 'XTick', []);
                ylabel(ax0, 'Pupil Size', 'Color', 'b');
            end
        else, y = y - 0.32; clear ax0;
        end
        ax = axes(fig, 'Position', [0.1 y+0.19 0.8 0.08]);
        gm = uDat.GM{i}(1:N);
        gm = (gm/mean(gm) - 1)*100; rg = ceil(std(gm)); % global mean
        % gm = [0; diff(gm)]/mean(gm) * 100; % delta GM
        plot(ax, gm, '.-m');
        set(ax, 'YTick', [], 'YLim', [-rg rg], 'XLim', [0 N+1], 'XTick', []);
        ylabel(ax, 'GM (%)', 'Color', 'm');
        try ax = ax0; catch, end
        title(ax, series, 'Interpreter', 'none', 'FontSize', 12);
        % ylabel(ax, [char(916) 'GM (%)']);
        ax = axes(fig, 'Position', [0.1 y+0.03 0.8 0.16]);
        yyaxis right; plot(ax, fd(1:N), '.-');
        ylabel(ax, 'FD (mm)'); xlabel(ax, 'Volume Number');
        set(ax, 'YTick', 0:0.6:1.8, 'YLim', [0 1.8], 'XLim', [0 N+1]);
        str = sprintf('meanFD=%.2g', mean(fd(2:N)));
        text(ax, 0.8, 0.9, str, 'Units', 'normalized', 'Color', [0.85 0.32 0.1]);
        yyaxis left; plot(ax, uDat.DV{i}(1:N), '.-');
        ylabel(ax, 'DVARS'); set(ax, 'YTick', 0:0.12:0.36, 'YLim', [0 0.3]);
    elseif contains(s.SequenceName, {'tfl3d' 'spc' 'tse2d' 'fm2d' 'epse2d'}) % T1/T2/fmap
        nams = dir([s.Filename(1:end-10) '*.dcm']);
        nams = strcat(nams(1).folder, '/', {nams.name});
        fh = nii_viewer(dicm2nii(nams, ' ', 'no_save'));
        nii_viewer('LocalFunc', 'nii_viewer_cb', [], [], 'center', fh);
        hsV = guidata(fh);
        drawnow; F = getframe(fh, hsV.frame.Position); close(fh); img = F.cdata;
        sz = size(img); sz = sz([2 1]) ./ [8.5 11]/96;
        if sz(1)>0.8, sz = sz / sz(1) * 0.8; end
        y = y - sz(2) - 0.08;
        if y<0.02 && y+sz(2)*0.2>0.02, sz = sz * 0.8; y = y + 0.2*sz(2); end
        if y<0.02, newPage(fig); y = 0.92-sz(2); end
        ax = axes(fig, 'Position', [(1-sz(1))/2 y+0.02 sz], 'Visible', 'off');
        imshow(img, 'Parent', ax);
        title(ax, series, 'Interpreter', 'none', 'FontSize', 12);
    end
end
newPage(fig); close(fig);

if ~isempty(el), save([rootDir 'RTMM_log/' subj '_eye.mat'], 'el'); end
pdfNam = rootDir+"RTMM_log/"+subj+"_"+dicmTime(s,1)+"_QC.pdf";
setenv('LD_LIBRARY_PATH', getenv('PATH'));
[~, ~] = system("pdfunite ./tmp_QC_*.pdf "+pdfNam); % for Linux

%% print fig to a new PDF, called by QC_report()
function newPage(fig)
nam = dir('./tmp_QC_*.pdf');
if isempty(nam), i = 1; else, i = str2double(nam(end).name(8:10)) + 1; end
print(fig, sprintf('./tmp_QC_%03i.pdf', i), '-dpdf');
clf(fig);

%% 3D to mos: for display purpose
function mos = vol2mos(vol)
[nr, nc, nSL] = size(vol);
nMos = ceil(sqrt(nSL));
mos = zeros([nr nc]*nMos, 'like', vol);
for i = 0:nMos-1
    r = i*nr + (1:nr);
    for j = 1:nMos
        iSL = i*nMos + j;
        if iSL>nSL, return; end
        c = (j-1)*nc + (1:nc);
        mos(r,c) = vol(:,:,iSL);
    end
end

%% return dicm dateTime in yymmddHHMMSS.fff
function dt = dicmTime(s, dateOnly)
if isfield(s, 'AcquisitionDateTime')
    dt = s.AcquisitionDateTime(3:end);
elseif isfield(s, 'AcquisitionDate')
    dt = [s.AcquisitionDate(3:8) s.AcquisitionTime];
else
    a = dir(s.Filename); dt = datestr(a.datenum, 'yymmddHHMMSS.fff');
end
if nargin>1 && dateOnly, dt = dt(1:6); end
%%

author: xiangruili


text:
classdef SmoothingExponentComputerFactory < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            switch cParams.type
                case 'Optimal'
                   obj = SmoothingExponentComputerOptimal(cParams);
                case 'Given'
                   obj = SmoothingExponentComputerGiven(cParams);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan