text:
function obj = export(obj, varargin)
% export Export gramm plot to vector or bitmap image
%
% Example syntax: gramm_object.export('file_name','my_figure','export_path','./test/','file_type','pdf')
% export() must be called after a draw() call and will save the plot to
% file using the specified 'name',value options:
% - 'file_name' Desired name for the file. Extension will be added
%   according to 'file_type'
% - 'export_path' Path of the destination folder (absolute or relative).
    %By default the file is saved in the current folder
% - 'file_type': supports 'svg','pdf' and 'eps' for vector output, and
%   'png' and 'jpg' for high-resolution bitmap output. The default is 'svg' on 2014b+ Matlab
%   versions, and 'pdf' on older versions.
% - 'width' and 'height' allow to set the desired size of the exported
%   figure. By default (recommended), the exported figure will look the same as the on-screen
%   plot. Setting this manually is tricky.
% - 'units' allows to set the units used for the 'width' and 'height'
%   input. Supports 'centimeters' and 'inches'. Default is 'centimeters'.

% Retrieve figure handle (all are the same so only the first is ok)
h_fig = obj(1).parent;

% calculating defaults
fig_pos = getpixelposition(h_fig);
%Resolution in pixels per cm (note: not the actual resolution, Matlab doesn't
%know about real screen size).
resolution = get(0, 'ScreenPixelsPerInch')/2.54;
width_fig = fig_pos(3)/resolution;
height_fig = fig_pos(4)/resolution;

%SVG is unsupported on older Matlab versions
if obj(1).handle_graphics
    default_file_type='svg';
else
    default_file_type='pdf';
end

% parse arguments and set defaults
p=inputParser;
my_addParameter(p,'file_name','gramm_export');
my_addParameter(p,'export_path','');
my_addParameter(p,'file_type', default_file_type);
my_addParameter(p,'width', width_fig);
my_addParameter(p,'height', height_fig);
my_addParameter(p,'units', 'centimeters');
parse(p,varargin{:});

% Turn off warning due to a custom defined figure resize function called on save
warning('off', 'MATLAB:print:CustomResizeFcnInPrint');

%Verify that the path is correct
if ~isempty(p.Results.export_path) && ~exist(p.Results.export_path,'dir')
    error('The specified folder does not exist')
end

% Building variables from (default) input data
file_path = fullfile(p.Results.export_path, p.Results.file_name);
width = p.Results.width;
height = p.Results.height;

% Verify that parent object is a figure object
if obj(1).handle_graphics
    if ~isa(h_fig, 'matlab.ui.Figure')
        error('Parent gramm object is not a matlab.ui.Figure.');
    end
end

% In case units are inches, we rescale to cm
if strcmp(p.Results.units, 'inches')
	if ~any(strcmp(p.UsingDefaults, 'width'))
		width = p.Results.width * 2.54;
	end

	if ~any(strcmp(p.UsingDefaults, 'height'))
		height = p.Results.height * 2.54;
	end
end

% Cropping paper to fit the figure
set(h_fig, 'Units', 'centimeters');
set(h_fig, 'PaperUnits', 'centimeters');
set(h_fig, 'PaperSize', [width, height]);
set(h_fig, 'PaperPosition', [0, 0, width, height]);

% Save figure in preferred file type
switch p.Results.file_type
	case 'pdf'
		print(h_fig, file_path, '-dpdf', '-painters');

	case 'eps'
		print(h_fig, file_path, '-depsc', '-painters', '-r300');

	case 'svg'
		print(h_fig, file_path, '-dsvg', '-painters');

	case 'png'
		print(h_fig, file_path, '-dpng', '-opengl', '-r300');

	case 'jpg'
		warning('JPEG is not recommanded for saving figures. Use PDF or SVG whenever possible');
		print(h_fig, file_path, '-djpeg', '-opengl', '-r300');

	otherwise
		error('Argument file_type not recognized. Available options are: pdf, eps, svg, png, jpg.')
end

% Turn warning back on
warning('on', 'MATLAB:print:CustomResizeFcnInPrint');

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function T_hat = run_tc(params)
  %T_hat = run_tc_full(params);
  T_hat = run_tc_blk(params);
end

author: andrewssobral


text:
function optionsDyn = vargplvmOptionsDyn(optionsDyn, X)
% VARGPLVMOPTIONSDYN Fill in an options structure with default parameters
% FORMAT
% DESC takes an VARGPLVM options structure amends it with default values
% wherever the field is empty.
% ARG optionsDyn : the VARGPLVM structure to fill in with default values.
% ARG X : a latent space initialised for the basic model (only used if
% times t are not given)
%
% COPYRIGHT : Michalis K. Titsias, 2011
% COPYRIGHT : Neil D. Lawrence, 2011
% COPYRIGHT : Andreas C. Damianou, 2011
%
% SEEALSO : vargplvmInitDynamics, vargplvmAddDynamics

% VARGPLVM

if isfield(optionsDyn, 'type') && ~isempty(optionsDyn.type)
    optionsDyn.type = 'vargpTime';
end

if strcmp(optionsDyn.type, 'vargpTime')
    if ~isfield(optionsDyn, 'initX')
        optionsDyn.initX = 'ppca';
    end

    if  ~isfield(optionsDyn, 'constrainType') || isempty(optionsDyn.constrainType)
        % If the dynamics kernel is of type invcmpnd, this cell will
        % explain what type of constraints each element will impose (if
        % other kernel is used, then this cell will only have one element).
        % Possible values (so far) include: "time" and "data".
        optionsDyn.constrainType = {'time'};
    end

    % Create time vector for each dimension; if t is not given, assume an
    % equally spaced time vector.
    if  ~isfield(optionsDyn, 't') || isempty(optionsDyn.t)
        fprintf(1, '# Time vector unknown; creating random, equally spaced time vector\n');
        t = linspace(0, 2*pi, size(X, 1)+1)';
        t = t(1:end-1, 1);
        optionsDyn.t = t;
    else
        t=optionsDyn.t;
    end

    if ~isfield(optionsDyn, 'kernType') || isempty(optionsDyn.kernType)
        optionsDyn.kernType = {'rbf','white', 'bias'};
    end

    % A better way to initialize the  kernel hyperparameter,
    % especially lengthscale, should be to learn it by running few iterations of
    % GP regression marginal likelihood maximization given as data the PCA output
    % (the default value is jsut reasonable and it sets the inverse lenthscale to quite
    % small value so the GP dynamic prior is weaker (less smoothed) at
    % initializations
    if  ~isfield(optionsDyn, 'kern') || isempty(optionsDyn.kern)
       % kern = kernCreate(t, {'rbf','white'});
        kern = kernCreate(t, optionsDyn.kernType);
        kern.comp{2}.variance = 1e-3; % 1e-1

        if  isfield(optionsDyn, 'inverseWidth')
            invWidth=optionsDyn.inverseWidth;
        else
            invWidth=5;
        end

        % The following is related to the expected number of zero-crossings.
        if kern.inputDimension == 1
            kern.comp{1}.inverseWidth = invWidth./(((max(t)-min(t))).^2);
        end
        kern.comp{1}.variance = 1;

        optionsDyn.kern = kern;
    end

    if ~isfield(optionsDyn,'seq')
        optionsDyn.seq=[];
    end

    % Set to 1 to reoptimise all inducing points during test time
     if ~isfield(optionsDyn, 'testReoptimise')
     	optionsDyn.testReoptimise = 1;  % !!! SET THE DEFAULT TO 1
     end

     % Set to 0 so that if you use a kernel from the exp. family for the
     % dynamics, then its variance is not learned and the lenghtscale
     % compensates this (we have one less free parameter then).
    if ~isfield(optionsDyn, 'learnVariance')
        optionsDyn.learnVariance = 0; % !!! SET THE DEFAULT TO 0
    end

    % If set to 1 (not recommended) then the initial means have some
    % smoothing in the first and last values (i.e. closer to zero). Useful
    % mostly for periodic kernels.
    if ~isfield(optionsDyn, 'regularizeMeans')
        optionsDyn.regularizeMeans = 0;
    end

end

author: SheffieldML

repo: vargplvm


text:
% The COBRAToolbox: testSingleProductionEnvelope.m
%
% Purpose:
%     - test the singleProductionEnvelope function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

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

% test variables
model = getDistributedModel('ecoli_core_model.mat');
deletions = model.rxns(21); % 'EX_acald(e)'
deletions2 = model.genes(21);
product = char(model.rxns(20)); % 'EX_ac(e)' - the created and deleted plot has the name of the reaction
biomassRxn = char(model.rxns(22));
refData_y1 = zeros(20,1);
refData_y2 = 20 * ones(20, 1);

% function outputs
% requires Global Optimization Toolbox
% this function's output is a curve, in tested case one of them y=0, the other y=20-2x
[x, y1, y2] = singleProductionEnvelope(model, deletions, product, biomassRxn, 'savePlot', 1);
singleProductionEnvelope(model, deletions2, product, biomassRxn, 'geneDelFlag', 1);

% tests
assert(isequal(refData_y1, y1));
assert(isequal(((refData_y2 - (2 * x') - y2) < 1e-4), ones(20, 1)));

% remove the created plots and the folder - .gitignore for windows because apparently 'rmdir' fails sometimes
if isunix
try
  rmdir([CBTDIR filesep 'tutorials' filesep 'additionalTutorials' filesep 'optKnock' filesep 'Results'], 's');
end
end

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function varargout = displayBVPinfo(handles, mode, varargin)
%DISPLAYBVPINFO   Show information on the CHEBGUI figure when solving BVPs.
%
% Calling sequence:
%   VARARGOUT = DISPLAYBVPINFO(HANDLES, MODE, VARARGIN)
% where
%   VARARGOUT when MODE = 'INIT':
%       VARARGOUT{1}:   Dummy output, required for consistency with the
%                       chebop/displayInfoInit method.
%       VARARGOUT{2}:   A tic/toc timer.
%   VARARGOUT when MODE = 'ITER':
%       VARARGOUT{1}:   A tic/toc timer.
%       VARARGOUT{2}:   Boolean variable, equal to 1 if the user pressed the
%                       'Stop' button on the GUI.
%   HANDLES:            The MATLAB handle to the CHEBGUI figure.
%   VARARGIN:           Useful input arguments for showing information, further
%                       described in 'help chebop/displayInfo'.
%
% See also: chebop/displayInfo

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

% Go into different modes, depending on where in the Newton iteration we are.
switch mode
    case 'exactInitial'
        displayBVPinfoExactInitial(handles, varargin{:});
    case 'init'
        [displayFig, displayTimer] = displayBVPinfoInit(handles, varargin{:});
        varargout{1} = displayFig;
        varargout{2} = displayTimer;
    case 'iter'
        [varargout{1}, varargout{2}] = displayBVPinfoIter(handles, varargin{:});
    case 'final'
        displayBVPInfoFinal(handles, varargin{:});
    case 'linear'
        displayBVPinfoLinear(handles, varargin{:});
end

end

function displayBVPinfoExactInitial(handles, pref)
%DISPLAYINFOEXACTINITIAL    Called if we pass CHEBOP/SOLVEBVP an initial guess
%                           that solves the BVP.

% Show that we got passed an initial guess that solves the BVP:
str = {'Initial guess appears to be a solution of the BVP.';
    'Returning initial guess as the solution.'};

% Update the information on the GUI:
set(handles.iter_list, 'String',  str);
set(handles.iter_list, 'Value', 1);

end

function [displayFig, displayTimer] = displayBVPinfoInit(handles, u, pref)
%DISPLAYBVPINFOINIT   Show information at the start of Newton iteration
%
% See also: chebop/displayInfoInit

% Update the iteration information header on the GUI
initString = [' Iter.       || du ||        Contraction      stepsize     ',...
    '   len(du)          len(u)'];
set(handles.iter_text, 'String', initString);

% Clear the update plot
axes(handles.fig_norm)
plot(0, 0)
xlim(handles.xLim)
cla

% Switch focus to the fig_sol plot on the GUI
axes(handles.fig_sol)

% Plot initial guess
plot(chebfun(u), '.-', 'linewidth', 2)
set(handles.panel_figSol, 'title', 'Initial guess of solution')

% Remove title from bottom plot panel
set(handles.panel_figNorm, 'title', 'Current correction step')

% Update the fontsize of plots
set(handles.fig_sol, 'fontsize', handles.fontsizePanels);
set(handles.fig_norm, 'fontsize', handles.fontsizePanels);

% Do different things for the axes depending on if the solution is real.
if ( isreal(chebfun(u)) )
    xlim(handles.xLim)
else
    axis equal
end

drawnow

% This only needs to be a dummy output, since the HANDLES variable contains
% pointers to the plots we want to draw on:
displayFig = 42;

% Start a timer, used to control pausing between plotting:
displayTimer = tic;

end

function [displayTimer, stopReq] = displayBVPinfoIter(handles, u, delta, ...
    iterNo, normDelta, cFactor, lenDelta, lambda, lenu, displayFig, ...
    displayTimer, pref) %#ok<INUSD>
%DISPLAYBVPINFOITER   Show information and plot during Newton iteration.
%
% See also: chebop/displayInfoIter

% Create a string for displaying information about the iteration.
if ( lambda == 1 )
    iterString = sprintf(' %2.2d %10.2e %10.2e %9.4f %7i %9i ', ...
        iterNo, normDelta, cFactor, lambda, lenDelta, lenu);
else
    iterString = sprintf(' %2.2d %10.2e %10.2e %10.2e %6i %9i ', ...
        iterNo, normDelta, cFactor, lambda, lenDelta, lenu);
end

% Update the iteration information on the GUI.
currString = get(handles.iter_list, 'String');
set(handles.iter_list, 'String', [currString ; iterString]);
set(handles.iter_list, 'Value', iterNo);

% Switch focus to the fig_sol plot on the GUI.
axes(handles.fig_sol)

% Plot current iterate
plot(chebfun(u), '.-', 'linewidth', 2)
set(handles.panel_figSol, 'title', 'Current solution')

% Do different things for the axes depending on if the solution is real.
if ( isreal(chebfun(u)) )
    xlim(handles.xLim)
else
    axis equal
end

% Now plot the Newton updates.

% Switch focus to the fig_norm plot on the GUI.
axes(handles.fig_norm)

% Plot Newton update
plot(chebfun(delta), '.-',  'linewidth', 2)
set(handles.panel_figNorm, 'title', 'Current correction step')

% Update the fontsize of plots
set(handles.fig_sol, 'fontsize', handles.fontsizePanels);
set(handles.fig_norm, 'fontsize', handles.fontsizePanels);

% Do different things for the axes depending on if the solution is real.
if ( isreal(chebfun(u)) )
    xlim(handles.xLim)
else
    axis equal
end

drawnow

% If the user has pressed the pause button on the GUI, we pause.
if ( ~isempty(handles) ...
            && strcmpi(get(handles.button_clear, 'String'), 'Continue') )
    waitfor(handles.button_clear, 'String')
end

% If the user has pressed the stop button on the GUI, we stop the Newton
% iteration.
if ( ~isempty(handles) ...
            && strcmpi(get(handles.button_solve, 'String'), 'Solve') )
    stopReq = true;
else
    stopReq = false;
end

% Restart the timer:
displayTimer = tic;

end

function displayBVPInfoFinal(handles, u, delta, iterNo, errEstDE, errEstBC, ...
	displayFig, displayTimer, pref) %#ok<INUSL>
%DISPLAYBVPINFOFINAL   Utility routine for displaying nonlinear solve progress.
%
%  This method prints out information after Newton iteration finishes when
%  problems are solved using CHEBGUI.
%
% See also: displayBVPInfo, chebop/displayInfoFinal

% How many steps did we need
if ( iterNo == 1 )
    finalStr = {'Newton''s method converged in 1 iteration\n'};
else
    finalStr = {sprintf('Newton''s method converged in %i iterations.',....
        iterNo)};
end

% Show what discretization was used
if ( strcmpi(func2str(pref.discretization), 'coeffs') || ...
        isequal(pref.discretization, @ultraS) || ...
        isequal(pref.discretization, @trigspec))
    discString = 'Coefficients';
else
    discString = 'Values';
end
finalStr = [finalStr; sprintf('Discretization basis used: %s.', ...
    discString)];

% Print info about the final error estimates...
finalStr = [finalStr; ...
    sprintf('Final error estimate: %.2e (differential equation) ', errEstDE)];
% ...and the error in the boundary conditions.
finalStr = [finalStr; ...
    sprintf('%30.2e (boundary conditions).', errEstBC)];

% Update the iteration information on the GUI
currString = get(handles.iter_list, 'String');
set(handles.iter_list,'String', [currString ; finalStr]);

% Set focus to bottom of printed list
set(handles.iter_list, 'Value', iterNo + 4);

end

function displayBVPinfoLinear(handles, u, nrmRes, pref)
%DISPLAYBVPINFOLINEAR   Display methods tailored to linear BVPs.

% Show that we detetected a linear equation:
str = {'Linear equation detected. Converged in one step.'};

% Show what discretization was used:
if ( strcmpi(func2str(pref.discretization), 'coeffs') || ...
        isequal(pref.discretization, @ultraS) || ...
        isequal(pref.discretization, @trigspec))
    discString = 'Coefficients';
else
    discString = 'Values';
end

% Concatenate strings:
str = [str ; sprintf('Discretization basis used: %s.',  discString)];
str = [str ; sprintf('Length of solution: %i.', length(chebfun(u)))];
str = [str ; sprintf('Norm of residual: %.2e.', nrmRes)];

% Update the information on the GUI:
set(handles.iter_list, 'String',  str);
set(handles.iter_list, 'Value', 1);

% Plot
axes(handles.fig_sol)
plot(u, 'linewidth', 2)
% Do different things depending on whether the solution is real or not
if ( isreal(u{1}) )
    axis tight
else
    axis equal
end

set(handles.panel_figSol, 'title', 'Solution')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [pfdata,pfunits,update] = readPFData(pfname,flies,params)

ss = fopen(pfname,'r');
qq = fread(ss,10,'*char')';
fclose(ss);
if nargin < 3
    params = struct;
    do_check = false;
else
    do_check = true;
end

% Doing matfile is slower for more than 4 flies
saved_params = struct;
update = false;
if strcmp(qq(8:10),'7.3') && numel(flies)<4
  m = matfile(pfname);
  pfdata = cell(1,numel(flies));
  for ndx = 1:numel(flies);
    pfdata(ndx) = m.data(1,flies(ndx));
  end
  pfunits = m.units;
  if isprop(m,'params')
      saved_params = m.params;
  else
      do_check = false;
  end
else
  m = load(pfname);
  pfdata = m.data(flies);
  pfunits = m.units;
  if isfield(m,'params')
      saved_params = m.params;
  else
      do_check = false;
  end
end

if ~isequal(saved_params,params) && do_check ...
    && ~isempty(saved_params) && ~isempty(params)
    update = true;
end

author: kristinbranson

repo: JAABA


text:
function [Q_mean]=imlabelMean(M,ML)

% function [Q_mean]=imlabelMean(M,ML)
% ------------------------------------------------------------------------
%
% This function takes the mean for each of the labeled groups (NaN's
% ignored) in ML according to the intensities in M.
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/04/01
%------------------------------------------------------------------------

defaultMethod=1; %Hard coded for now
switch defaultMethod
    case 1 %SPARSE ARRAY BASED (more efficient for large arrays)
        labelSet=unique(ML(:));
        labelSet=labelSet(~isnan(labelSet)); %Remove the nan labelled group

        numLabels=max(labelSet);

        logic_ML_not_nan=~isnan(ML);
        nnzQ=nnz(logic_ML_not_nan);
        Iq=ML(logic_ML_not_nan);
        Jq=1:nnzQ;
        Sq=M(logic_ML_not_nan);
        sizQ=[numLabels,nnzQ];
        nnzQ=nnz(logic_ML_not_nan);
        Q=sparse(Iq,Jq,Sq,sizQ(1),sizQ(2),nnzQ);

        Q_sum=full(sum(Q,2));
        Q_voxelCount=full(sum(spones(Q),2));
        Q_mean=Q_sum./Q_voxelCount;
    case 2 %REGIONPROPS BASED
        ML(isnan(ML))=0;
        A=regionprops(ML,M,'MeanIntensity');
        Q_mean=[A.MeanIntensity];
        Q_mean=Q_mean(:);
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [xi, gamma, Phi] = mfldainit (W, K, w)

  % Get the number of documents.
  D = length(w);

  % Get the size of the corpus.
  M = 0;
  for d = 1:D
    M = M + length(w{d});
  end

  % Generate a random starting point for xi.
  xi = rand(W,K);

  % Generate a random starting point for gamma.
  gamma = rand(K,D);

  % Generate a random starting point for phi.
  Phi = rand(K,M);
  Phi = Phi ./ repmat(sum(Phi),K,1);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [err, Info] = Info_1D(v, fname)
%
% [err, Info] = Info_1D(M, name);
% an Info structure for a 1D file. To use with WriteBrik
%
% Specify name if you want Info.RootName and Info.Extension_1D to be filled out
% otherwise, you can do that on your own later.
%
% see also Read_1D, BrikLoad, and BrikInfo
%     Ziad S. Saad saadz@mail.nih.gov, SSCC/NIMH/NIH/USA

FuncName = 'Info_1D';
err = 1;
Info = [];

if (isempty(v)),
   fprintf(2,'Error %s:\nEmpty M.\n', FuncName);
   return;
end

%some fake Info stuff
if (nargin == 2),
   [Info.RootName, Info.Extension_1D] = Remove1DExtension(fname);
else
   fname = '';
   Info.RootName = '';
   Info.Extension_1D = '';
end

Info.TypeName = '';
Info.TypeBytes = 0;
Info.Orientation = '';
Info.ByteOrder = '';
Info.FileFormat = '1D';
Info.DATASET_DIMENSIONS = [size(v, 1) 1 1 0 0];
Info.DATASET_RANK = [3 size(v,2) 0 0 0 0 0 0];
Info.BRICK_TYPES = [];
Info.BRICK_STATS = [];
Info.BRICK_FLOAT_FACS = '';
Info.BYTEORDER_STRING = '';
Info.ORIENT_SPECIFIC = [];
Info.ORIGIN = [0.0 0.0 0.0];
Info.DELTA = [1 1 1];
Info.BRICK_LABS = '';
Info.BRICK_KEYWORDS = '';
Info.SCENE_DATA = [];
Info.TYPESTRING = '';
Info.IDCODE_STRING = '';
Info.IDCODE_DATE = '';
Info.BRICK_STATAUX = [];
Info.STAT_AUX = [];
Info.HISTORY_NOTE = [];
Info.IDCODE_ANAT_PARENT = '';

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p = prior_sqinvunif(varargin)
%PRIOR_SQINVUNIF  Uniform prior structure for the square inverse of the parameter
%
%  Description
%    P = PRIOR_SQINVUNIF creates uniform prior structure for the
%    square inverse of the parameter.
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010,2012 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQINVUNIFORM';
  ip.addOptional('p', [], @isstruct);
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'SqInv-Uniform';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'SqInv-Uniform')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqinvunif_pak;
    p.fh.unpak = @prior_sqinvunif_unpak;
    p.fh.lp = @prior_sqinvunif_lp;
    p.fh.lpg = @prior_sqinvunif_lpg;
    p.fh.recappend = @prior_sqinvunif_recappend;
  end

end

function [w, s, h] = prior_sqinvunif_pak(p, w)
  w=[];
  s={};
  h=[];
end

function [p, w] = prior_sqinvunif_unpak(p, w)
  w = w;
  p = p;
end

function lp = prior_sqinvunif_lp(x, p)
  lJ = -log(x)*3 + log(2);  % log(-2/x^3) log(|J|) of transformation
  lp = sum(lJ);
end

function lpg = prior_sqinvunif_lpg(x, p)
  lJg = -3./x;              % gradient of log(|J|) of transformation
  lpg = lJg;
end

function rec = prior_sqinvunif_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function cellInfo = GCS_Cmd_Bus(varargin)
% GCS_CMD_BUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'GCS_Cmd_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'mode', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'cmd_1', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Operation channel 1')}; ...
{'cmd_2', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Operation channel 2')}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function gf=gradloglikGaPExp(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt_r, cx, cy, Wxk] = loglikGaPreadparamHfix(x,varargin);

Hkt=exp(Hkt_r); %nonnegativity constrains
[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);

P=Wxkbg*Hktbg; %current approximation

%linear grasdient shifted by cx
xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx');
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk;
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% gfHkt=(alpha-1)*1./Hkt - 1/beta*Hkt + Wxk'*(Vxt./P)-1;
gfHkt= Hkt.*(Wxk'*(Vxt./P)-1); %without background (->not Wxkgb) and d(log(L)/dcx)
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');
% d(log(L))/dcy:
gfcy=diag(Wxtcy'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');

% gf = [reshape(gfHkt,1,peval.nt*peval.ncomp), gfcx', gfcy'];
% gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];
gf = [gfcx', gfcy'];
gf=-gf; %conjugate gradient is minimizing!
end

function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
        xxp=double(xx(sizevec, 'corner')); %linear function - pixels
    case 'yy'
        xxp=double(yy(sizevec, 'corner')); %linear function - pixels
    otherwise
        error('Wrong dir')
end
xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function  [plx] = get_all_from_plx(fileName)

plx = [];
if (nargin == 0 | isempty(fileName))
   [fname, pathname] = uigetfile('*.plx', 'Select a .plx file');
   if isequal(fname,0)
     error 'No file was selected'
     return
   end
   fileName = fullfile(pathname, fname);
end

[plx.tscounts, plx.wfcounts, plx.evcounts, plx.slowcounts] = plx_info(fileName,1);
[OpenedFilePath, plx.Version, plx.Freq, plx.Comment, plx.Trodalness, plx.NPW, plx.PreThresh, plx.SpikePeakV, plx.SpikeADResBits, plx.SlowPeakV, plx.SlowADResBits, plx.Duration, plx.DateTime] = plx_information(fileName);
[dirstr, name, ext] = fileparts(OpenedFilePath);
plx.FileName = [name ext];

% try valid and invalid channels and units
% max unit is 26
% max spike channel number is 128
for iunit = -1:30
    for ich = -1:130
        [plx.ts.n{iunit+2,ich+2}, plx.ts.ts{iunit+2,ich+2}] = plx_ts(fileName, ich , iunit );
        [plx.wf.n{iunit+2,ich+2}, plx.wf.npw{iunit+2,ich+2}, plx.wf.ts{iunit+2,ich+2}, plx.wf.wf{iunit+2,ich+2}] = plx_waves(fileName, ich , iunit );
        [plx.wf_v.n{iunit+2,ich+2}, plx.wf_v.npw{iunit+2,ich+2}, plx.wf_v.ts{iunit+2,ich+2}, plx.wf_v.wf{iunit+2,ich+2}] = plx_waves_v(fileName, ich , iunit );
     end
end

[nspk,plx.spk_filters] = plx_chan_filters(fileName);
[nspk,plx.spk_gains] = plx_chan_gains(fileName);
[nspk,plx.spk_threshs] = plx_chan_thresholds(fileName);
[nspk,plx.spk_names] = plx_chan_names(fileName);

% try valid and invalid channels
% ona 128-channel system, there could be 128*3 a/d channels
for ich = -1:400
	[plx.ad_gi.freq{ich+2}, plx.ad_gi.nad{ich+2}, plx.ad_gi.ts{ich+2}, plx.ad_gi.frag{ich+2}] = plx_ad_gap_info(fileName, ich);
	[plx.ad.freq{ich+2}, plx.ad.nad{ich+2}, plx.ad.ts{ich+2}, plx.ad.frag{ich+2}, plx.ad.val{ich+2}] = plx_ad(fileName, ich);
	[plx.ad_v.freq{ich+2}, plx.ad_v.nad{ich+2}, plx.ad_v.ts{ich+2}, plx.ad_v.frag{ich+2}, plx.ad_v.val{ich+2}] = plx_ad_v(fileName, ich);
	[plx.ad_span.freq{ich+2}, plx.ad_span.nad{ich+2}, plx.ad_span.val{ich+2}] = plx_ad_span(fileName, ich, 10,20);
	[plx.ad_span_v.freq{ich+2}, plx.ad_span_v.nad{ich+2} plx.ad_span_v.val{ich+2}] = plx_ad_span_v(fileName, ich, 10,20);
    % try invalid start/end combinations
	[plx.ad_span_inv1.freq{ich+2}, plx.ad_span_inv1.nad{ich+2} plx.ad_span_inv1.val{ich+2}] = plx_ad_span_v(fileName, ich, 10,-5);
    [plx.ad_span_inv2.freq{ich+2}, plx.ad_span_inv2.nad{ich+2} plx.ad_span_inv2.val{ich+2}] = plx_ad_span_v(fileName, ich, 10,9);
    [plx.ad_span_inv3.freq{ich+2}, plx.ad_span_inv3.nad{ich+2} plx.ad_span_inv3.val{ich+2}] = plx_ad_span_v(fileName, ich, -2,3);
end

[nad,plx.adfreqs] = plx_adchan_freqs(fileName);
[nad,plx.adgains] = plx_adchan_gains(fileName);
[nad,plx.adnames] = plx_adchan_names(fileName);

% the events
[u,plx.evnames] = plx_event_names(fileName);
[u,plx.evchans] = plx_event_chanmap(fileName);
for iev = -1:310
	[plx.nevs{iev+2}, plx.tsevs{iev+2}, plx.svdummy{iev+2}] = plx_event_ts(fileName, iev);
	if ( iev > 0 && iev <= length(plx.evchans) && plx.evcounts(iev) > 0 )
	    evch = plx.evchans(iev);
	    if ( evch == 257 )
			[plx.strobed.n, plx.strobed.ts, plx.strobed.v] = plx_event_ts(fileName, evch);
            [plx.nCoords, plx.nDim, plx.nVTMode, plx.c] = plx_vt_interpret(plx.strobed.ts, plx.strobed.v);
        end
    end
end

plx_close(fileName);

% try to get a/d channel after closing
for ich = -1:400
	[plx.ad2.freq{ich+2}, plx.ad2.nad{ich+2}, plx.ad2.ts{ich+2}, plx.ad2.frag{ich+2}, plx.ad2.val{ich+2}] = plx_ad(fileName, ich);
end
plx_close('');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function output = ...
  AnalyzeAperiodicity(x, fs, f0, frame_time, analysis_condition)
% F0 refinement and aperiodicity analysis
%
% output = ...
%   AnalyzeAperiodicity(x, fs, f0, frame_time, analysis_condition)
%
%  Input argument
%    x: speech signal, one dimensional column vector
%    fs: sampling frequency (Hz)
%    f0: initial F0 estimate (Hz)
%    frame_time: time location of the center of the analysis window
%    analysis_condition: structure variable with following fields
%      n_harmonics_for_H : number of harmonic components to use in
%        harmonics based F0 refinement
%      window_magnification_for_H : time window stretching facotor from
%        the permissible shortest length
%      enable_downsampling_for_H : 1: enables down sampling before analysis
%      conditions_for_T : structure variable array for defining condition
%        for time-warping based F0 refinement with the following fields
%        n_harmonics : number of harmonics to be uses
%        window_magnification : time window stretching facotor from
%          the permissible shortest length
%        enable_downsampling : 1: enables down sampling before analysis
%
%  Return value
%    output : structure variable with the following fields
%      aperiodicity_matrix : aperiodicity of each harmonic component
%      refined_f0 : refined F0 (Hz)
%      frame_time : analysis window center location at each frame (s)

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

narginchk(4, 5);
if nargin < 5, analysis_condition = GenerateAperiodicityAnalysisCondition; end;
start_tic = tic;
n_harmonics_for_H = analysis_condition.n_harmonics_for_H;
window_magnification_for_H = analysis_condition.window_magnification_for_H;
enable_downsampling_for_H = analysis_condition.enable_downsampling_for_H;
% Since at the very beginning F0 information is not very reliable,
% it is necessary to refine it without using F0 adaptive time warping
% Effects of time warping propagates to all the following process.
% The fundamental component is usually damaged by microphone low-cut for
% preventing proximity effect and background noise is stronger in lower
% frequency region. In addition, time warping procedure is fragile.
% It should be used carefully.
[refined_f0, ~] = RefineF0byHarmonics(x, fs, f0, frame_time, ...
  n_harmonics_for_H, window_magnification_for_H, enable_downsampling_for_H);
% Recursively apply time-warping based F0 refinement
conditions_for_T = analysis_condition.conditions_for_T;
n_steps = size(conditions_for_T, 2);
f0_last = refined_f0;
for ii = 1:n_steps
  [refined_f0, aperiodicity_matrix] = ...
    RefineF0byTimeWarping(x, fs, f0_last, frame_time, ...
                          conditions_for_T(ii).n_harmonics , ...
                          conditions_for_T(ii).window_magnification, ...
                          conditions_for_T(ii).enable_downsampling);
  f0_last = refined_f0;
end;
output = struct('aperiodicity_matrix', aperiodicity_matrix, ...
                'refined_f0', f0_last, 'sampling_frequency', fs, ...
                'frame_time', frame_time, 'elapsed_time', toc(start_tic));
end

author: google


text:
%
% subcortical_labeling_table.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.4 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

[Dvolume,Ivol]=subcortical_labeling_dir_afd('/space/neo/2/recon/buckner', 0.01)
% load('/space/okapi/3/data/laurence/ADF/subcortical_labeling/Percent_Vol0517.mat') %loads Dvolume
% load('/space/okapi/3/data/laurence/ADF/subcortical_labeling/Subj_nb_Vol0517.mat') %loads Ivol
files=dir('/space/neo/2/recon/buckner');
fp=fopen('~/subcortlab_table.txt', 'w');
sz=size(Dvolume);
fprintf(fp,'# FSAFD SubcorticalLabelingCheck 1\n# date 20050516\n# $Id: subcortical_labeling_table.m,v 1.4 2011/03/02 00:04:13 nicks Exp $\n# Buckner data set\n# ncols 20\n# nrows %d\n', sz(1));
ROI_Labels=[4 43 17 53 10 49 11 50 12 51 13 52 18 54 26 58 14 15 5 44];
labmapfile=('/space/lyon/1/fsdev/freesurfer_dev');
labmap=strcat(labmapfile,'/','tkmeditColorsCMA');
[label name val1 val2 val3 val4]=textread(labmap,'%d %s %d %d %d %d',89);
nn=1;
for u=1:length(ROI_Labels)
    fprintf(fp, '# label_col %d %s\n', u, char(name(ROI_Labels(u)+1)));
end
for i=1:sz(1)
     if(Dvolume(i,1)>0)
        nb=Ivol(i);
        fprintf(fp, '# label_row %d %s\n',nn, char(files(nb).name));
         nn=nn+1;
     end
end
for j=1:sz(1)
    for k=1:20
        fprintf(fp,'%f  ', Dvolume(j,k));
    end
    fprintf(fp, '\n');
end
fclose(fp)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef PTKContextSet
    % PTKContextSet. An enumeration used to specify the regions of interest
    %     over which a plugin can operate
    %
    %     A context set represents a set of contexts. A context is a region
    %     of interest for which the plugin can be executed. Therefore the
    %     context type repesents the set of contexts for which the plugin
    %     can be run.
    %
    %     For example, the "SingleLung" type represents two contexts: "LeftLung"
    %     and "RightLung". A plugin with the "SingleLung" type can only be run
    %     with one of these two input contexts. If a result is requested for
    %     another context, the Framework will try to convert the contexts
    %     appropriately (this happens in the MimContextHierarchy class).
    %
    %     For example, if a result is requested for the conext "LungROI" but the
    %     plugin being called has type "SingleLung", then the plugin will be
    %     called twice (once for "RightLung" and once for "LeftLung"). The
    %     final result will contain both results.
    %
    %     Plugins specify their PTKContextSet by setting their "Context"
    %     class property.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        OriginalImage,  % The full image
        LungROI,        % A region containing the lung and airways
        Lungs,          % Left and right lungs
        SingleLung,     % Left lung or right lung
        Lobe,
        Segment,

        Any             % Valid for any context
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%
for i = 1:p
    X_poly(:, i) = (X .^ i)';
end

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

end

author: ecmadao

repo: Coding-Guide


text:
% SIMPLE ALGORITHM TO FOLLOW A LINE IN SPACE. Error correction based on a P
% controller on the closest point to the line vector.
%
% Copyright (C) 2019, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function path_planning_line_simple
close all
% velocidad lineal entre puntos consecutivos
abs_linear_speed = 0.5; % (m/s)
delta_time = 0.05;
epsilon = 0.02;
kp = 10;

robot = load_robot('ABB', 'IRB52')
q = [0 0 0 0 0 0 0];

fprintf('\nPRESS ANY KEY TO CONTINUE...')
pause

%NOA matrix initial point
T1=[1 0 0 0.8;
    0 1 0 -0.5;
    0 0 1 0.3;
    0 0 0  1];
%NOA matrix end point
T2=[1 0 0 0.5;
    0 1 0 0.4;
    0 0 1 0.7;
    0 0 0  1];

start_point = T1(1:3,4);
end_point = T2(1:3,4);
% vector velocidad en la direcci�n de la trayectoria
v = (end_point-start_point);
v = abs_linear_speed*v/norm(v); %vector normalizado en la direcci�n de la recta
w = [0 0 0]';
xd = [v; w];

% Solve inverse kinematics at first position
qinv = inversekinematic(robot, T1, q);

%Select arbitrarily the first solution
q = qinv(:,1);
qs = [];
qds = [];

errors_line = [];
% using a standard inverse
while 1
   % Se calcula la Jacobiana del manipulador
   J = manipulator_jacobian(robot, q);
   qd1 = % hallar las velocidades articulares

   % nueva posición articular calculada
   q = q + % ...;
   % hallar la posición actual del extremo y guardarla en p
   T = %....;

   % Encuentra el error en el seguimiento de la recta
   % se le da esta parte realizada al alumno
   [delta_end, error_line, error_line_vector] = find_errors(start_point, end_point, p);
   % Esta parte se le da al alumno resuelta igualmente
   % se realiza una pequeña corrección para que el robot no se aleje del
   % seguimiento de la recta. Esta parte corrije el error inherente al
   % seguimiento de la trayectoria.
   % el alumno debe visualizar los errores
   % a) con las dos líneas siguientes comentadas
   % b) con las dos líneas siguientes sin comentar.
   %qd2 = inv(J)*[error_line_vector' 0 0 0]';
   %q = q + kp*delta_time*qd2;

   % El alumno debe comprobar si se ha llegado al final de la trayectoria y
   % salir del bucle

   % guardar datos para hacer plots posteriormente
   errors_line = [errors_line error_line];
   qs = [qs q];
   qds = [qds qd1];
   % pintar al robot
   drawrobot3d(robot, q)
   line([start_point(1) end_point(1)], [start_point(2) end_point(2)] , [start_point(3) end_point(3)] )
   % una pequeña pausa para permitir que el plot se actualice
   pause(0.1)
end

% ploteo de informacion
figure, plot(errors_line), title('Minimo error con la recta'), xlabel('Num. de movimiento'), ylabel('Error (m)')
figure, plot(qs'),  title('Coordenadas articulares'), xlabel('Movement number'), ylabel('Position (rad)')
figure, plot(qds'),  title('Velocidades articulares'), xlabel('Movement number'), ylabel('Speed (rad/s)')

% find:
% error_end: the error with respect to the end point
% error_line: the error of p to the line defined b point a and vector n.
% the line is defined by points a and b.
function [error_end, error_line, error_line_vector]=find_errors(a, b, p)
error_end = sqrt((p-b)'*(p-b));

% define the line as a, n
n = (b-a);
n = n/norm(n);

error_line = (a-p)-((a-p)'*n)*n;
error_line = norm(error_line);

% Now obtain current point minus initial point a
% project
w = p-a;
p_ = n*dot(w, n);
% add the origin since the point p_ is referred to the line
p_ = p_ + a;
% find a vector connecting the current point p and the point belonging to
% the line p_
error_line_vector = p_ - p;

repo: ARTE


text:
function f =  quantity(q, which)

%tstoolbox/@unit/quantity
%   returns quantity name of unit q. If argument which is omitted, the
%   english quantity name will be returned.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

if isa(q, 'unit')
	if nargin > 1
		if strcmp(which, 'ger')
			f = q.quantity.ger;
		else
			f = q.quantity.eng;
		end
	else
		f = q.quantity.eng;
	end
else
	f = '';
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% checkerBoard3D
% Below is a demonstration of the features of the |checkerBoard3D| function

%%
clear; close all; clc;

%% Syntax
% |M=checkerBoard3D(siz);|

%% Description
% This function creates a checkboard image of the size siz whereby elements
% are either black (0) or white (1). The first element is white.

%% Examples

%%
% Plot settings
fontSize=15;

%% Example creating a 2D checkerboard image

siz=[4 6]; %Image size
M=checkerBoard3D(siz); %Create checkerboard image

%%
% Plotting results

cFigure;
title('A 2D checkerboard pattern','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;
imagesc(M);

colormap gray;
axis equal; view(2); axis tight;
set(gca,'FontSize',fontSize);
drawnow;

%% Example creating a 3D checkerboard image

siz=[6 6 3]; %Image size
M=checkerBoard3D(siz); %Create checkerboard image

%%
% Plotting results

[Fv,Vv,Cv]=im2patch(M,1:numel(M),'vb'); %Create patch data for plotting

cFigure; hold on;
title('A 3D checkerboard pattern','FontSize',fontSize);
gpatch(Fv,Vv,Cv);

camlight('headlight');
colormap gray;
axisGeom;
set(gca,'FontSize',fontSize);
drawnow;

%% Example, changing checkerboard block size

siz=[50 50]; %Image size
blockSize=10; %Block size in pixel units
M=checkerBoard3D(siz,blockSize); %Create checkerboard image

%%
% Plotting results

[Fv,Vv,Cv]=im2patch(M,1:numel(M),'vb'); %Create patch data for plotting

cFigure; hold on;
title('A 2D checkerboard pattern','FontSize',fontSize);
gpatch(Fv,Vv,Cv,'r',1,1);

colormap gray;
axisGeom; view(2);
set(gca,'FontSize',fontSize);
drawnow;

%%

siz=[12 12 6]; %Image size
blockSize=3; %Block size in pixel units
M=checkerBoard3D(siz,blockSize); %Create checkerboard image

%%
% Plotting results

[Fv,Vv,Cv]=im2patch(M,1:numel(M),'vb'); %Create patch data for plotting

cFigure; hold on;
title('A 3D checkerboard pattern','FontSize',fontSize);
gpatch(Fv,Vv,Cv,'r',1,1);

camlight('headlight');
colormap gray;
axisGeom;
set(gca,'FontSize',fontSize);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program 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/>.


text:
% test the sparse updates of augmented Lagrange.
function [] = testMatrixMultiply2()
clc;

data  = [2 3  4  5  6 ];
Known = [1 11 13 20 44];

m = 5; n = 10; r = 3;

[Ik, Jk] = ind2sub([m n],Known);
S_data   = sparse(Ik, Jk, data, m, n);
rho = 0.5;

% initially S_comp = empty. X0 = U0*V0 + S_comp
S_comp =  sparse(Ik, Jk, eps * ones(size(data)), m, n);

rng(100);
U0 = rand(m, r); U0s = U0;
V0 = rand(r, n); V0s = V0;

check_func = @(NameStr, X, Xs) fprintf('[%s] Related Error %.4g\n', NameStr, norm(X-Xs, 'fro')/norm(X, 'fro'));
%check_func = @(NameStr, X, Xs) fprintf('[%s] Related Error %.4g\n', NameStr, norm(X-Xs, 'fro'));

X0 = (U0 * V0 + S_comp);

disp('init lambda');
[X1  U1  V1]         = update_initLambda (X0, U0, V0, Known, data, rho);
[X1s U1s V1s S_comp] = update_initLambdaSparse (S_comp, U0s, V0s, Ik, Jk, data, rho, true);

check_func('U+', U1, U1s);
check_func('V+', V1, V1s);
check_func('X+', X1, X1s);
for i = 1: 2
    [X1  U1  V1]         = update_initLambda (X1, U1, V1, Known, data, rho);
    [X1s U1s V1s S_comp] = update_initLambdaSparse (S_comp, U1s, V1s, Ik, Jk, data, rho, false);
    check_func('U+', U1, U1s);
    check_func('V+', V1, V1s);
    check_func('X+', X1, X1s);
end

disp('given sparse lambda');
U0 = rand(m, r); U0s = U0;
V0 = rand(r, n); V0s = V0;
X0 = (U0 * V0 + S_comp);

Lambda = sparse(Ik, Jk, ones(size(data)), m, n);
[X1 U1 V1] = update_givenLambda (X0, U0, V0, Known, data, rho, Lambda);
[X1s U1s V1s S_comp] = update_givenLambdaSparse (S_comp, U0s, V0s, Ik, Jk, data, rho, Lambda, true);
check_func('U+', U1, U1s);
check_func('V+', V1, V1s);
check_func('X+', X1, X1s);
for i = 1: 2
    [X1  U1  V1]         = update_givenLambda (X1, U1, V1, Known, data, rho, Lambda);
    [X1s U1s V1s S_comp] = update_givenLambdaSparse (S_comp, U1s, V1s, Ik, Jk, data, rho, Lambda, false);
    check_func('U+', U1, U1s);
    check_func('V+', V1, V1s);
    check_func('X+', X1, X1s);
end

    function [X1 U1 V1] = update_initLambda (X0, U0, V0, Known, data, rho)
        Lambda0 = ones(m, n);
        % X, V => U
        T = X0 - Lambda0/rho;
        U1  = T * V0';
        % X, U => V
        V1  = U1' * T;
        % U, V => X
        T2 = Lambda0/rho + U1 * V1;
        T3 = (2 * data + rho * T2(Known))/(2+ rho);
        X1 = T2; X1(Known) = T3;
    end

    function [X1s U1s V1s S_comp] = update_initLambdaSparse (S_comp, U0, V0, Ik, Jk, data, rho, first_iter)
        % Lambda is given by Lambda = e * et
        e  = ones(m, 1);
        et = ones(1, n);
        if first_iter
            % X, V => U
            U1s =  U0 * (V0 * V0') + S_comp  * V0' -  e * (et * V0') /rho;
            % X, U => V
            V1s = (U1s' * U0) * V0 + U1s' * S_comp  - (U1s' * e) * et / rho;
            % U, V => X
            s_comp_val = 2/(2+rho) * ( data - ones(size(data))/ rho - sparse_inp(U1s', V1s, Ik, Jk));
        else
            % X, V => U
            U1s =  U0 * (V0 * V0') + S_comp  * V0';
            % X, U => V
            V1s = (U1s' * U0) * V0 + U1s' * S_comp;
            % U, V => X
            s_comp_val = 2/(2+rho) * ( data  - sparse_inp(U1s', V1s, Ik, Jk));
        end
        sparse_update(S_comp, s_comp_val);
        X1s = e * et/rho + U1s * V1s + S_comp;
    end

    function [X1 U1 V1] = update_givenLambda (X0, U0, V0, Known, data, rho, Lambda)
        % X, V => U
        T = X0 - Lambda/rho;
        U1  = T * V0';
        % X, U => V
        V1  = U1' * T;
        % U, V => X
        T2 = Lambda/rho + U1 * V1;
        T3 = (2 * data + rho * T2(Known))/(2+ rho);
        X1 = T2; X1(Known) = T3;
    end

    function [X1s U1s V1s S_comp] = update_givenLambdaSparse (S_comp, U0, V0, Ik, Jk, data, rho, LamSMat, first_iter)
        % Lambda is a sparse matrix here.
        if first_iter
            % X, V => U
            U1s =  U0 * (V0 * V0') + S_comp  * V0' -  (LamSMat * V0') /rho;
            % X, U => V
            V1s = (U1s' * U0) * V0 + U1s' * S_comp  - (U1s' * LamSMat) / rho;
            % U, V => X
            [~, ~, ll] = find(LamSMat);
            s_comp_val = 2/(2+rho) * ( data - ll'/ rho - sparse_inp(U1s', V1s, Ik, Jk));
        else
            % X, V => U
            U1s =  U0 * (V0 * V0') + S_comp  * V0';
            % X, U => V
            V1s = (U1s' * U0) * V0 + U1s' * S_comp;
            % U, V => X
            s_comp_val = 2/(2+rho) * ( data  - sparse_inp(U1s', V1s, Ik, Jk));
        end
        sparse_update(S_comp, s_comp_val);
        X1s = LamSMat/rho + U1s * V1s + S_comp;
    end
end

author: andrewssobral

repo: lrslibrary


text:
%kMarker 'Marker Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Marker.pane file
%
% Parameters:
% Integer: graphicsMarkerScale 'Marker Scale', default: 1: 'marker scale'
%
% Example: kMarker( {'graphicsMarkerScale',1})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Marker Visual Object"
% A marker visual object supports the display of a marker annotation.
% The (x, y) location of the marker may be specified in device coordinates
% or world coordinates, as determined by the application.  The menuform,
% however, displays the world coordinate location of the marker,
% regardless of whether the application used world or device coordinates to
% specify that location.  Attributes of the marker object that can be
% interactively specified include its type, scale, and color, as well as its
% world coordinate location.
%
% .item "Marker Type"
% The indicator may use a variety of marker types, including Arc, Bow Tie, Box,
% Caret, Circle, Cross, Dagger, Diamond, Dot, Hexagon, Point, Square, Triangle,
% V, and X.
%
% .item "Marker Scale"
% Markers can be scaled so that they are displayed in a larger size.  The scale
% is an integer value;  for a scale of 2, the marker will appear twice its
% normal size, for a scale of 3, the marker will appear three times its normal
% size, and so on.
%
% .item "Marker Color"
% This stringlist selection lets you set the foreground color
% of the marker object, using the color name.
%
% .item "X Placement, Y Placement"
% This pair of double values specifies the (x, y) location of the marker object
% on its parent in world coordinates.  The world coordinate range of the marker
% is specified by the application, and thus varies depending on context.
%
% .end tagged

function varargout = kMarker(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kMarker(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'graphicsMarkerScale', 1};
maxval={100};
minval={1};
istoggle=[0];
was_set=istoggle * 0;
paramtype={'Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'Marker" '],Inputs,narglist);


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

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

% Last Modified by GUIDE v2.5 15-Jul-2020 14:54:25

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

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

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

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

handles.fields = varargin{1};
handles.FieldFlags = varargin{2};

%populate available and selected fields
Count1 = 0;
Count2 = 0;
AvailFields = {};
SelFields = {};
for i=1:length(handles.fields)
    if handles.FieldFlags(i)
        Count2 = Count2 + 1;
        SelFields{Count2} = handles.fields{i};
    else
        Count1 = Count1 + 1;
        AvailFields{Count1} = handles.fields{i};
    end
end

set(handles.AvailableFields,'String',AvailFields);
set(handles.SelectedFields,'String',SelFields);

if isempty(AvailFields)
    set(handles.AddField,'Enable','off');
    set(handles.AddAll,'Enable','off');
end

if isempty(SelFields)
    set(handles.RemoveField,'Enable','off');
    set(handles.RemoveAll,'Enable','off');
end

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.fields;
varargout{2} = handles.FieldFlags;
delete(handles.figure1);

% --- Executes on selection change in AvailableFields.
function AvailableFields_Callback(hObject, eventdata, handles)
% hObject    handle to AvailableFields (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns AvailableFields contents as cell array
%        contents{get(hObject,'Value')} returns selected item from AvailableFields

% --- Executes during object creation, after setting all properties.
function AvailableFields_CreateFcn(hObject, eventdata, handles)
% hObject    handle to AvailableFields (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in SelectedFields.
function SelectedFields_Callback(hObject, eventdata, handles)
% hObject    handle to SelectedFields (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns SelectedFields contents as cell array
%        contents{get(hObject,'Value')} returns selected item from SelectedFields

% --- Executes during object creation, after setting all properties.
function SelectedFields_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SelectedFields (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in AddField.
function AddField_Callback(hObject, eventdata, handles)
% hObject    handle to AddField (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%get selection, remove from available fields and add to selected fields
Strings = get(handles.AvailableFields,'String');
Val = get(handles.AvailableFields,'Value');
StringsNew = {};
count = 0;
for i=1:length(Strings)
    if i~=Val
        count = count + 1;
        StringsNew{count} = Strings{i};
    end
end
set(handles.AvailableFields,'Value',1);
set(handles.AvailableFields,'String',StringsNew);
if isempty(StringsNew)
    set(handles.AddField,'Enable','off');
    set(handles.AddAll,'Enable','off');
else
    set(handles.AddField,'Enable','on');
    set(handles.AddAll,'Enable','on');
end

set(handles.RemoveField,'Enable','on');
set(handles.RemoveAll,'Enable','on');

Strings2 = get(handles.SelectedFields,'String');
Strings2{end+1} = Strings{Val};
set(handles.SelectedFields,'String',Strings2);

% --- Executes on button press in RemoveField.
function RemoveField_Callback(hObject, eventdata, handles)
% hObject    handle to RemoveField (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%get selection, remove from available fields and add to selected fields
Strings = get(handles.SelectedFields,'String');
Val = get(handles.SelectedFields,'Value');
StringsNew = {};
count = 0;
for i=1:length(Strings)
    if i~=Val
        count = count + 1;
        StringsNew{count} = Strings{i};
    end
end
set(handles.SelectedFields,'Value',1);
set(handles.SelectedFields,'String',StringsNew);
if isempty(StringsNew)
    set(handles.RemoveField,'Enable','off');
    set(handles.RemoveAll,'Enable','off');
else
    set(handles.RemoveField,'Enable','on');
    set(handles.RemoveAll,'Enable','on');
end
set(handles.AddField,'Enable','on');
set(handles.AddAll,'Enable','on');

Strings2 = get(handles.AvailableFields,'String');
Strings2{end+1} = Strings{Val};
set(handles.AvailableFields,'String',Strings2);

% --- Executes on button press in AddAll.
function AddAll_Callback(hObject, eventdata, handles)
% hObject    handle to AddAll (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Strings = get(handles.AvailableFields,'String');
Strings2 = get(handles.SelectedFields,'String');

set(handles.AvailableFields,'Value',1);
set(handles.AvailableFields,'String',[]);
set(handles.SelectedFields,'String',vertcat(Strings2,Strings));

set(handles.AddField,'Enable','off');
set(handles.AddAll,'Enable','off');

set(handles.RemoveField,'Enable','on');
set(handles.RemoveAll,'Enable','on');

% --- Executes on button press in RemoveAll.
function RemoveAll_Callback(hObject, eventdata, handles)
% hObject    handle to RemoveAll (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Strings = get(handles.AvailableFields,'String');
Strings2 = get(handles.SelectedFields,'String');

set(handles.SelectedFields,'Value',1);
set(handles.SelectedFields,'String',[]);
set(handles.AvailableFields,'String',vertcat(Strings,Strings2));

set(handles.RemoveField,'Enable','off');
set(handles.RemoveAll,'Enable','off');

set(handles.AddField,'Enable','on');
set(handles.AddAll,'Enable','on');

% --- Executes on button press in OK.
function OK_Callback(hObject, eventdata, handles)
% hObject    handle to OK (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%update fields and FieldFlags based on the list box contents
SelFields = get(handles.SelectedFields,'String');
AvailFields = get(handles.AvailableFields,'String');

handles.fields = vertcat(SelFields,AvailFields);
handles.FieldFlags = [ones(length(SelFields),1); zeros(length(AvailFields),1)];

guidata(hObject, handles);

uiresume(handles.figure1);

% --- Executes on button press in Cancel.
function Cancel_Callback(hObject, eventdata, handles)
% hObject    handle to Cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

uiresume(handles.figure1);

author: ngageoint


text:
function D = agreement_weighted(CI,Wts)
%WEIGHTED_AGREEMENT     weights agreement matrix
%
%   D = AGREEMENT_WEIGHTED(CI,WTS) is identical to AGREEMENT, with the
%   exception that each partitions contribution is weighted according to
%   the corresponding scalar value stored in the vector WTS. As an example,
%   suppose CI contained partitions obtained using some heuristic for
%   maximizing modularity. A possible choice for WTS might be the Q metric
%   (Newman's modularity score). Such a choice would add more weight to
%   higher modularity partitions.
%
%   NOTE: Unlike AGREEMENT, this script does not have the input argument
%   BUFFSZ.
%
%   Inputs:     CI,     set of partitions
%               WTS,    relative weight of importance of each paritition
%
%   Outputs:    D,      weighted agreement matrix
%
%   Richard Betzel, Indiana University, 2013

Wts = Wts./sum(Wts);
[N,M] = size(CI);
D = zeros(N);
for i = 1:M
    d = dummyvar(CI(:,i));
    D = D + (d*d')*Wts(i);
end


text:
function pass = test_abs( )
% Test abs in SPHEREFUN

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

f = spherefun(@(x,y,z) -(x.^2 + y.^2 + z.^2));
pass(1) = norm(abs(f) + f, inf) < tol;

end


text:
% Fit DCM-phase model to MEG data
%
% With the DCM-GUI one is restricted to models
% with phase interactions that are sin functions.
% With scripts, this can be relaxed to arbitrary order
% Fourier series
%
% Additionally with the DCM-GUI one is
% restricted to only one or two experimental conditions
% This is relaxes in a script

% MEG data file
DCM.xY.Dfile='DIRECTORY\fdas_delay_run1.mat';

NFourier=1;  % Order of Fourier expansion for PIFs

% Conditions 1 and 3; 'control' and 'hard'
DCM.options.trials=[1,3];
DCM.xU.X=[0 1]';

% Just R-MTL, R-VIS, R-IFG
source_pos(1,:)=[27,-18,-27];
source_pos(2,:)=[10,-100,0];
source_pos(3,:)=[39,28,-12];
DCM.Sname{1}='Right-MTL';
DCM.Sname{2}='Right-VIS';
DCM.Sname{3}='Right-IFG';
Nr=length(DCM.Sname);

DCM.Lpos=source_pos';
DCM.options.spatial='ECD';
DCM.xY.modality='MEG';

D=spm_eeg_load(DCM.xY.Dfile);
Ic = D.indchantype(DCM.xY.modality, 'GOOD');
DCM.xY.Ic       = Ic;

% forward model (spatial)
%------------------------------------
DCM = spm_dcm_erp_dipfit(DCM,1);
DCM.M.dipfit.type='ECD';

% Define time points to model
DCM.options.Tdcm=[5000 6000];

% Project MEG data onto source locations and average
DCM = spm_dcm_phase_data(DCM);

% Define DCM structure
A=[0 1 0; 0 0 0; 0 1 0];
for n=1:NFourier,
    DCM.As(:,:,n)=A;
    DCM.Ac(:,:,n)=A;
    DCM.Bs{1}(:,:,n)=A;
    DCM.Bc{1}(:,:,n)=A;
end

Ns=size(DCM.xY.y{1},1);

% Define Frequency Range
DCM.options.Fdcm=[4 8];

DCM.xU.u=zeros(Ns,1);

% Fit model
tic;
DCM=spm_dcm_phase(DCM);
toc

% Show results
disp('');
disp('Freqs:');
DCM.Ep.df

disp('Control parameters');
disp('Fitted');
DCM.Ep.As
DCM.Ep.Ac

disp('Memory parameters');
DCM.Ep.As+DCM.Ep.Bs{1}
DCM.Ep.Ac+DCM.Ep.Bc{1}

author: spm


text:
function [irf_record,D_record,gamma_record,struct_irf_record,irf_estimates,D_estimates,gamma_estimates,strshocks_record,strshocks_estimates]=panel4irf_mean(Ymat,Xmat,beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,IRFband,N,n,m,p,k,T,IRFt,signrestable,signresperiods)

% because there is only one model, there is no need to loop over units
% hence, estimate the IRFs for the model, as with a standard normal-Wishart
% first run the Gibbs sampler to obtain posterior draws
[irf_record]=bear.irf(beta_gibbs,It,Bu,IRFperiods,n,m,p,k);

% If IRFs have been set to an unrestricted VAR (IRFt=1):
if IRFt==1
% run a pseudo Gibbs sampler to obtain records for D and gamma (for the trivial SVAR)
[D_record gamma_record]=bear.irfunres(n,It,Bu,sigma_gibbs);
struct_irf_record=[];
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(irf_record,n,IRFperiods,IRFband,IRFt,[],[]);

% If IRFs have been set to an SVAR with Choleski identification (IRFt=2):
elseif IRFt==2
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record D_record gamma_record]=bear.irfchol(sigma_gibbs,irf_record,It,Bu,IRFperiods,n);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record);

% If IRFs have been set to an SVAR with triangular factorisation (IRFt=3):
elseif IRFt==3
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record D_record gamma_record]=bear.irftrig(sigma_gibbs,irf_record,It,Bu,IRFperiods,n);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record);

% if IRFs have been set to an SVAR with sign restrictions
elseif IRFt==4
%if Magres==1
%[struct_irf_record D_record gamma_record]=bear.irfres_relmagnitude_panel(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,m,p,k,signrestable,signresperiods,relmagrestable, relmagresperiods);
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
%end
[struct_irf_record,D_record,gamma_record]=bear.irfsignrespanel(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,p,m,k,signrestable,signresperiods);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record);
end

% also, if a s structural identification was implemented, compute structural shocks
strshocks_record={};
strshocks_estimates={};
if IRFt~=1
   % because shocks have to be computed for each unit, loop over units
   for ii=1:N
   % run the Gibbs sampler
   strshocks_record(:,:,ii)=bear.strshocks(beta_gibbs,D_record,Ymat(:,:,ii),Xmat(:,:,ii),n,k,It,Bu);
   % obtain point estimates and credibility intervals
   strshocks_estimates(:,:,ii)=bear.strsestimates(strshocks_record(:,:,ii),n,T,IRFband);
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function C = horzcat(varargin)
%HORZCAT   Horizontally concatenate chebmatrices.
%   Z = [A,B,C,...] horizontally combines the chebmatrices, operator
%   blocks, chebfuns, and scalars given in the call, if their row sizes
%   are compatible.
%
% See also CHEBMATRIX.CAT, CHEBMATRIX.VERTCAT.

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

% Remove empty arguments.
isemp = cellfun(@isempty,varargin);
varargin(isemp) = [];
nargin = length(varargin);

% First, process the blocks of the inputs, disregarding sizes. This
% simplifies the logic. The entries of the results are cells holding the
% relevant blocks to be concatenated.
blocks = cell(1, nargin);
for n = 1:nargin
    item = varargin{n};
    if ( isa(item, 'chebmatrix') )   % already a cell
        blocks{n} = item.blocks;
    elseif ( isa(item, 'linBlock') || isa(item,'chebfun') )
        blocks{n} = { item };        % encase in a cell
    elseif ( ~isempty(item) )
        blocks{n} = num2cell(item);  % encase in cell(s)
    end
end

% Now check the row sizes. Create an empty vessel if OK.
[m,n] = cellfun(@size,blocks);
if ( all(m==m(1)) )
    B = cell( m(1), sum(n) );
else
    error('CHEBFUN:CHEBMATRIX:horzcat:sizeMismatch', ...
        'Incompatible row sizes.')
end

% Now we need to flatten out all the inner nested cell divisions, leaving
% just a cell of blocks.
cs = [0 cumsum(n)];
for j = 1:nargin
    B(:,cs(j)+(1:n(j))) = blocks{j};
end

% This step will perform domain compatibility checking.
C = chebmatrix( B );

% Finally, check block size compatibility.
[row, col] = blockSizes(C);
if ( any( any( bsxfun(@ne, row(:,1), row) ) ) )
    error('CHEBFUN:CHEBMATRIX:horzcat:blockSizeMismatch', ...
        'Block row sizes must be the same.')
end

end

author: chebfun


text:
%% Band-Pass filter visualization

[b, a] = butter(4, [0.5 50] / 100, 'bandpass');
y = filter(b, a, cnt.x);
cnt.x=cnt.x(:,1:34);
cnt.clab=cnt.clab(:,1:34);

f_sz=ceil(length(cnt.x)/2);
f=100*linspace(0,1,f_sz);
f_X=fft(cnt.x);
f_y=fft(y);
subplot(2,1,1)
stem(f,abs(f_X(1:f_sz)));
title('Original signal');
xlabel('frequency');
xlim([0 50]);
ylabel('power');
subplot(2,1,2)
stem(f,abs(f_y(1:f_sz)));
xlim([0 50]);

title('Application of the beta (13-30Hz) bandpass filter')
xlabel('frequency');
ylabel('power');

author: PatternRecognition

repo: OpenBMI


text:
function asa226_test ( )

%*****************************************************************************80
%
%% ASA226_TEST tests ASA226.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA226_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA226 library.\n' );

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function parameters = cohenimcra_estimation(ns_ps,parameters);
%         parameters = struct('n',2,'len',len_val,'noise_cap',ns_ps,'noise_tild',ns_ps,'gamma',ones(len_val,1),'Sf',Sf_val,...
%             'Smin',Sf_val,'S',Sf_val,'S_tild',Sf_val,'GH1',ones(len_val,1),'Smin_tild',Sf_val,'Smin_sw',Sf_val,'Smin_sw_tild',Sf_val,...
%             'stored_min',max(ns_ps)*ones(len_val,U_val),'stored_min_tild',max(ns_ps)*ones(len_val,U_val)','u1',1,'u2',1,'j',2,...
%             'alpha_d',0.85,'alpha_s',0.9,'U',8,'V',15,'Bmin',1.66,'gamma0',4.6,'gamma1',3,'psi0',1.67,'alpha',0.92,'beta',1.47,...
%             'b',b_val,'Sf_tild',Sf_tild_val);

n = parameters.n;
len = parameters.len;
gamma0 = parameters.gamma0;
psi0 = parameters.psi0;
alpha = parameters.alpha;
beta = parameters.beta;
b = parameters.b;
U = parameters.U;
V = parameters.V;
Bmin = parameters.Bmin;
gamma = parameters.gamma;
gamma1 = parameters.gamma1;
GH1 = parameters.GH1;
Sf = parameters.Sf;
Smin = parameters.Smin;
S = parameters.S;
S_tild = parameters.S_tild;
Sf_tild = parameters.Sf_tild;
Smin_tild = parameters.Smin_tild;
Smin_sw = parameters.Smin_sw;
Smin_sw_tild = parameters.Smin_sw_tild;
stored_min = parameters.stored_min;
stored_min_tild = parameters.stored_min_tild;
u1 = parameters.u1;
u2 = parameters.u2;
j = parameters.j;
alpha_d = parameters.alpha_d;
alpha_s = parameters.alpha_s;

noise_cap = parameters.noise_ps;
noise_tild = parameters.noise_tild;

gamma_old=gamma;
gamma=ns_ps./noise_cap;                                     %Eq 3 to compute a posteriori SNR
eps_cap=alpha*(GH1.^2).*gamma_old+(1-alpha)*max(gamma-1,0); %Eq 32 to compute a priori SNR
v=gamma.*eps_cap./(1+eps_cap);
exp_int=expint(v); %ei(v);
GH1=eps_cap.*exp(.5*exp_int)./(1+eps_cap);                    %Eq 33 to compute the value of GH1
Sf(1)=ns_ps(1);Sf(end)=ns_ps(end);
for f=2:len-1
    Sf(f)=sum(b.*[ns_ps(f-1); ns_ps(f); ns_ps(f+1)]);
end

S=alpha_s*S+(1-alpha_s)*Sf;                                 %Eq 14 and 15 for computing S(k,l)
Smin=min(Smin,S);Smin_sw=min(Smin_sw,S);
gamma_min=ns_ps./(Bmin*Smin);
psi=S./(Bmin*Smin);                                         %Eq 18 to compute gamma_min  and Psi
I=zeros(len,1);
index=find(gamma_min<gamma0 & psi<psi0);
I(index)=1;                                                 %Eq 21 to compute I(k,l)
for f=2:len-1
    if (I(f-1)+I(f)+I(f+1))==0
        Sf_tild(f)=S_tild(f);
    else
        Sf_tild(f)=sum(b.*[I(f-1); I(f); I(f+1)].*[ns_ps(f-1); ns_ps(f); ns_ps(f+1)])/sum(b.*[I(f-1); I(f); I(f+1)]);
    end                                                                   %Eq 26 for updating Sf_tild
end
if I(1)==0
    Sf_tild(1)=S_tild(1);
    Sf_tild(end)=S_tild(end);
else
    Sf_tild(1)=ns_ps(1);
    Sf_tild(end)=ns_ps(end);
end

S_tild=alpha_s*S_tild+(1-alpha_s)*Sf_tild;                              %Eq 27 for updating S_tild
Smin_tild=min(Smin_tild,S_tild);Smin_sw_tild=min(Smin_sw_tild,S_tild);
gamma_min_tild=ns_ps./(Bmin*Smin_tild);
psi_tild=S./(Bmin*Smin_tild);
q=zeros(len,1);                                                  %Eq 29 to find a priori speech absence probability
index=find(gamma_min_tild<=1 & psi_tild<psi0);
index1=setdiff([1:len],index);
if (~isempty(index))
    q(index)=1;
end
index=find(gamma_min_tild>1 & gamma_min_tild<gamma1 & psi_tild<psi0);
if (~isempty(index))
    q(index)=(gamma1-gamma_min_tild(index))/(gamma1-1);
end
%        p=1./(1+((q./(1-q)).*(1+eps_cap).*exp(-v)));
p=zeros(len,1);
if (~isempty(index1))
    temp1 = q(index1)./(1-q(index1));
    temp2 = 1 + eps_cap(index1);
    temp3 = exp(-v(index1));
    p(index1) = (1 + temp1.*temp2.*temp3).^-1;
end
%             p(index1)=1./(1+((q(index1)./(1-q(index1))).*(1+eps_cap(index1)).*exp(-v(index1))));    %Eq 7 to find conditional speech presence probability
alpha_d_tild=alpha_d+(1-alpha_d)*p;                                 %Eq 11 to update the time and frequency dependent smoothing factor
noise_tild=alpha_d_tild.*noise_tild+(1-alpha_d_tild).*ns_ps;        %Eq 10 to update noise estimate
noise_cap=beta*noise_tild;                                          %Eq 12 to update noise estimate for correction factor
j=j+1;
if j==V

    stored_min(:,u1)=Smin_sw;
    u1=u1+1;if u1==U+1; u1=1;end
    Smin=min(stored_min,[],2);
    Smin_sw=S;
    stored_min_tild(:,u2)=Smin_sw_tild;
    u2=u2+1;if u2==U+1; u2=1;end
    Smin_tild=min(stored_min_tild,[],2);
    Smin_sw_tild=S_tild;
    j=0;
end
ix=find(p>1 | q>1 | p<0 | q<0);
if (~isempty(ix))
    keyboard;
end
noise_ps = noise_cap;

parameters.n = n+1;
parameters.gamma = gamma;
parameters.GH1 = GH1;
parameters.Sf = Sf;
parameters.Smin = Smin;
parameters.S = S;
parameters.S_tild = S_tild;
parameters.Smin_tild = Smin_tild;
parameters.Smin_sw = Smin_sw;
parameters.Smin_sw_tild = Smin_sw_tild;
parameters.stored_min = stored_min;
parameters.stored_min_tild = stored_min_tild;
parameters.u1 = u1;
parameters.u2 = u2;
parameters.j = j;
parameters.noise_tild=noise_tild;
parameters.noise_ps = noise_ps;

author: jtkim-kaist

repo: Speech-enhancement


text:
function prob_test005 ( )

%*****************************************************************************80
%
%% TEST005 tests ANGLIT_MEAN, ANGLIT_SAMPLE, ANGLIT_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST005\n' );
  fprintf ( 1, '  For the Anglit PDF:\n' );
  fprintf ( 1, '  ANGLIT_MEAN computes the mean;\n' );
  fprintf ( 1, '  ANGLIT_SAMPLE samples;\n' );
  fprintf ( 1, '  ANGLIT_VARIANCE computes the variance.\n' );

  mean = anglit_mean ( );
  variance = anglit_variance ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF mean =     %14f\n', mean );
  fprintf ( 1, '  PDF variance = %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = anglit_sample ( seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
%% patch2obj
% Below is a demonstration of the features of the |patch2obj| function

%%
clear; close all; clc;

%% Syntax
% |patch2obj(objFileName,F,V);|
% |patch2obj(objFileName,F,V,C);|
% |patch2obj(objFileName,F,V,C,cMap);|
% |patch2obj(objFileName,F,V,C,cMap,cLim);|
% |patch2obj(objFileName,F,V,C,cMap,cLim,mtlStruct);|

%% Description
% This function exports the patch data defined by the faces (F), vertices
% (V) and the color data (C) to the OBJ (Wavefront .obj) format. The
% function generates a .obj file, a .mtl file, and a .jpg image file. The
% .obj file contains the geometry information and texture/color coordinates
% to use. The .mtl file contains the material information and refers to the
% image to use to look up the colors based on the texture coordinates in
% the .obj file.
% The color data C should ideally define either the vertex or face colors
% in the form of an nx1 array. If face colors are provided these are
% re-sampled (averaged) to vertex colors which is the required format for
% OBJ files. Colors are obtained from the input color map as well as the
% color limits. The input structure mtlStruct defines the MTL file
% components. With the default entries:
%
% mtlStruct.Ka=[1 1 1]; %Ambient color
% mtlStruct.Kd=[1 1 1]; %Diffuse color
% mtlStruct.Ks=[0 0 0]; %Specular color, black=off
% mtlStruct.Ns=0; %Specular component [0-1000]
% mtlStruct.Ni=1.45; %Optical density/index of refraction
% mtlStruct.d=1; %"dissolved"/transparancy [0-1]
% mtlStruct.Tr=0; %1 - d, used instead of d by some software
% mtlStruct.illum=1; %Illumination model
%
% Illumination models:
% 0. Color on and Ambient off
% 1. Color on and Ambient on
% 2. Highlight on
% 3. Reflection on and Ray trace on
% 4. Transparency: Glass on, Reflection: Ray trace on
% 5. Reflection: Fresnel on and Ray trace on
% 6. Transparency: Refraction on, Reflection: Fresnel off and Ray trace on
% 7. Transparency: Refraction on, Reflection: Fresnel on and Ray trace on
% 8. Reflection on and Ray trace off
% 9. Transparency: Glass on, Reflection: Ray trace off
% 10. Casts shadows onto invisible surfaces
%
%
% For more information on the OBJ file format see:
% https://en.wikipedia.org/wiki/Wavefront_.obj_file
% http://paulbourke.net/dataformats/obj/minobj.html

%% EXAMPLES

clear; close all; clc;

%% Example 1: Export colored patch data to the OBJ format

%Define patch data
testCase=1;
switch testCase
    case 1 %David
        [F,V]=graphicsModels(9);
        t=V(:,1)-min(V(:,1));
        t=t./max(t(:));
        C=sin(2*t*2*pi);
        C=abs(C);
        cMap=gjet(250); %Define colormap
    case 2
        [X,Y,Z]=peaks(25);
        [F,V,~]=grid2patch(X,Y,Z,Z);
        C=V(:,3);
        cMap=turbo(250); %Define colormap
    case 3 %Femur
        [F,V]=graphicsModels(5);
        C=V(:,1);
        cMap=turbo(250); %Define colormap
    case 4
        [F,V]=stanford_bunny;
        C=V(:,1);
        cMap=viridis(250); %Define colormap
end

%Define file name
gibbonFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(gibbonFolder,'data','OBJ');

fileName=fullfile(savePath,'test.obj');

%%
% Visualiza patch data

cFigure; hold on;
title('MATLAB patch','FontSize',25);
hp=gpatch(F,V,C,'none'); hp.FaceColor='interp';
axisGeom;
colormap(cMap);
camlight headlight;
gdrawnow;

%%
% Export to obj

patch2obj(fileName,F,V,C,cMap);

%%
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program 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/>.

author: gibbonCode

repo: GIBBON


text:
%FEATSELF Trainable mapping for forward feature selection
%
%   [W,R] = FEATSELF(A,CRIT,K,T)
%   [W,R] = A*FEATSELF([],CRIT,K,T)
%   [W,R] = A*FEATSELF(CRIT,K,T)
%   [W,R] = FEATSELF(A,CRIT,K,N)
%   [W,R] = A*FEATSELF([],CRIT,K,N)
%   [W,R] = A*FEATSELF(CRIT,K,N)
%
% INPUT
%   A    Training dataset
%   CRIT Name of the criterion or untrained mapping
%        (default: 'NN', i.e. the LOO 1-Nearest Neighbor error)
%   K    Number of features to select (default: K = 0, return optimal set)
%   T    Tuning dataset (optional)
%   N    Number of cross-validations (optional)
%
% OUTPUT
%   W    Output feature selection mapping
%   R    Matrix with step-by-step results
%
% DESCRIPTION
% Forward selection of K features using the dataset A. CRIT sets the
% criterion used by the feature evaluation routine FEATEVAL. If the
% dataset T is given, it is used as test set for FEATEVAL. Alternatvely a
% a number of cross-validation N may be supplied. For K = 0, the optimal
% feature set (corresponding to the maximum value of FEATEVAL) is returned.
% The result W can be used for selecting features using B*W.
% The selected features are stored in W.DATA and can be found by +W.
% In R, the search is reported step by step as:
%
% 	R(:,1) : number of features
% 	R(:,2) : criterion value
% 	R(:,3) : added / deleted feature
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELLR, FEATSEL,
% FEATSELO, FEATSELB, FEATSELI, FEATSELP, FEATSELM

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: featself.m,v 1.3 2007/04/21 23:06:46 duin Exp $

function [w,r] = featself(varargin)

  varargin = shiftargin(varargin,{'char','prmapping'});
  argin = setdefaults(varargin,[],'NN',0,[],[]);
  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained','Forward FeatSel');
    return
  end

  [a,crit,ksel,t,fid] = deal(argin{:});
  [w,r] = featsellr(a,crit,ksel,1,0,t);

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit