text:
function [ f,resL2,qualMeasOut] = OS_AwPCSD(proj,geo,angles,maxiter,varargin)
%OS_AwPCSD solves the reconstruction problem using adaptive-weighted
%projection-controlled steepest descent method
%
%   OS_AwPCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using
%   the projection data PROJ taken over ALPHA angles, corresponding to the
%   geometry described in GEO, using NITER iterations.
%
%   OS_AwPCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter for the SART iterations.
%                  Default is 1
%
%   'lambdared':   Reduction of lambda.Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%       'init':    Describes diferent initialization techniques.
%                   •  'none'     : Initializes the image to zeros(default)

%                   •  'FDK'      : intializes image to FDK reconstrucition
%
%   'TViter':      Defines the amount of TV iterations performed per SART
%                  iteration. Default is 20
%
%   'maxL2err'     Maximum L2 error to accept an image as valid. This
%                  parameter is crucial for the algorithm, determines at
%                  what point an image should not be updated further.
%                  Default is 20% of the FDK L2 norm.
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%   'delta'        Defines Parameter to control the amount of smoothing
%                  for pixels at the edges. A large 'delta' is not able to
%                  differentiate image gradients at different pixels. A
%                  small 'delta' give low weights to almost every pixels,
%                  making the algorithm inefficient in removing noise or
%                  straking artifacts. Default is -0.00055
%   'QualMeas'     Asks the algorithm for a set of quality measurement
%                  parameters. Input should contain a cell array of desired
%                  quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                  These will be computed in each iteration.
%   'BlockSize':   Sets the projection block size used simultaneously. If
%                  BlockSize = 1 OS-SART becomes SART and if  BlockSize = length(alpha)
%                  then OS-SART becomes SIRT. Default is 20.
% 'OrderStrategy'  Chooses the subset ordering strategy. Options are
%                  'ordered' :uses them in the input order, but divided
%                  'random'  : orders them randomply
%                  'angularDistance': chooses the next subset with the
%                                     biggest angular distance with the ones used.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri and Manasavee Lohvithee
%--------------------------------------------------------------------------
%% parse inputs
[beta,beta_red,f,ng,verbose,epsilon,delta,blocksize,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);

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

[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy);
%
% angles=cell2mat(alphablocks);
% index_angles=cell2mat(orig_index);

if measurequality
    qualMeasOut=zeros(length(QualMeasOpts),maxiter);
end

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

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

% Projection weight, W
W=computeW(geo,angles,gpuids);

% Back-Projection weigth, V
V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids);

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

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

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

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

    %Enforcing ART along all projections if squared delta_p > epsilon
    if (delta_p^2)>epsilon
        for jj=1:length(alphablocks)
            if size(offOrigin,2)==size(angles,2)
                geo.offOrigin=offOrigin(:,jj);
            end
            if size(offDetector,2)==size(angles,2)
                geo.offDetector=offDetector(:,jj);
            end
            if size(rotDetector,2)==size(angles,2)
                geo.rotDetector=rotDetector(:,jj);
            end
            if size(DSD,2)==size(angles,2)
                geo.DSD=DSD(jj);
            end
            if size(DSO,2)==size(angles,2)
                geo.DSO=DSO(jj);
            end

            f=f+beta* bsxfun(@times,1./sum(V(:,:,jj),3),Atb(W(:,:,jj).*(proj(:,:,orig_index{jj})-Ax(f,geo,alphablocks{:,jj},'gpuids',gpuids)),geo,alphablocks{:,jj},'gpuids',gpuids));

        end
    end

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

    if measurequality
        qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts);
    end

    % Compute L2 error of actual image. Ax-b
    dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2');

    % Compute change in the image after last SART iteration
    dp_vec=(f-f0);

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

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

    % Reduce SART step
    beta=beta*beta_red;

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

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

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

    if computeL2
        geo.offOrigin=offOrigin;
        geo.offDetector=offDetector;
        resL2(ii)=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2');                       % Compute error norm2 of b-Ax
        % If the error is not minimized.
        if  iter~=1 && resL2(ii)>resL2(ii-1)
            if verbose
                disp(['Convergence criteria met, exiting on iteration number:', num2str(iter)]);
            end
            return;
        end

    end

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

end
end

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

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

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

            else
                if strcmp(val,'FDK')
                    f0=FDK(proj, geo, angles);
                else
                    error('TIGRE:MLEM:InvalidInput','Invalid init')
                end
            end
            % Number of iterations of TV
            %  =========================================================================
        case 'tviter'
            if default
                ng=20;
            else
                ng=val;
            end
            %  Maximum L2 error to have a "good image"
            %  =========================================================================
        case 'maxl2err'
            if default
                epsilon=im3Dnorm(FDK(proj,geo,angles),'L2')*0.2; %heuristic
            else
                epsilon=val;
            end
            %  =========================================================================
            %  Parameter to control the amount of smoothing for pixels at the
            %  edges
            %  =========================================================================
        case 'delta'
            if default
                delta=-0.005;
            else
                delta=val;
            end
            %  =========================================================================
            %  Block size for OS-SART
            %  =========================================================================
        case 'blocksize'
            if default
                block_size=20;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:OS_AwASD_POCS:InvalidInput','Invalid BlockSize')
                end
                block_size=val;
            end
            %  Order strategy
            %  =========================================================================
        case 'orderstrategy'
            if default
                OrderStrategy='random';
            else
                OrderStrategy=val;
            end
            %  =========================================================================
            %  Image Quality Measure
            %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:OS_AwPCSD:InvalidInput','Invalid quality measurement parameters');
                end
            end
            %  Non negative
            %  =========================================================================
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
            %  GPU Ids
            %  =========================================================================
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:OS_AwPCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_AwPCSD()']);

    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function b=rowswap(a,i,j)
%i and j are row number%
a([i,j],:)=a([j,i],:);
b=a;

author: Sable

repo: mcbench-benchmarks


text:
function indxV = getAssociatedStr(uid)
%"getAssociatedStr"
%   Returns a structure index for corresponding structure uid passed. uid can be a
%   cell array of strUID's in that case output will be an array of structure indices
%
%	DK 12/07/2006
%
%Usage:
%   indxV = getAssociatedStr(uid)

% copyright (c) 2001-2008, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).

global planC
indexS = planC{end};

if length(planC{indexS.structures})< 1
    indxV = [];
    return
end

[strUID{1:length(planC{indexS.structures})}] = deal(planC{indexS.structures}.strUID);

if length(uid)==0
    indxV = [];
else
    [jnk,indxV] = ismember(uid,strUID);
    if indxV == 0;
        indxV = [];
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function jaccard = Jaccard_Index(SEG, GT)
    % SEG, GT are the binary segmentation and ground truth areas, respectively.
    % jaccard index
    jaccard = double(sum(uint8(SEG(:) & GT(:)))) / double(sum(uint8(SEG(:) | GT(:))));
end

repo: cluster


text:
 function varargout = subsref_mm(ob, args)
%function varargout = subsref_mm(ob, args)
% handle subscript references like ob.ref or ob.ref(arg[s])
% Copyright 2006-1-19, Jeff Fessler, University of Michigan
% 2011-10-25 modified by Madison McGaffin to support multiple outputs

st = struct(ob);
varargout = cell(nargout,1);

if args(1).type ~= '.'
	error 'only ob.data and ob.meth(args) done'
end

%
% ob.name
%
name = args(1).subs;
if isfield(st, name)
	if length(args) == 1
%		out = st.(name);
		varargout = wrap_vararg(nargout, 'st.(name)', st, ob, args);
	else
%		out = subsref(st.(name), args(2:end));
		varargout = wrap_vararg(nargout, ...
			'subsref(st.(args(1).subs), args(2:end));', ...
			st, ob, args );
	end

%
% ob.name, where name is in data field
%
elseif isfield(st.data, name)
%	out = st.data.(name);
	varargout = wrap_vararg(nargout, 'st.data.(args(1).subs);', st, ob, args);

	if length(args) > 1

		% ob.name{:} does not work properly because there seems
		% to be no way to return an arbitrary number of arguments
		% (in a list) to the caller.  this is unfortunate because
		% a.b = {1,2} and a.b{:} works for ordinary structures.
% todo: see this solution involving numel()
% http://www.mathworks.com/support/solutions/data/1-1ABOD.html?solution=1-1ABOD
% n.b., link broken; will need to rewrite this for multiple return arguments.

%		if iscell(out) && ...
		if length(varargout) > 0 && iscell(varargout{1}) && ...
		(isequal(args(2).subs, {':'}) || length(args(2).subs{1}) > 1)
			warn 'instead of strum.data{:} etc., use two steps:'
			warn 'tmp = strum.data; tmp{:} or tmp{1:end} etc.'
			fail 'strum.data{:} etc. unsupported thanks to matlab'
		end

		% fix: strum1.strum2.plot() does not work here
%		out = subsref(out, args(2:end));
		varargout = {subsref(varargout{1}, args(2:end))};
	end

%
% ob.name or ob.name(arg[s]), where name is in meth field
%
elseif isfield(st.meth, name)
	fun = st.meth.(name); % function handle for method
	if length(args) == 1
%		if isfreemat || nargout(fun) % freemat: does not like nargout()
		if isfreemat || nargout % freemat: does not like nargout()
%			out = fun(ob);
			varargout = wrap_vararg(nargout, 'fun(ob);', st, ob, ...
				args, fun);
		else
			fun(ob);
		end
	elseif length(args) == 2
%		if isfreemat || nargout(fun) % freemat: does not like nargout()
		if isfreemat || nargout % freemat: does not like nargout()
%			out = fun(ob, args(2).subs{:});
			varargout = wrap_vararg(nargout, 'fun(ob, args(2).subs{:});', ...
					st, ob, args, fun );
		else
			fun(ob, args(2).subs{:});
		end
	else
		printm('unknown subs')
		keyboard
	end

else
	disp(st)
	fail(['unknown field name: ' name])
end

end % subsref

% wrap_vararg()
% return cell array with each output
function out = wrap_vararg(nout, script, st, ob, args, fun)

switch nout
case 0
	eval(script);
	out = {};
case 1
	v = eval(script);
	out = {v};
case 2
	[v1, v2] = eval(script);
	out = {v1, v2};

otherwise % general case using eval
	head = sprintf('v%d, ', 1:(nout-1));
	full = sprintf('[ %s vlast ] = %s', head, script);
	out = cell(nout,1);
	for i=1:(nout-1)
		out{i} = eval(sprintf('v%d', i));
	end
	out{end} = vlast;
end

end % wrap_vararg

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [mtexFig,newFigure] = newMtexFigure(varargin)
% set up a new plotting figure
%
% Syntax
%
%   newMtexFigure('figSize','large')
%   newMtexFigure('layout',[2,3])
%
% Options
%  figSize - huge, large, normal, small, tiny
%  layout  - [numberRows, numberColumns]
%

mtexFig = gcm;

% check hold state
newFigure = isempty(mtexFig) || check_option(varargin,'newFigure') || ...
  (strcmp(getHoldState,'off') && ~check_option(varargin,{'hold','parent','add2all'}) ...
  && ~isempty(get(mtexFig.currentAxes,'Children')));

% check tag
if ~newFigure && check_option(varargin,'ensureTag') && ...
    ~any(strcmpi(get(gcf,'tag'),get_option(varargin,'ensureTag')))
  newFigure = true;
  warning('MTEX:newFigure','Plot type not compatible to previous plot! I''going to create a new figure.');
end

% check appdata
ad = get_option(varargin,'ensureAppdata');
if ~newFigure
  try
    for i = 1:length(ad)
      if isappdata(gcf,ad{i}{1}) && ~isempty(ad{i}{2})
        ad_stored = getappdata(gcf,ad{i}{1});
        if isa(ad{i}{2},'symmetry')
          ad{i}{2} = ad{i}{2}.Laue;
          ad_stored = ad_stored.Laue;
        end
        if ~all(ad_stored == ad{i}{2})
          newFigure = true;
          if ishold
            warning('MTEX:newFigure','Plot properties not compatible to previous plot! I''going to create a new figure.');
          end
          break
        end
      elseif ~isappdata(gcf,ad{i}{1})
        newFigure = true;
        if ishold
          warning('MTEX:newFigure','Plot properties not compatible to previous plot! I''going to create a new figure.');
        end
        break
      end
    end
  catch %#ok<CTCH>
    newFigure = true;
  end
end

% set up a new figure
if newFigure

  if check_option(varargin,'parent')

    fig = get_option(varargin,'parent');
    mtexFig.gca = fig;
    while ~isempty(fig) && ~strcmp('figure', get(fig,'type'))
      fig = get(fig,'parent');
    end
    mtexFig.parent = fig;
    newFigure = false;

  else

    mtexFig = mtexFigure(varargin{:});

    % set tag
    if check_option(varargin,'ensureTag','char')
      set(gcf,'tag',get_option(varargin,'ensureTag'));
    end

    % set appdata
    if check_option(varargin,'ensureAppdata')
      for i = 1:length(ad)
        setappdata(gcf,ad{i}{1},ad{i}{2})
      end
    end
  end
else % use an existing figure

  % get existing mtexFigure
  mtexFig = getappdata(gcf,'mtexFig');

  holdState = getHoldState;
  % distribute hold state over all axes
  for i=1:numel(mtexFig.children)
    hold(mtexFig.children(i),holdState);
  end

  % set current axis
  if check_option(varargin,'parent')
    p =  get_option(varargin,'parent');
    if isgraphics(p,'axes')
      mtexFig.currentAxes = p;
    else
      mtexFig.currentAxes = get(p,'parent');
    end
  elseif check_option(varargin,'add2all')
    mtexFig.currentId = 1;
  end
end

if nargout == 0, clear mtexFig; end

end

author: mtex-toolbox

repo: mtex

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


text:
clear;
clc;
close all;
format short

eul_312 = [30, 40, 50]; % [pitch(绕X轴) roll(绕Y轴)  yaw(绕Z轴)]

fprintf("按312顺序(先转Z-然后X-最后Y)旋转，其中X,Y,Z旋转角度为%.3f° %.3f° %.3f° 得到对应的坐标变换矩阵：\n", eul_312(1), eul_312(2), eul_312(3));

eul_312rad = deg2rad(eul_312);
[Cb2n_312, ~] = ch_eul2m(eul_312rad);
Cb2n_312

fprintf("将Cb2n_312转回欧拉角:\n");
[eul_312, ~]  = ch_m2eul(Cb2n_312);
eul_312 = rad2deg(eul_312);

fprintf("坐标变换矩阵转欧拉角:%.3f°(Pitch) %.3f°(Roll) %.3f°(Yaw)\n", eul_312(1), eul_312(2), eul_312(3));

eul_321 = [30, 40, 50]; % [roll(绕X轴) pitch(绕Y轴)  yaw(绕Z轴)]

fprintf("\n按321顺序(先转Z-然后X-最后Y)旋转，其中X,Y,Z旋转角度为%.3f° %.3f° %.3f° 得到对应的坐标变换矩阵：\n", eul_321(1), eul_321(2), eul_321(3));

eul_321rad = deg2rad(eul_321);
[~, Cb2n_321] = ch_eul2m(eul_321rad);
Cb2n_321

fprintf("将Cb2n_312转回欧拉角:\n");
[~, eul_321]  = ch_m2eul(Cb2n_321);
eul_321 = rad2deg(eul_321);

fprintf("坐标变换矩阵转欧拉角:%.3f°(Roll) %.3f°(Pitch) %.3f°(Yaw)\n", eul_321(1), eul_321(2), eul_321(3));

author: yandld


text:
function lagrange_nd_test07 ( )

%*****************************************************************************80
%
%% LAGRANGE_ND_TEST07 tests LAGRANGE_ND_COMPLETE in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_ND_TEST07\n' );
  fprintf ( 1, '  LAGRANGE_COMPLETE determines\n' );
  fprintf ( 1, '  the Lagrange interpolating polynomials L(x)\n' );
  fprintf ( 1, '  for ND points in D dimensions, assuming that\n' );
  fprintf ( 1, '  the number of points exactly coincides with\n' );
  fprintf ( 1, '  R = Pi(D,N), the number of monomials of degree N or less\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data points are the grid nodes of a tetrahedron.\n' );

  d = 3;
  n = 2;
  r = mono_upto_enum ( d, n );
  nd = 10;
  xd = [ 0.0,  0.0,  0.0; ...
         1.0,  0.0,  0.0; ...
         2.0,  0.0,  0.0; ...
         0.0,  1.0,  0.0; ...
         1.0,  1.0,  0.0; ...
         0.0,  2.0,  0.0; ...
         0.0,  0.0,  1.0; ...
         1.0,  0.0,  1.0; ...
         0.0,  1.0,  1.0; ...
         0.0,  0.0,  2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension D = %d\n', d );
  fprintf ( 1, '  Maximum degree N = %d\n', n );
  fprintf ( 1, '  Number of monomials R = %d\n', r );
  fprintf ( 1, '  Number of data points ND = %d\n', nd );

  r8mat_transpose_print ( d, nd, xd, '  Data points XD:' );

  [ po, pc, pe ] = lagrange_complete ( d, n, r, nd, xd );
%
%  Print the polynomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lagrange polynomials for XD data points:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nd
    o = po(i);
    label = sprintf ( '  P(%d)(x) =', i );
    polynomial_print ( d, o, pc(i,1:o), pe(i,1:o), label );
  end
%
%  Evaluate the polynomials at XD.
%
  value = zeros ( nd, nd );

  for j = 1 : nd
    o = po(j);
    label = sprintf ( '  P(%d)(x) =', j );
    value(1:nd,j) = polynomial_value ( d, o, pc(j,1:o), pe(j,1:o), nd, xd );
  end

  err = r8mat_is_identity ( nd, value );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of Lagrange matrix error = %g\n', err );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = remap(F,old,new)

F = flatten(F);
for i = 1:length(F.clauses)
    X = F.clauses{i}.data;
    X = sdpvarremap(X,old,new);
    F.clauses{i}.data = X;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% OPTI Toolbox AMPL Demo
%
% This file loads a number of supplied examples and shows how to solve
% them using the OPTI Toolbox. You should read and complete Basic_demo.m
% BEFORE running the below examples.
%
% The underlying parser uses Netlib's AMPL Solver Library (ASL) Interface.
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)

%% Loading an AMPL Problem
% OPTI Toolbox is supplied with a number of example AMPL problems ranging
% from LP to NLP including continuous and discrete problems. To load an
% AMPL model it must be in .NL format, see the supplied user's guide page
% on AMPL Interfacing in order to generate the .NL file from your model. To
% load an AMPL problem, simply use the command below. Returned will be a
% optiprob structure containing the data in the file.

prob = amplRead('diet.nl')

%% Example 1 - Solving a Loaded AMPL Problem
% Solving a loaded AMPL problem is simple, just pass it to the opti
% constructor and call solve. You will note that the command amplRead will
% automatically determine the problem type (LP, QP, NLP, etc) based on
% properties supplied by the ASL.

Opt = opti(prob)

x = solve(Opt)

%% Loading an AMPL Nonlinear Program
% The AMPL format can easily specify nonlinear problems unlike the LP or
% MPS formats.

prob = amplRead('ch3.nl')

%% Example 2 - Solving a Nonlinear Problem from AMPL
% In order to solve a nonlinear problem the model .nl file
% remains open after the amplRead function is called to allow callbacks for
% the objective, gradient, etc. Calling solve() will automatically close
% the interface for linear and quadratic problems, however you must manually
% call close for nonlinear problems. Alternatively it will be automatically
% called once Matlab is closed to clean up ASL memory.

Opt = opti(prob)

x = solve(Opt)

asl('close')

%% Loading an Integer Problem
% AMPL will reorder the variables from the original model based on internal
% rules (see the supplied pdfs). Based on this the ASL indicates which
% variables have binary or integer constraints, which are passed through
% the amplRead interface:

prob = amplRead('multmip1.nl')

%% Example 3 - Solving an AMPL Integer Problem
% As the optiprob structure can fully define an integer problem, no change
% is required to build and solve the MIP

Opt = opti(prob)

[x,fval] = solve(Opt);
fval

%% A Note On Hessians
% As of OPTI v1.58 the Hessian of an NLP will automatically be added to the
% problem structure. However note it is the Hessian of the Lagrangian, with
% the calling form:
%
%   H = hessian(x,sigma,lambda)
%
% This means the Hessian includes second derivatives of both the objective
% and the constraints. Sigma is a scalar scaling factor on the objective
% derivatives only, and lambda are the lagrange multipliers at the current
% point x (scaling the constraint derivatives).

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [P,E,logw,acc,traj] = spm_mci_ais_single (mcmc,M,U,Y)
% Produce a single independent sample using AIS
% FORMAT [P,E,logw,acc,traj] = spm_mci_ais_single (mcmc,M,U,Y)
%
% mcmc      Sampling settings
% M         Model structure
% U         Input structure
% Y         Data
%
% P         [Np x 1] sample
% E         Negative log joint
% logw      Contribution to model evidence
% acc       acc(j) is acceptance rate at temperature j
% traj      traj(p,j) is value of parameter p at temperature j
%           (only set if mcmc.rec_traj=1)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_ais_single.m 6548 2015-09-11 12:39:47Z will $

beta=mcmc.beta;
nprop=mcmc.nprop;
prop=mcmc.prop;
scale=mcmc.scale;

J=length(beta);

if isstruct(M.pC)
    pC = full(diag(spm_vec(M.pC)));
else
    pC = M.pC;
end
% prior cov in subspace
pC=M.V'*pC*M.V;

Np=size(pC,1);

% Sample from prior (in subspace)
x(:,1) = spm_normrnd(zeros(Np,1),pC,1);
[L(1),L2] = spm_mci_joint (x(:,1),M,U,Y);
Lsum = (beta(2)-beta(1))*L2;
acc(1) = 1;
for j=2:J,
    xs=x(:,j-1); Ls=L(:,j-1);acc(j)=0;
    switch prop,
        case 'mh',
            % Generate sample at next temperature by composing
            % Metropolis moves at different scales
            for s=1:nprop,
                dx=spm_normrnd(zeros(Np,1),scale(s)^2*pC,1);
                xcand=xs+dx;
                [Lcand,L2cand] = spm_mci_joint (xcand,M,U,Y,beta(j));
                dL = Lcand-Ls;
                r = exp(dL);
                alpha = min(1,r);
                test_prob = rand(1);
                if alpha > test_prob
                    % Accept
                    xs = xcand;
                    Ls = Lcand;
                    L2 = L2cand;
                    acc (j) = 1;
                end
            end

        case 'lmc',
            % Generate sample at next temperature using
            % Langevin Monte Carlo
            M.beta=beta(j);
            lgv_mcmc.init=xs;
            lgv_mcmc.maxits=nprop;
            lgv_mcmc.update_obs_noise=0;

            [tmp,stats] = spm_mci_lgv (lgv_mcmc,M,U,Y);

            xs = stats.P(:,end);
            Ls = -stats.E(:,end);
            L2 = stats.L2(end);
            acc (j) = any(stats.acc(2:end));

        otherwise
            disp('Unknown proposal type in spm_mci_ais_single.m');
    end
    x(:,j)=xs;
    L(j)=Ls;
    if j < J
        Lsum = Lsum + (beta(j+1)-beta(j))*L2;
    end
end
logw=Lsum;
P=x(:,J);
E=-L(J);

if mcmc.rec_traj
    nj=size(x,2);
    traj=spm_vec(M.pE)*ones(1,nj)+M.V*x;
else
    traj=[];
end

author: spm


text:
%Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
%Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
%
%THE BSD LICENSE
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions
%are met:
%
%1. Redistributions of source code must retain the above copyright
%   notice, this list of conditions and the following disclaimer.
%2. Redistributions in binary form must reproduce the above copyright
%   notice, this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%
%THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
%IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
%OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
%IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
%INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
%NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
%THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function test_flann
data_path = './';
outcome = {'FAILED!!!!!!!!!', 'PASSED'};

failed = 0;
passed = 0;
cnt = 0;
ok = 1;

    function assert(condition)
        if (~condition)
            ok = 0;
        end
    end

    function run_test(name, test)
        ok = 1;
        cnt = cnt + 1;
        tic;
        fprintf('Test %d: %s...',cnt,name);
        test();
        time = toc;
        if (ok)
            passed = passed + 1;
        else
            failed = failed + 1;
        end
        fprintf('done (%g sec) : %s\n',time,cell2mat(outcome(ok+1)))
    end

    function status
        fprintf('-----------------\n');
        fprintf('Passed: %d/%d\nFailed: %d/%d\n',passed,cnt,failed,cnt);
    end

    dataset = [];
    testset = [];
    function test_load_data
        % load the datasets and testsets
        % use single precision for better memory efficiency
        % store the features one per column because MATLAB
        % uses column major ordering
        % The dataset.dat and testset.dat files can be downloaded from:
        % http://people.cs.ubc.ca/~mariusm/uploads/FLANN/datasets/dataset.dat
        % http://people.cs.ubc.ca/~mariusm/uploads/FLANN/datasets/testset.dat
        dataset = single(load([data_path 'dataset.dat']))';
        testset = single(load([data_path 'testset.dat']))';

        assert(size(dataset,1) == size(testset,1));
    end
    run_test('Load data',@test_load_data);

	match = [];
    dists = [];
    function test_linear_search
        [match,dists] = flann_search(dataset, testset, 10, struct('algorithm','linear'));
        assert(size(match,1) ==10 && size(match,2) == size(testset,2));
    end
    run_test('Linear search',@test_linear_search);

    function test_kdtree_search
        [result, ndists] = flann_search(dataset, testset, 10, struct('algorithm','kdtree',...
                                                          'trees',8,...
                                                          'checks',64));
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('kd-tree search',@test_kdtree_search);

    function test_kmeans_search
        [result, ndists] = flann_search(dataset, testset, 10, struct('algorithm','kmeans',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'checks',120));
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('k-means search',@test_kmeans_search);

    function test_composite_search
        [result, ndists] = flann_search(dataset, testset, 10, struct('algorithm','composite',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'trees', 1,...
                                                          'checks',64));
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('composite search',@test_composite_search);

    function test_autotune_search
        [result, ndists] = flann_search(dataset, testset, 10, struct('algorithm','autotuned',...
                                                          'target_precision',0.95,...
                                                          'build_weight',0.01,...
                                                          'memory_weight',0));
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('search with autotune',@test_autotune_search);

    function test_index_kdtree_search
        [index, search_params ] = flann_build_index(dataset, struct('algorithm','kdtree', 'trees',8,...
                                                          'checks',64));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index kd-tree search',@test_index_kdtree_search);

    function test_index_kmeans_search
        [index, search_params ] = flann_build_index(dataset, struct('algorithm','kmeans',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'checks',120));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index kmeans search',@test_index_kmeans_search);

    function test_index_kmeans_search_gonzales
        [index, search_params ] = flann_build_index(dataset, struct('algorithm','kmeans',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'checks',120,...
                                                          'centers_init','gonzales'));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index kmeans search gonzales',@test_index_kmeans_search_gonzales);

    function test_index_kmeans_search_kmeanspp
        [index, search_params ] = flann_build_index(dataset, struct('algorithm','kmeans',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'checks',120,...
                                                          'centers_init','kmeanspp'));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index kmeans search kmeanspp',@test_index_kmeans_search_kmeanspp);

    function test_index_composite_search
        [index, search_params ] = flann_build_index(dataset,struct('algorithm','composite',...
                                                          'branching',32,...
                                                          'iterations',3,...
                                                          'trees', 1,...
                                                          'checks',64));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index composite search',@test_index_composite_search);

   function test_index_autotune_search
        [index, search_params, speedup ] = flann_build_index(dataset,struct('algorithm','autotuned',...
                                                          'target_precision',0.95,...
                                                          'build_weight',0.01,...
                                                          'memory_weight',0));
        [result, ndists] = flann_search(index, testset, 10, search_params);
        n = size(match,2);
        precision = (n-sum(abs(result(1,:)-match(1,:))>0))/n;
        assert(precision>0.9);
        assert(sum(~(match(1,:)-result(1,:)).*(dists(1,:)-ndists(1,:)))==0);
    end
    run_test('index autotune search',@test_index_autotune_search);

    status();
end

repo: gspbox


text:
function outputImage = edge(this, method, thresh)
% Computes edges of a 3D or 4D image per 2D-slice
%
%   Y = MrImage()
%   outputImage = edge(Y, method, thresh)
%
% This is a method of class MrImage.
%
% IN
%   method      'sobel', 'roberts', 'prewitt' See also edge
%   thresh      custom threshold for edge detection; default: [] determines
%               threshold automatically
% OUT
%   outputImage binary image with detection edges == 1
%
% EXAMPLE
%   edgeY = Y.edge('prewitt', 300);
%
%   See also MrImage perform_unary_operation

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-25
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public 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/>.

if nargin < 2
    method = 'sobel';
end

if nargin < 3
    thresh = [];
end

if isreal(this)
    outputImage = this.perform_unary_operation(...
        @(X) edge(X, method, thresh), '2D');
else % perform on abs for complex data
    outputImage = this.abs.perform_unary_operation(...
        @(X) edge(X, method, thresh), '2D');
end

author: translationalneuromodeling

repo: tapas


text:
classdef KTUEAImpedanceFile

% Plots a graphical representation of the impedance data collected by
% Cerebus and saved in a text file.
%
% Usage:
%
%   plotImpedances(mapfile)
%
%   WHERE:
%
%   mapfile:    Is the mapfile structure obtained with mapfileLoader
%               function that maps the electrodes to channel numbers. If
%               mapfile is not provided then a default map will be used.
%
%   To see a list of methods type methods(plotImpedances).
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.1.0
%
%   Ideas borrowed from Sergey Stavisky, The BrainGate Project

    properties (Hidden, SetAccess = private, GetAccess = private)
        chanCount    = 96;
        impedanceDataValues = NaN(1, 96);
        Mapfile
        fileHandle
        pathHandle
    end
    methods(Hidden)
        function obj = KTUEAImpedanceFile(Mapfile)
            if ~exist('Mapfile', 'var')
                obj.Mapfile = KTUEAMapFile;
            else
                obj.Mapfile = Mapfile;
            end
            if ~obj.Mapfile.isValid; return; end;

            if exist('getFile.m', 'file') == 2
                [obj.fileHandle obj.pathHandle] = getFile('*.txt', 'Open an impedance file...');
            else
                [obj.fileHandle obj.pathHandle] = uigetfile('*.txt', 'Open an impedance file...');
            end
            if ~obj.fileHandle; disp('No file was selected.'); return; end;
            impedanceDataCell = importdata([obj.pathHandle obj.fileHandle], ' ', 200);
            impedanceDataCell(1:9) = [];
            impedanceDataCell(97:end) = [];
            impedanceDataCellParsed = regexp(impedanceDataCell, '\t', 'split');
            for i = 1:size(impedanceDataCellParsed, 1)
                impedanceSingleCell = impedanceDataCellParsed{i,:}(2);
                impedanceSingleText = impedanceSingleCell{:};
                obj.impedanceDataValues(i) = str2double(impedanceSingleText(1:end-4));
            end
        end
    end
    methods
        function FilePath = PathName(obj)
            FilePath = obj.pathHandle;
        end
        function FileName = Filename(obj)
            FileName = obj.fileHandle;
        end
        function validFlag = isValid(obj)
            if any(isnan(obj.getChannelImpedances))
                validFlag = 0;
            else
                validFlag = 1;
            end
        end
        function impedanceDataValues = getImpedances(obj)
            impedanceDataValues = obj.impedanceDataValues;
        end
        function impedanceDataValue = getChannelImpedance(obj, chanNum)
            if ~exist('chanNum', 'var')
                disp('Channel number is a required argument.');
                return;
            end
            impedanceDataValue = obj.impedanceDataValues(chanNum);
        end
        function analyzedStat = getImpedanceMean(obj)
            analyzedStat = mean(obj.impedanceDataValues);
        end
        function analyzedStat = getImpedanceMedian(obj)
            analyzedStat = median(obj.impedanceDataValues);
        end
        function analyzedStat = getImpedanceRange(obj)
            analyzedStat = obj.getImpedanceMax-obj.getImpedanceMin;
        end
        function analyzedStat = getImpedanceSTD(obj)
            analyzedStat = std(obj.impedanceDataValues);
        end
        function analyzedStat = getImpedance975Perc(obj)
            sortedImpedances = sort(obj.impedanceDataValues, 'ascend');
            analyzedStat = sortedImpedances(round(obj.chanCount*0.975));
        end
        function analyzedStat = getImpedance750Perc(obj)
            sortedImpedances = sort(obj.impedanceDataValues, 'ascend');
            analyzedStat = sortedImpedances(round(obj.chanCount*0.75));
        end
        function analyzedStat = getImpedance250Perc(obj)
            sortedImpedances = sort(obj.impedanceDataValues, 'ascend');
            analyzedStat = sortedImpedances(round(obj.chanCount*0.25));
        end
        function analyzedStat = getImpedance25Perc(obj)
            sortedImpedances = sort(obj.impedanceDataValues, 'ascend');
            analyzedStat = sortedImpedances(round(obj.chanCount*0.025));
        end
        function analyzedStat = getImpedanceMin(obj)
            analyzedStat = min(obj.impedanceDataValues);
        end
        function analyzedStat = getImpedanceMax(obj)
            analyzedStat = max(obj.impedanceDataValues);
        end
        function plotStatistics(obj)
            figure('Name', 'Array Statistics');
            title('Array Impedance Values Statistics');
            set(gcf, 'Position', [657 540 288 287]);
            axis off;
            text(0,0.9, 'Mean');                                                                                text(0.6,0.9, num2str(obj.getImpedanceMean, '%5.1f'));
            text(0,0.8, 'Median');            text(0.6,0.8, num2str(obj.getImpedanceMean, '%5.1f'));
            text(0,0.7, 'Range');             text(0.6,0.7, num2str(obj.getImpedanceRange, '%5.1f'));
            text(0,0.6, 'STD');               text(0.6,0.6, num2str(obj.getImpedanceSTD, '%5.1f'));
            text(0,0.5, '97.5th Percentile'); text(0.6,0.5, num2str(obj.getImpedance975Perc, '%5.1f'));
            text(0,0.4, '75th Percentile');   text(0.6,0.4, num2str(obj.getImpedance750Perc, '%5.1f'));
            text(0,0.3, '25th Percentile');   text(0.6,0.3, num2str(obj.getImpedance250Perc, '%5.1f'));
            text(0,0.2, '2.5th Percentile');  text(0.6,0.2, num2str(obj.getImpedance25Perc, '%5.1f'));
            text(0,0.1, 'Min');               text(0.6,0.1, num2str(obj.getImpedanceMin, '%5.1f'));
            text(0,0.0, 'Max');               text(0.6,0.0, num2str(obj.getImpedanceMax, '%5.1f'));
        end
        function plotHistogram(obj, binCount, maxRange)
            if ~exist('binCount', 'var')
                binCount = 20;
            end
            if ~exist('maxRange', 'var')
                if max(obj.impedanceDataValues) > 10000
                    maxRange = 10000;
                else
                    maxRange = max(obj.impedanceDataValues);
                end
            end
            hist(obj.impedanceDataValues, min(obj.impedanceDataValues):...
                                          maxRange/binCount:...
                                          maxRange,1);
        title('Impedance Values Histogram');
        xlabel('Impedance Values (kOhm)');
        ylabel('Impedance Occurances');
        end
        function plotBoxPlot(obj, boxStyle)
            if ~exist('boxStyle', 'var');
                boxStyle = 'traditional';
            end
            sortedImpedances = sort(obj.impedanceDataValues);

            boxplot(obj.impedanceDataValues);
            title(['Box Plot (' boxStyle ')']);
            ylabel('Impedance Values (kOhm)');
            set(gca, 'xTickLabel', ' ');
            xlim([0.9, 1.1]);
            set(gcf, 'Position', [632 356 235 398]);
        end
        function plotImpedances(obj, yelThreshold, redThreshold)
            if ~exist('yelThreshold', 'var')
                yelThreshold = 100;
            end
            if ~exist('redThreshold', 'var')
                redThreshold = 800;
            end
            plotFigure = KTFigure;
            plotFigure.EnlargeFigure;
            plotFigure.MakeBackgroundWhite;
            hold on;
            for channelIDX = 1:obj.chanCount
                if obj.impedanceDataValues(channelIDX) > redThreshold
                    spColor = [1,0,0];
                    fColor  = [1,1,1];
                elseif obj.impedanceDataValues(channelIDX) < redThreshold && obj.impedanceDataValues(channelIDX) > yelThreshold
                    spColor = [0,0,0];
                    fColor  = [1,1,1];
                elseif obj.impedanceDataValues(channelIDX) < yelThreshold
                    spColor = [1,1,0];
                    fColor  = [0,0,0];
                end
                obj.Mapfile.GenerateChannelSubplot(channelIDX);
                obj.Mapfile.GenerateChannelSubplotNames(channelIDX, fColor);
                text(0.05, 0.7, [num2str(obj.impedanceDataValues(channelIDX)) ' kOhm'], 'Color', fColor, 'FontSize', 10, 'FontWeight', 'bold');
                axis on;
                obj.Mapfile.setAxisBackgroundColor(spColor);
                obj.Mapfile.setAxesColor(spColor);
            end
            hold off;
            plotFigure.ShowFigure;
        end
        function plotImpedancesColorBackground(obj)
            for channelIDX = 1:obj.chanCount
                obj.Mapfile.GenerateChannelSubplot(channelIDX);
                axis on;
                obj.Mapfile.setAxisBackgroundColor([0,obj.impedanceDataValues(channelIDX)/max(max(obj.impedanceDataValues)),0]);
            end
        end
        function plotImpedancesComparison(obj)
            changeThreshold = 0.2;
            yelThreshold = 100;
            secImpedanceFile = KTUEAImpedanceFile(obj.Mapfile);
            secImpedanceValues = secImpedanceFile.getChannelImpedances;
            plotFigure = KTFigure;
            plotFigure.EnlargeFigure;
            plotFigure.MakeBackgroundWhite;
            hold on;
            for channelIDX = 1:obj.chanCount
                percentChange = (obj.impedanceDataValues(channelIDX) - secImpedanceValues(channelIDX))/obj.impedanceDataValues(channelIDX);
                if abs(percentChange)  > changeThreshold && percentChange > 0
                    % Lime Green -- If impedance is increased by more
                    % than threshold.
                    if abs(percentChange) > 2 * changeThreshold
                        spColor = [34,139,34] ./255;
                        fColor  = [1,1,1];
                    else
                        spColor = [154,205,50] ./255;
                        fColor  = [0,0,0];
                    end
                elseif abs(percentChange) > changeThreshold && percentChange < 0
                    % Chocolate orange -- If impedance is dropped by more
                    % than threshold.
                    if abs(percentChange) > 2 * changeThreshold
                        spColor = [210,105,30] ./255;
                        fColor  = [1,1,1];
                    else
                        spColor = [245,222,179] ./255;
                        fColor  = [0,0,0];
                    end
                else
                    % Honeydew -- If threshold change is less than
                    % threshold.
                    spColor = [240,255,240] ./255;
                    fColor  = [0,0,0];
                end
                obj.Mapfile.GenerateChannelSubplot(channelIDX);
                patchHandle = patch([0.5, 0.5, 0, 0], [1, 0, 0, 1], spColor);
                obj.Mapfile.GenerateChannelSubplotNames(channelIDX, fColor);
                text(0.05, 0.7, ['1: ' num2str(obj.impedanceDataValues(channelIDX)) ' kOhm'], 'Color', fColor, 'FontSize', 10, 'FontWeight', 'bold');
                text(0.05, 0.5, ['2: ' num2str(secImpedanceValues(channelIDX)) ' kOhm'], 'Color', fColor, 'FontSize', 10, 'FontWeight', 'bold');
            end
            hold off;
            plotFigure.ShowFigure;
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ierror = perm_check2 ( n, p, base )

%*****************************************************************************80
%
%% PERM_CHECK2 checks that a vector represents a permutation.
%
%  Discussion:
%
%    The routine verifies that each of the integers from BASE to
%    to BASE+N-1 occurs among the N entries of the permutation.
%
%    Set the input quantity BASE to 0, if P is a 0-based permutation,
%    or to 1 if P is a 1-based permutation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer P(N), the array to check.
%
%    Input, integer BASE, the index base.
%
%    Output, integer IERROR, error flag.
%    0, the array represents a permutation.
%    nonzero, the array does not represent a permutation.
%
  ierror = 0;

  for seek = base : base + n - 1

    ierror = 1;

    for find = 1 : n
      if ( p(find) == seek )
        ierror = 0;
        break;
      end
    end

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PERM_CHECK2 - Fatal error!\n' );
      fprintf ( 1, '  The input array does not represent\n' );
      fprintf ( 1, '  a proper permutation.\n' );
      error ( 'PERM_CHECK2 - Fatal error!' );
      return
    end

  end

  return
end


text:
function out = isDiag(~)
%ISDIAG   Returns 1 for a SPINOP3.
%
% See also SPINOP3.

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

out = 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function set_label_fontsize(hax, size)

for i=1:length(hax)
  h_xlabel = get(hax(i),'XLabel');
  set(h_xlabel,'FontSize',size);
  h_ylabel = get(hax(i),'YLabel');
  set(h_ylabel,'FontSize',size);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [im_alt] = altertxt(im_tx1,im_tx2)
% This function arrange the pixels alternatively from Text image 1 and 2, and
% returns the interleaved Text image
% im_tx1 - Input Text image 1,3,5 to be encoded
% im_tx2 - Input Text image 2,4,6 to be encoded
% im_alt - Output Interleaved Text image

% Locate Text images 1 and 2 in Odd pixel columns
im_alt(1:2:size(im_tx1,1)-1,1:2:size(im_tx1,2)-1,:) = im_tx1(1:2:size(im_tx1,1)-1,1:2:size(im_tx1,2)-1,:);
im_alt(2:2:size(im_tx1,1)  ,1:2:size(im_tx1,2)-1,:) = im_tx2(2:2:size(im_tx2,1)  ,1:2:size(im_tx2,2)-1,:);
% Locate Text images 1 and 2 in Even pixel columns
im_alt(1:2:size(im_tx1,1)-1,2:2:size(im_tx1,2)  ,:) = im_tx2(1:2:size(im_tx2,1)-1,2:2:size(im_tx2,2)  ,:);
im_alt(2:2:size(im_tx1,1)  ,2:2:size(im_tx1,2)  ,:) = im_tx1(2:2:size(im_tx1,1)  ,2:2:size(im_tx1,2)  ,:);


text:
function [ff,xx,xScatter,yScatter] = BF_ViolinPlot(dataCell,addMeans,doveTail,makeFigure,extraParams)
% Plots a scatter of a set of distributions with data points offset randomly (jittered)
%
%---INPUTS:
% dataCell, a cell where each element is a vector of numbers specifying a distribution
% addMeans, a flag for whether to add a strip for the mean of each group
% doveTail, whether to show a kernel smoothed distributions
% makeFigure, a flag to specify whether to generate a new figure to plot into
% extraParams, a structure of additional custom parameters
%
%---EXAMPLE USAGE:
% dataCell = {randn(1000,1),randn(500,1)*2+1};
% BF_ViolinPlot(dataCell,true,true,true);

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(addMeans)
    addMeans = true; % Add strip for mean of each group
end
if nargin < 3 || isempty(doveTail)
    doveTail = true; % Add kernel distribution
end
if nargin < 4 || isempty(makeFigure)
    makeFigure = true; % Generate a new figure
end
if nargin < 5 || isempty(extraParams)
    extraParams = struct;
end

% ------------------------------------------------------------------------------
% Set extra parameters

numGroups = length(dataCell);

% Custom marker for points within the distribution:
if ~isfield(extraParams,'customSpot')
    customSpot = '.';
else
    customSpot = extraParams.customSpot;
end

% Custom x-offset
if ~isfield(extraParams,'customOffset')
    customOffset = 0;
else
    customOffset = extraParams.customOffset;
end

% Custom offset range; width of each jitter:
% (proportion of each consecutive unit interval to use for the plot)
if ~isfield(extraParams,'offsetRange')
    offsetRange = 0.5;
else
    offsetRange = extraParams.offsetRange;
end

% Horizontal or vertical
if ~isfield(extraParams,'makeHorizontal')
    makeHorizontal = false;
else
    makeHorizontal = extraParams.makeHorizontal;
end

% Custom HitTest settings
if ~isfield(extraParams,'dontHitMe')
    dontHitMe = false;
else
    dontHitMe = extraParams.dontHitMe;
end

% Custom colormap
if ~isfield(extraParams,'theColors')
    if numGroups <= 3
        theColors = BF_GetColorMap('set1',numGroups,1);
    else
        theColors = BF_GetColorMap('accent',numGroups,1);
    end
    if length(theColors) < numGroups
        theColors = arrayfun(@(x) zeros(3,1),1:numGroups,'UniformOutput',false);
    end
else
    theColors = extraParams.theColors;
    % fprintf(1,'Using custom colors\n');
end

%-------------------------------------------------------------------------------
% Reset random number generator for reproducibility:
rng('default');

if makeFigure
    figure('color','w');
end
hold('on');
box('on');

% ------------------------------------------------------------------------------
% Add kernel distribution
% ------------------------------------------------------------------------------
if doveTail
    ff = cell(numGroups,1);
    xx = cell(numGroups,1);
    for i = 1:numGroups
        if isempty(dataCell{i}) || all(isnan(dataCell{i}))
            continue
        end
        if any(isnan(dataCell{i}))
            warning('NaNs in dataCell')
        end
        [f,x] = ksdensity(dataCell{i},'npoints',500);
        f = f/max(f);

        % Only keep range where data exists:
        minKeep = max([1,find(x >= min(dataCell{i}),1,'first')]);
        maxKeep = min([length(x),find(x >= max(dataCell{i}),1,'first')]);
        keepR = minKeep:maxKeep;
        x = x(keepR);
        f = f(keepR);
        if makeHorizontal
            plot(x,customOffset + i + f*offsetRange/2,'-','color',theColors{i},'LineWidth',2)
            plot(x,customOffset + i - f*offsetRange/2,'-','color',theColors{i},'LineWidth',2)

            % Plot top and bottom
            plot(min(x)*ones(2,1),customOffset + i + [-f(1),f(1)]*offsetRange/2,'-','color',theColors{i},'LineWidth',0.1)
            plot(max(x)*ones(2,1),customOffset + i + [-f(end),f(end)]*offsetRange/2,'-','color',theColors{i},'LineWidth',0.1)
        else
            plot(customOffset + i + f*offsetRange/2,x,'-','color',theColors{i},'LineWidth',2)
            plot(customOffset + i - f*offsetRange/2,x,'-','color',theColors{i},'LineWidth',2)

            % Plot top and bottom
            plot(customOffset + i + [-f(1),f(1)]*offsetRange/2,min(x)*ones(2,1),'-','color',theColors{i},'LineWidth',0.1)
            plot(customOffset + i + [-f(end),f(end)]*offsetRange/2,max(x)*ones(2,1),'-','color',theColors{i},'LineWidth',0.1)
        end

        % Keep for dovetailing the jittered scatter points:
        xx{i} = x;
        ff{i} = f;
    end
end

% ------------------------------------------------------------------------------
% Plot jittered scatter for each group
% ------------------------------------------------------------------------------
xScatter = cell(numGroups,1);
yScatter = cell(numGroups,1);
if ~isempty(customSpot)
    pHandles = cell(numGroups,1);
    for i = 1:numGroups
        if doveTail
            xRand = zeros(length(dataCell{i}),1);
            for j = 1:length(dataCell{i})
                try
                    xRand(j) = (rand(1)*offsetRange - offsetRange/2)*ff{i}(find(xx{i} >= dataCell{i}(j),1));
                end
            end
        else
            xRand = rand([length(dataCell{i}),1])*offsetRange - offsetRange/2;
        end
        if makeHorizontal
            xScatter{i} = dataCell{i};
            yScatter{i} = customOffset + i + xRand;
        else
            xScatter{i} = customOffset + i + xRand;
            yScatter{i} = dataCell{i};
        end
        pHandles{i} = plot(xScatter{i},yScatter{i},customSpot,'color',theColors{i});
    end
    if dontHitMe
        for i = 1:length(pHandles)
            pHandles{i}.HitTest = 'off';
        end
    end
end

% ------------------------------------------------------------------------------
% Add strips for means and stds:
% ------------------------------------------------------------------------------
brightenAmount = -0.5;
if any(cellfun(@(x)any(isnan(x)),dataCell)), warning('NaNs in data'); end
for i = 1:numGroups
    try
        brightColor = brighten(theColors{i},brightenAmount);
    catch
        brightColor = theColors{i};
    end

    if addMeans
        if makeHorizontal
            plot(nanmean(dataCell{i})*ones(2,1),[customOffset + i - offsetRange/2,customOffset + i + offsetRange/2],'-',...
                            'color',brightColor,'LineWidth',3)
        else
            plot([customOffset + i - offsetRange/2,customOffset + i + offsetRange/2],nanmean(dataCell{i})*ones(2,1),'-',...
                            'color',brightColor,'LineWidth',3)
        end
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function pathcosts(ch,~,~,controls,~)
F  = controls.F;
ch.add( 1e-3*F^2 );

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function f1 = p01_f1 ( x )

%*****************************************************************************80
%
%% P01_F1 evaluates the first derivative for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = 2.0 * ( x - 2.0 );

  return
end


text:
function c=comp_atrousfilterbank_td(f,g,a,offset)
%COMP_ATROUSFILTERBANK_TD   Uniform filterbank by conv2
%   Usage:  c=comp_atrousfilterbank_fft(f,g,a,skip);
%
%   Input parameters:
%         f   : Input data - L*W array.
%         g   : Filterbank filters - filtLen*M array.
%         a   : Filter upsampling factor - scalar.
%         offset: Delay of the filters - scalar or array of length M.
%
%   Output parameters:
%         c  : L*M*W array of coefficients
%

%input data length
L=size(f,1);
%input channel number
W=size(f,2);
%filter number
M=size(g,2);
g = comp_ups(g,a,1);

%length of filters
filtLen = size(g,1);
skip = -offset;
% Allow filter delay only in the filter support range
if(all(skip>=filtLen) || all(skip<0))
  error('%s: The filter zero index position outside of the filter support.', upper(mfilename));
end

if(numel(skip)==1)
    skip = skip*ones(M,1);
end

% Output memory allocation
c=zeros(L,M,W,assert_classname(f,g));

% Explicitly extend the input. length(fext) = length(f) + 2*(filtLen-1)
fext = comp_extBoundary(f,filtLen-1,'per','dim',1);
% CONV2 does 2-D linear convolution. 'valid' option crops tails
% length(fextconv2) = length(f) + (filtLen-1)
% length(c(:,m,:)) = N
% W channels done simultaneously by conv2
for m=1:M
  c(:,m,:) = comp_downs(conv2(fext,g(:,m),'valid'),1,skip(m),L);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [nQb, pos, vel, q] = ins(w_b, f_b, nQb, pos, vel, gravity, dt)
%% 捷联更新
rotate_vector = w_b*dt;
rotate_vector_norm = norm(rotate_vector);
if(rotate_vector_norm <1e-10) % fix nan issue
    q = [1 0 0 0];
else
    q = [cos(rotate_vector_norm/2); rotate_vector/rotate_vector_norm*sin(rotate_vector_norm/2)]';
end

% 姿态更新
nQb = ch_qmul(nQb, q); %四元数更新（秦永元《惯性导航（第二版）》P260公式9.3.3）
nQb = ch_qnormlz(nQb); %单位化四元数

% 速度更新
f_n = ch_qmulv(nQb, f_b);
dv = (f_n + [0; 0; -gravity]); %比力方程
vel = vel + dv*dt;

% 位置更新
pos = pos + vel*dt;


text:
function [Z] = backProjection(Y, X)
% Back projection technique for fixing the frequency-wise scales of the
% signals estimated by ICA-based blind source separation techniques.
% Both monaural and multichannel outputs are supported.
%
% Coded by D. Kitamura (d-kitamura@ieee.org)
%
% Copyright 2018 Daichi Kitamura
%
% These programs are distributed only for academic research at
% universities and research institutions.
% It is not allowed to use or modify these programs for commercial or
% industrial purpose without our permission.
% When you use or modify these programs and write research articles,
% cite the following references:
%
% # Original paper (The algorithm was called "Rank-1 MNMF" in this paper)
% D. Kitamura, N. Ono, H. Sawada, H. Kameoka, H. Saruwatari, "Determined
% blind source separation unifying independent vector analysis and
% nonnegative matrix factorization," IEEE/ACM Trans. ASLP, vol. 24,
% no. 9, pp. 1626-1641, September 2016.
%
% # Book chapter (The algorithm was renamed as "ILRMA")
% D. Kitamura, N. Ono, H. Sawada, H. Kameoka, H. Saruwatari, "Determined
% blind source separation with independent low-rank matrix analysis,"
% Audio Source Separation. Signals and Communication Technology.,
% S. Makino, Ed. Springer, Cham, pp. 125-155, March 2018.
%
% See also:
% http://d-kitamura.net
% http://d-kitamura.net/demo-ILRMA_en.html
%
% [syntax]
%   [Z, D] = backProjection(Y, X)
%
% [inputs]
%   Y: estimated (separated) signals (frequency bin x time frame x source)
%   X: reference channel of observed (mixture) signal (frequency bin x time frame x 1)
%      or observed multichannel signals (frequency bin x time frame x channels)
%
% [outputs]
%   Z: scale-fitted estimated signals (frequency bin x time frame x source)
%      or scale-fitted estimated source images (frequency bin x time frame x source x channel)
%

% check errors
if (nargin<2)
    error('Too few input arguments.\n');
end
[I,J,M] = size(Y); % frequency bin x time frame x source

% back projection
if (size(X,3)==1) % calculate scale-fixed estimated signals using X(:,:,1)
    A = zeros(1,M,I);
    Z = zeros(I,J,M);
    for i=1:I
        Yi = squeeze(Y(i,:,:)).'; % channels x frames (M x J)
        A(1,:,i) = X(i,:,1)*Yi'/(Yi*Yi');
    end
    A(isnan(A)) = 0;
    A(isinf(A)) = 0;
    for m=1:M
        for i=1:I
            Z(i,:,m) = A(1,m,i)*Y(i,:,m);
        end
    end
elseif (size(X,3)==M) % calculate scale-fixed source images of estimated signals
    A = zeros(M,M,I);
    Z = zeros(I,J,M,M); % frequency bin x time frame x source x channel
    for i=1:I
        for m=1:M
            Yi = squeeze(Y(i,:,:)).'; % channels x frames (M x J)
            A(m,:,i) = X(i,:,m)*Yi'/(Yi*Yi');
        end
    end
    A(isnan(A)) = 0;
    A(isinf(A)) = 0;
    for n=1:M
        for m=1:M
            for i=1:I
                Z(i,:,n,m) = A(m,n,i)*Y(i,:,n);
            end
        end
    end
else
    error('The number of channels in X must be 1 or equal to that in Y.\n');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: d-kitamura

repo: ILRMA

save_path: github-repos/MATLAB/d-kitamura-ILRMA


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function Qw = matVecQw(w,Q,flag)
%
% computes y = Q*w or Q'*w
% where Q =kron(Q3,kron(Q2,Q1)) is efficiently stored as Q= {Q1,Q2,Q3}
% see also splineTransformation2Dsparse
%==============================================================================

function Qw = matVecQw(w,Q,flag)

if nargin == 0,
  help(mfilename);
  return;
end;

% default is Q*w
if ~exist('flag','var'), flag = 'Qw';  end;

if strcmp(flag,'QTw'),
  % build transposes, note: small matrices
  for i=1:length(Q),    Q{i}=Q{i}';  end;
  Qw = matVecQw(w,Q);
  return;
end;

% reconstruct sizes m and p from Q and w, Q{j} is m{j}-by-p{j}

dim = length(Q);
for j=1:length(Q),
  p(j) = size(Q{j},2);
  m(j) = size(Q{j},1);
end;

switch dim,
  case 1,
    Qw = Q{1}*w;
  case 2,
    w  = reshape(w,[p,dim]);
    Qw = zeros(m(1),m(2),dim);
    for i =1:dim,
      Qw(:,:,i) = Q{1}*w(:,:,i)*Q{2}';
    end;
  case 3,
    w  = reshape(w,[p,dim]);
    Qw = zeros(m(1),m(2),m(3),dim);
    for i =1:dim,
        Qw(:,:,:,i) = tensorProd(Q,w(:,:,:,i));
    end;
end;

Qw = reshape(Qw,[],1);
%------------------------------------------------------------------------------

function Qw = tensorProd(Q,w);
for i=1:length(Q),
  p(i,:) = size(Q{i});
end;
Qw = zeros(p(:,1)');

for i1=1:size(Qw,1),
  for i2=1:size(Qw,2),
    for i3=1:size(Qw,3),
      for j1=1:size(w,1),
        for j2=1:size(w,2),
          for j3=1:size(w,3),
            Qw(i1,i2,i3) = Q{1}(i1,j1)*Q{2}(i2,j2)*Q{3}(i3,j3)*w(j1,j2,j3);
          end;
        end;
      end;
    end;
  end;
end;
%==============================================================================


text:
function [I,k] = verifyquad(funfcn,a,b,tol,see,mmax,varargin)
%QUAD         Verified quadrature using Romberg's scheme - rudimentary implementation
%
%The input function must be real and at least 4 times differentiable. The quadrature
%  does NOT work if f or its derivatives have poles in or near the interval [a,b].
%  Not much optimization of the parameters is done; if the routine works well, there
%  is at least a verified inclusion of the integral available.
%
%   [ I , k ] = verifyquad(f,a,b,tol,see,mmax,varargin)
%
%The univariate real function f is integrated from a to b; must at least be 4 times differentiable.
%
% optional input    tol     anticipated (relative) tolerance (default 1e-6)
%                   see     see intermediate results
%                   mmax    maximal order of Romberg scheme <=16; note that derivatives
%                             up to 2*mmax+2 are computed (default mmax=4)
%                   P1,...  extra parameters for function evaluation
%                           f(x,P1,P2,...)
% optional output   k       number of subintervals
%
%A simple example:
%
%   Q = verifyquad('exp(x*sin(x))',0,10)
%
%Routine verifyquad is by no means optimized, just written straightforwardly. But sometimes
%it is even faster (and much more accurate) than the Matlab built-in function quad:
% f = @(x)(sinh(exp(x))), a = 0; b = 5;
%   tic, Approx = quad(f,a,b), toc,
%   tic, Incl = verifyquad(f,a,b), toc
%produces
% f =
%     @(x)(sinh(exp(x)))
% Warning: Maximum function count exceeded; singularity likely.
% > In quad at 100
% Approx =
%   3.3124e+032
% Elapsed time is 0.281967 seconds.
% intval Incl =
%   1.0e+061 *
%     9.6709
% Elapsed time is 0.161995 seconds.
%
%Note that verifyquad is almost twice as fast, and the approximate value is off
%  by several orders of magnitude. But a warning is given.
%It may also happen that no warning is given:
% f = @(x)(sin(x+exp(x))), a = 0; b = 8;
%   tic, Approx = quad(f,a,b), toc,
%   tic, Incl = verifyquad(f,a,b), toc
% f =
%     @(x)(sin(x+exp(x)))
% Approx =
%     0.2511
% Elapsed time is 0.178671 seconds.
% intval Incl =
%     0.3474
% Elapsed time is 1.125084 seconds.
%
%Note, however, that the approximate value is only correct to one figure, and no warning is given.
%
%Based on Romberg's rule with error term by
%   F. L. Bauer, H. Rutishauser, and E. Stiefel: New Aspects in Numerical Quadrature,
%   Proc. Symp. Appl. Math. Vol. XV, AMS 1963, pp. 198-218.
%

% written  05/29/09     S.M. Rump
% modified 11/30/09     S.M. Rump  function check and vectorize
%

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

  % store standard function exception mode
  RealStdFctsExcptnMode = intvalinit('RealStdFctsExcptn',0);
  intvalinit('RealStdFctsExcptnNaN',0);

  % Convert to inline function as needed
  try
    if isa(funfcn,'function_handle')
      strfun = fcnchk(eval(vectorize(funfcn)),length(varargin));
    else
      strfun = fcnchk(vectorize(funfcn),length(varargin));
    end
  catch
    strfun = fcnchk(funfcn,length(varargin));
  end

  if ( nargin<4 ) | isempty(tol)
    tol = 1e-6;                         % default (absolute) tolerance 1e-6
  end
  tol = max(tol,1e-16);                 % avoid tiny tolerances
  D = intval(b)-a;                      % inclusion of b-a
  if ( nargin<5 ) | isempty(see)
    see = 0;
  end
  if ( nargin<6 ) | isempty(mmax)
    mmax = 4;                           % up to 2*mmax+2 derivatives
  end
  if mmax>16
    warning('maximum order of Romberg scheme is mmax=16.')
    mmax = 16;
  end

  % store warning mode
  wng = warning;
  warning off

  I = infsup(-inf,inf);
  relerrold = inf;
  kmin = max(5,mmax);                   % start with at least 2^kmin subintervals

  % bounds for taylor^(2*m+2) for m=1..mmax
  xi = linspace(a,b,17);
  TF_ = feval(strfun,taylorinit(infsup(xi(1:end-1),xi(2:end)),2*mmax+2),varargin{:});
  if any(isnan(TF_.t(:))) | any(isinf(TF_.t(:)))  % second try with many subintervals
    xi = linspace(a,b,1025);
    TF_ = feval(strfun,taylorinit(infsup(xi(1:end-1),xi(2:end)),2*mmax+2),varargin{:});
  end
  if any(isnan(TF_.t(:))) | any(isinf(TF_.t(:)))  % derivative calculation failed
    X = feval(strfun,xi,varargin{:});
    if any(isnan(X(:))) | any(isinf(X(:)))
      error('function evaluation failed over interval [a,b]')
    else
      error('evaluation of higher derivatives failed, try smaller maximum order')
    end
  end
  for m=1:mmax
    TF{m} = infsup(min(TF_{2*m+2}.inf),max(TF_{2*m+2}.sup));
  end
  if see
    disp(['Inclusions of ' int2str(mmax) '-th derivative'])
    TF{mmax}
  end

  % first T{kmin..kmin-mmax+1}(m) for m=1
  k = kmin+1;
  h = 2^(-k)*D;                         % inclusion of (b-a)/2^k
  Xi = a + ( (0:2^k)*(2^(-k)) ) * D;    % inclusion of grid points a+i*(b-a)/2^k
  w = [ 1 repmat([4 2],1,2^(k-1)) ];    % no rounding errors
  w(end) = 1;
  y = feval(strfun,Xi,varargin{:})/3;
  T{k-1}(1) = h * sum(w.*y);
  for j=1:m-1
    w = w(1:(2^(k-j)+1));
    w(end) = 1;
    h = 2*h;                            % inclusion of (b-a)/2^k
    T{k-j-1}(1) = h * sum(w.*y(1:2^j:end));
  end

  % create tableaux up to mmax
  for m=2:mmax
    for k=kmin-mmax+1:kmin-m+1
      T{k}(m) = ( 4^m*T{k+1}(m-1) - T{k}(m-1) ) / ( 4^m-1 );
    end
  end

  % first error term for m=1    [ B(m) = +/- Bernoulli(2*m+2) ]
  B = intval([1 1 1 5 691 7 3617 43867 174611 854513 236364091 8553103 ...
                 23749461029 8615841276005 7709321041217 2577687858367]) ./ ...
             [30 42 30 66 2730 6 510 798 330 138 2730 6 870 14322 510 6];
  m = mmax;
  k = kmin-m+1;
  E = 2^(-2*k*(m+1)-m*(m+1)) * TF{m}*D^(2*m+3);
  err = B(m)*E;
  I = T{k}(m) - err;
  if see
    disp('First inclusion')
    I
  end

  % initialize Romberg iteration
  if in(0,I)
    relerrI = diam(I);
  else
    relerrI = relerr(I);
  end
  j = 0;

  while ( relerrI>tol ) & ( relerrI<relerrold ) & ( kmin+j<20 )
    j = j+1;
    Iold = I;
    relerrold = relerrI;
    k = kmin+j+1;
    % create entry T{kmin+j}(1)
    h = 2^(-k)*D;                         % inclusion of (b-a)/2^k
    Xi = a + ( (0:2^k)*(2^(-k)) ) * D;    % inclusion of grid points a+i*(b-a)/2^k
    w = [ 1 repmat([4 2],1,2^(k-1)) ];    % no rounding errors
    w(end) = 1;
    y = feval(strfun,Xi,varargin{:})/3;
    k = k-1;
    T{k}(1) = h * sum(w.*y);
    % update last row of Romberg table
    for m=2:mmax
      T{k-m+1}(m) = ( 4^m*T{k-m+2}(m-1) - T{k-m+1}(m-1) ) / ( 4^m-1 );
    end
    % new inclusion
    m = mmax;
    k = k-m+1;
    E = 2^(-2*k*(m+1)-m*(m+1)) * TF{m}*D^(2*m+3);
    err = B(m)*E;
    I = T{k}(m) - err;
    if see
      disp('New inclusion')
      I
    end
    if in(0,I)
      relerrI = diam(I);
    else
      relerrI = relerr(I);
    end
  end

  if relerrI>relerrold
    I = Iold;
    k = k-1;
  end

  % restore warning and exception mode
  warning(wng)
  % restore out-of-range exception mode
  intvalinit(RealStdFctsExcptnMode,0);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [o, others] = slover(params, others, varargin)
% class constructor for slice overlay (slover) object
% FORMAT [o, others] = slover(params, others, varargin)
%
% Inputs
% params    - either:
%             - action string implementing class methods (see below)
%             - array of image names / vol structs to display
%             - structure with some fields for object (see below)
% others    - structure, containing extra fields for object (or children)
% varargin  - maybe some other parameters for action calls (see below)
%
% Outputs
% o       - slover object
% others  - any unrecognized fields from params, others
%
% Object fields are:
%  - img - array of structs with information for images to display
%        - img structs contain fields
%             type - one of {'truecolour' 'split', 'contour'};
%                   truecolour - displays transparent (see prop) image
%                      overlaid with any previous
%                   split - in defined area, replaces image present (SPM
%                      type activation display)
%                   contour - contour map of image overlaid.  See help
%                      for contours function in matlab
%             vol - vol struct info (see spm_vol)
%                   can also be vol containing image as 3d matrix
%                   set with add_blobs method
%             cmap - colormap for this image
%             nancol - color for NaN. If scalar, this is an index into
%                    the image cmap.  If 1x3 vector, it's a colour
%             prop - proportion of intensity for this cmap/img
%             func - function to apply to image before scaling to cmap
%                    (and therefore before min/max thresholding. E.g. a func of
%                    'i1(i1==0)=NaN' would convert zeros to NaNs
%             range - 2x1 vector of values for image to distribute colormap across
%                    the first row of the colormap applies to the first
%                    value in 'range', and the last value to the second
%                    value in 'range'
%             outofrange - behavior for image values to the left and
%                    right of image limits in 'range'.  Left means
%                    colormap values < 1, i.e for image values <
%                    range(1), if (range(1)<range(2)), and image values >
%                    range(1) where (range(1)>range(2)). If missing,
%                    display min (for Left) and max (for Right) value from colormap.
%                    Otherwise should be a 2 element cell array, where
%                    the first element is the colour value for image values
%                    left of 'range', and the second is for image values
%                    right of 'range'.  Scalar values for
%                    colour index the colormap, 3x1 vectors are colour
%                    values.  An empty array attracts default settings
%                    appropriate to the mode - i.e. transparent colour (where
%                    img(n).type is truecolour), or split colour.  Empty cells
%                    default to 0. 0 specifies that voxels with this
%                    colour do not influence the image (split =
%                    background, true = black)
%            hold  - resampling order for image (see spm_sample_vol) -
%                    default 1
%            background - value when resampling outside image - default
%                    NaN
%            linespec - string, applies only to contour map,
%                    e.g. 'w-' for white continuous lines
%            contours - vector, applies to contour map only, defines
%                    values in image for which to show contours
%                    (see help contours)
%            linewidth - scalar, width in points of contour lines
%
% - transform - either - 4x4 transformation to apply to image slice position,
%             relative to mm given by slicedef, before display
%               or     - text string, one of axial, coronal, sagittal
%                        These orientations assume the image is currently
%                        (after its mat file has been applied) axially
%                        oriented
% - slicedef - 2x3 array specifying dimensions for slice images in mm
%             where rows are x,and y of slice image, and cols are neg max dim,
%             slice separation and pos max dim
% - slices   - vector of slice positions in mm in z (of transformed image)
% - figure    - figure handle for slice display figure
%               The object used for the display is attached as 'UserData'
%               to this figure
% - figure_struct - stored figure parameters (in case figure dies and
%               needs to be recreated)
% - refreshf  - flag - if set or empty, refresh axis info for figure
%             else assume this is OK
% - clf       - flag, non zero -> clear figure before display.  Redundant
%               if refreshf == 0
% - resurrectf - if not zero, and figure (above) does not exist, will
%               attempt to recreate figure with same area properties.
%               Otherwise painting will give an error.
% - userdata  - flag, non zero -> attaches object to figure when ploting,
%               for use by callbacks (default is 1)
% - area      - struct with fields
%                  position - bottom left, x size y size 1x4 vector of
%                      area in which to display slices
%                  units    - one of
%                    inches,centimeters,normalized,points,{pixels}
%                  halign - one of left,{center},right
%                  valign - one of top,{middle},bottom
% - xslices  - no of slices to display across figure (defaults to an optimum)
% - cbar      - if empty, missing, no colourbar.  If an array of integers, then
%             indexes img array, and makes colourbar for each cmap for
%             that img.  Cbars specified in order of appearance L->R
% - labels - struct can be:
%                  - empty (-> default numerical labels)
%                  - 'none' (string) (no labels)
%                  - or contain fields:
%                  colour - colour for label text
%                  size - font size in units normalized to slice axes
%                  format - if = cell array of strings =
%                  labels for each slice in Z.  If is string, specifies
%                  sprintf format string for labelling in distance of the
%                  origin (Xmm=0, Ymm=0) of each slice from plane containing
%                  the AC, in mm, in the space of the transformed image
% - callback - callback string for button down on image panels. THe
%              following examples assume that you have the 'userdata'
%              field set to 1, giving you access to underlying object
%              To print to the matlab window the equivalent position in
%              mm of the position of a mouse click on one of the image
%              slices, set callback to:
%                   'get_pos(get(gcf, ''UserData''))'
%              To print the intensity values of the images at the clicked point:
%                   ['so_obj = get(gcf, ''UserData''); ' ...
%                    'point_vals(so_obj, get_pos(so_obj))']
% - printstr - string for printing slice overlay figure window, e.g.
%              'print -dpsc -painters -noui' (the default)
% - printfile - name of file to print output to; default 'slices.ps'
%
% Action string formats:
% FORMAT [cmap warnstr] = slover('getcmap', cmapname)
% Gets colormap named in cmapname string
%
% FORMAT [mx mn] = slover('volmaxmin', vol)
% Returns maximum and minimum finite values from vol struct 'vol'
%
% FORMAT vol = slover('blobs2vol', XYZ, vals, mat)
% returns (pseudo) vol struct for 3d blob volume specified
% in matrices as above
%
% FORMAT vol = slover('matrix2vol', mat3d, mat)
% returns (pseudo) vol struct for 3d matrix
% input matrices as above
%
% FORMAT obj = slover('basic_ui' [,dispf])
% Runs basic UI to fetch some parameters, does display, returns object
% If optional dispf parameter = 0, supresses display
%__________________________________________________________________________

% Matthew Brett
% $Id: slover.m 6623 2015-12-03 18:38:08Z guillaume $

myclass = 'slover';

% Default object structure
defstruct = struct('img', [], ...
    'transform', 'axial', ...
    'slicedef', [], ...
    'slices', [], ...
    'figure', [], ...
    'figure_struct', [], ...
    'refreshf', 1, ...
    'clf', 1, ...
    'resurrectf', 1, ...
    'userdata', 1, ...
    'area', [], ...
    'xslices', [], ...
    'cbar', [], ...
    'labels', [], ...
    'callback', ';', ...
    'printstr', 'print -dpsc -painters -noui', ...
    'printfile', 'slices.ps');

if nargin < 1
    o = class(defstruct, myclass);
    others = [];
    return
end
if nargin < 2
    others = [];
end

% parse out string action calls (class functions)
if ischar(params)
    switch params
        case 'getcmap'
            if nargin < 2
                error('Need colormap name');
            end
            o = pr_getcmap(others);
            return
        case 'volmaxmin'
            if nargin < 2
                error('Need volume to calculate max/min');
            end
            [o,others] = pr_volmaxmin(others);
            return
        case 'blobs2vol'
            if nargin < 4
                error('Need XYZ, vals, mat');
            end
            o = pr_blobs2vol(others, varargin{:});
            return
        case 'matrix2vol'
            if nargin < 3
                error('Need matrix and mat');
            end
            o = pr_matrix2vol(others, varargin{:});
            return
        case 'basic_ui'
            o = pr_basic_ui(others, varargin{:});
            if ~isempty(o), o = paint(o); end
            return
    end

    % if not action string, must be filename(s)
    params = spm_vol(params);
end

% Could these just be image vol structs?
if isfield(params, 'fname')
    for i = 1:numel(params)
        obj.img(i).vol = params(i);
    end
    params = obj;
end

% Deal with passed objects of this (or child) class
if isa(params, myclass)
    o = params;
    % Check for simple form of call
    if isempty(others), return, end

    % Otherwise, we are being asked to set fields of object
    [p,others] = mars_struct('split', others, defstruct);
    o = mars_struct('ffillmerge', o, p);
    return
end

% fill params with defaults, parse into fields for this object, children
params = mars_struct('fillafromb', params, others);
[params, others] = mars_struct('ffillsplit', defstruct, params);

% set the slover object
o  = class(params, myclass);

% refill with defaults
o = fill_defaults(o);


text:
% EEG_EVENTTABLE - returns all events contained in the EEG structure (and
%                    optionally exports them to a CSV file)
%
% Usage:
%   >> events = eeg_eventtable(EEG);
%
% Inputs:
%   EEG          - EEG structure
%
% Optional inputs:
%   'unit'       - Unit of events containing time samples, can be either
%                  'samples' (default) or 'seconds'
%   'dispTable'  - Display an overview of all events if set to true (default);
%                  do not display if set to false
%   'exportFile' - Exports events as a CSV file (using tabs as delimiters); set
%                  this parameter to the file name you wish to export the events
%                  to (string); default: no file is exported
%   'indexCol'   - Adds a column with the event indices as the first row; note
%                  that this is not an event field; default: true
%
% Outputs:
%   events       - event structure (cell array)
%
% Examples:
%   Get all events contained in the EEG structure and display an overview:
%     >> events = eeg_eventtable(EEG);
%
%   In addition to displaying an overview, export into a CSV file:
%     >> events = eeg_eventtable(EEG, 'exportFile', 'test.csv');

% Copyright by Clemens Brunner <clbrunner@ucsd.edu>, 2011
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% Revision: 0.15
% Date: 10/13/2011

% Revision history:
%   0.15: Added option to add/remove index column
%   0.11: Changed function name and updated documentation
%   0.10: Initial version

function events = eeg_eventtable(EEG, varargin)

% Default parameters, can be overwritten by varargin
unit = 'samples';  % Unit of latencies is 'samples'
dispTable = true;  % Display the event table
exportFile = [];  % Export events to text file
indexCol = true;  % Add an index column as first column

if ~isempty(varargin)  % Are there optional parameters available?
    k = 1;
    while k <= length(varargin)
        if strcmp(varargin{k}, 'dispTable')
            dispTable = varargin{k + 1};
            k = k + 2;
        elseif strcmp(varargin{k}, 'unit')
            unit = varargin{k + 1};
            k = k + 2;
        elseif strcmp(varargin{k}, 'exportFile')
            exportFile = varargin{k + 1};
            k = k + 2;
        elseif strcmp(varargin{k}, 'indexCol')
            indexCol = varargin{k + 1};
            k = k + 2;
        else  % Ignore unknown parameters
            k = k + 2;
        end
    end
end

fields = fieldnames(EEG.event);  % Get all event fields
nFields = length(fields);  % Number of event fields
nEvents = length(EEG.event);  % Number of events

events = cell(nFields, 1);
for k = 1:nFields
    events{k}.name = fields{k};  % Event name
    events{k}.values = cell(nEvents, 1);  % Will contain values for this event

    for l = 1:nEvents  % Write values
        if strcmp(events{k}.name, 'latency') && strcmp(unit, 'seconds')
            events{k}.values{l} = getfield(EEG.event(l), fields{k})/EEG.srate;
        elseif strcmp(events{k}.name, 'duration') && strcmp(unit, 'seconds')
            events{k}.values{l} = getfield(EEG.event(l), fields{k})/EEG.srate;
        else
            events{k}.values{l} = getfield(EEG.event(l), fields{k});
        end
    end

end

if dispTable || ~isempty(exportFile)  % Display overview or export to CSV file?

    nr = cell(nEvents, 1);  % Consecutive event numbers
    for l = 1:nEvents
        nr{l} = l;  end

    eventsTable = [];
    if indexCol  % Add index column
        titleTable{1} = 'number';
        titleOffset = 1;
    else  % Do not add index column
        titleOffset = 0;
    end
    for k = 1:nFields
        eventsTable = [eventsTable, events{k}.values];
        titleTable{k + titleOffset} = events{k}.name;
    end
    if indexCol
        eventsTable = [nr, eventsTable];  % Add event numbers in first column
    end
    finalTable = [titleTable; eventsTable];  % Add title in first row

    if dispTable
        format('shortG');
        disp(finalTable);
        format;
    end

    if ~isempty(exportFile)
        fid = fopen(exportFile, 'w');
        for rows = 1:size(finalTable, 1)
            for cols = 1:size(finalTable, 2)
                dataValue = finalTable(rows, cols);
                if cols < size(finalTable, 2)  % If we're not in the last column, print the delimiter
                    if iscellstr(dataValue)
                        fprintf(fid, '%s\t', dataValue{:});
                    elseif iscellnumeric(dataValue)
                        fprintf(fid, '%f\t', dataValue{:});
                    end
                else  % If we're in the last column, do not print the delimiter
                    if iscellstr(dataValue)
                        fprintf(fid, '%s', dataValue{:});
                    elseif iscellnumeric(dataValue)
                        fprintf(fid, '%f', dataValue{:});
                    end
                end
            end
            fprintf(fid, '\n');
        end
        fclose(fid);
    end
end

function b = iscellnumeric(C)
% Return 1 if all elements of cell array are numeric

b = all(cellfun(@(x) isnumeric(x),C));

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [mi] = mm2mi(mm)
% Convert length from millimeters to miles.
% Chad A. Greene 2012
mi = mm*6.213711922373e-7;

author: Sable

repo: mcbench-benchmarks