text:
function status=mrInitRetFromFSL_NIC(procScans,SessionName,subjectName)
%   status=mrInitRetFromFSL(foundScans,SessionName,subjectName)
%   creates a MLR session out of a mixture of Dicom and analyze files
%   this mixture is the result of the prepareUCSF script that preprocesses
%   the functional data from the UCSF scanner, but leaves the inplanes
%   as Dicom.
%   2006.04.20 Mark Schira (mark@ski.org) wrote it.
% 061407: ARW Modified it to use the new pipeline for processing data from
% the NIC.
% Found scans will contain a sparse list of processed scans that can be
% found in ./Reconned.
% The functional data will be in _4d.hdr files. The inplanes will be in a
% _3d file. There may be other data sets there (T1s etc taken during the
% same session).
% First we will identify the inplane scans. Then we will do the
% functionals.

mrGlobals
HOMEDIREC = pwd;

doCrop=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the inplane anatomies
% We can read the volume siye from the header. Maybe other info as well.
% The inplane data should be in the dataPath named anat.hdr/img

InplaneDir=0;
% find the InplaneDir
for thisScan=1:length(procScans)
    if (~isempty(procScans{thisScan}))
    if strcmp(procScans{thisScan}.Action,'Inplanes')
        InplaneDir=thisScan; %take the last if several are set
    end
    end

end
%get the anatomy
%anatDir=[HOMEDIREC,filesep,procScans{InplaneDir}.DirName];
if (~InplaneDir)
    % Prompt for an analyze file that serves as an inplane ref
    [fname,pathname]=uigetfile('*.nii','Pick Analyze/Nifti format 3d file');
    fname=fullfile(pathname,fname);
else
    fname=[procScans{InplaneDir}.Filename3d,'.nii'];
end

[anat, inplanes] = GetAnatomyFromAnalyze(fname);
% %since cropping is anoying
inplanes.cropSize=inplanes.fullSize;
inplanes.crop=cat(1,[1 1],inplanes.fullSize);

%AnatomySize=size(anat);
% Save anat
anatFile = fullfile(HOMEDIREC, 'Inplane', 'anat.mat');
save(anatFile, 'anat', 'inplanes');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the functionals
functScans=[];
for thisScan=1:length(procScans)
    if (~isempty(procScans{thisScan}))
        if strcmp(procScans{thisScan}.Action,'Functional')
            functScans=cat(2,functScans,thisScan); %take the last if several are set

        end
    end

end

% put information about the functionals here. I tried to make this independent of the DICOM format but in the end I needed to get at least the TR from that header.
% I know that the individual analyze vols should have the TR as the 4th
% dimension. Something to work on...

for thisScan=1:length(functScans)
    scanParams(thisScan).hdr=spm_vol([procScans{functScans(thisScan)}.Filename4d,'.nii']);
    scanParams(thisScan).imatrix=spm_imatrix(scanParams(thisScan).hdr(1).mat); % imatrix decomposes the analyze (nifti) header into its components.
    scanParams(thisScan).fullSize=scanParams(thisScan).hdr(1).dim; %#ok<AGROW>

    scanParams(thisScan).junkFirstFrames=procScans{functScans(thisScan)}.SkipVols;
    scanParams(thisScan).nFrames=procScans{functScans(thisScan)}.Volumes-procScans{functScans(thisScan)}.SkipVols;
% We don't crop these days. Perhaps we could get rid of that crop param
% altogether?
    scanParams(thisScan).cropSize=scanParams(thisScan).fullSize(1:2);
    scanParams(thisScan).crop=cat(1,[1 1],scanParams(1).fullSize(1:2));
    scanParams(thisScan).TR=procScans{functScans(thisScan)}.hdr{1}.RepetitionTime;
    scanParams(thisScan).PfileName=procScans{functScans(thisScan)}.hdr{1}.SeriesNumber;
    scanParams(thisScan).framePeriod=scanParams(thisScan).TR/1000; % this is in seconds
    scanParams(thisScan).slices=1:procScans{functScans(thisScan)}.Slices;
	 scanParams(thisScan).effectiveResolution = abs(scanParams(thisScan).imatrix(7:9));		% voxel dimensions (mm)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create/load mrSESSION and dataTYPES, modify them, and save

% If mrSESSION already exits, load it.
sessionFile = fullfile(HOMEDIREC, 'mrSESSION.mat');
if exist(sessionFile, 'file')
    loadSession;
    % if docrop, make sure that the mrSESSION is up-to-date
    if doCrop
        mrSESSION.inplanes = inplanes;
        mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams);
        saveSession;
    end
end

% If we don't yet have a session structure, make a new one.
if (~exist('mrSESSION','var'))
    mrSESSION = CreateNewSession(HOMEDIREC, inplanes, mrLoadRetVERSION);
end

if isempty(mrSESSION)
    mrSESSION = CreateNewSession(HOMEDIREC, inplanes, mrLoadRetVERSION);
end

% Update mrSESSION.functionals with scanParams corresponding to any new Pfiles.
% Set mrSESSION.functionals(:).crop & cropSize fields
mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams,0);

% this lines replace then EditSession command. One gui-window less to bother
mrSESSION.description=SessionName;
mrSESSION.subject=subjectName;

saveSession;

% Create/edit dataTYPES
if ieNotDefined('dataTYPES')
    dataTYPES = CreateNewDataTypes(mrSESSION);
else
    dataTYPES = UpdateDataTypes(dataTYPES,mrSESSION);
end

% Save any changes that may have been made to mrSESSION & dataTYPES
saveSession;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract time series
mrVista('inplane');
v=getSelectedInplane;

for thisScan=1:length(functScans)
    analyzeNum=functScans(thisScan);
    inFile=procScans{analyzeNum}.Filename4d; % This has a trailing .hdr part
    %make sure there is a .hdr at the end...
    [pathstr, name, ext, versn]=fileparts(inFile);
    inFile=[pathstr,filesep,name,'.hdr'];
    disp(inFile)
    volsToSkip=mrSESSION.functionals(thisScan).junkFirstFrames;
    scaleFact=[1];
    v=read4dAnalyzeToTseries(v,inFile,thisScan,volsToSkip,0,1,0); % No slice flip, 3 90 deg rotations
end

for thisScan=1:length(dataTYPES(1).blockedAnalysisParams)
    dataTYPES(1).blockedAnalysisParams(thisScan).temporalNormalization=0;
    dataTYPES(1).blockedAnalysisParams(thisScan).nCycles=procScans{functScans(thisScan)}.Cycles;
    dataTYPES(1).eventAnalysisParams(thisScan).eventAnalysis=0;
end
saveSession;

status=1;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Copyright (C) 2016-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun fresnels (@var{x})
%% Numerical fresnels function.
%%
%% Example:
%% @example
%% @group
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% fresnels (1.1)
%%   @result{} ans = 0.5365
%% @end group
%% @end example
%%
%% @strong{Note} this function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{fresnels}.
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @seealso{@@sym/fresnels}
%% @end defun

function y = fresnels (x)
  if (nargin ~= 1)
    print_usage ();
  end
  cmd = { 'L = _ins[0]'
          'A = [complex(mpmath.fresnels(x)) for x in L]'
          'return A,' };
  c = pycall_sympy__ (cmd, num2cell (x(:)));
  y = reshape (cell2mat (c), size (x));
end

%!error fresnels (1, 2)

%!test
%! x = 1.1;
%! y = sym(11)/10;
%! A = fresnels (x);
%! B = double (fresnels (y));
%! assert (A, B, -4*eps);

%!test
%! y = [2 3 sym(pi); exp(sym(1)) 5 6];
%! x = double (y);
%! A = fresnels (x);
%! B = double (fresnels (y));
%! assert (A, B, -4*eps);

%!test
%! % maple:
%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
%! % > for a in A do evalf(FresnelS(a)) end do;
%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
%! A = [  36.725464883991438430 + 15.587751104404587334*1i
%!        0.12213736710980573217e13 - 0.47688568479924574726e12*1i
%!        0.49681690114783755327
%!       -0.46816997858488224040*1i
%!       -0.52344169596561937423e-12 + 0.15707439669173367248e-13*1i
%!        0.75738824160998910399e24 + 0.15391592966931193100e26*1i ];
%! B = fresnels (x);
%! assert (A, B, -eps)

%!assert (fresnels (inf), 0.5)
%!assert (fresnels (-inf), -0.5)

repo: octsympy


text:
function r = and(p,q)
% this function compares the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVector')&&isa(q, 'CellVector'))
    r=p;
    r.xvalue = p.xvalue&q.xvalue;
    r.yvalue = p.yvalue&q.yvalue;
    r.zvalue = p.zvalue&q.zvalue;
elseif isa(p, 'CellVector')
    r=p;
    r.xvalue = p.xvalue&q;
    r.yvalue = p.yvalue&q;
    r.zvalue = p.zvalue&q;
else
    r=q;
    r.xvalue = p&q.xvalue;
    r.yvalue = p&q.yvalue;
    r.zvalue = p&q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function cycle_floyd_test03 ( )

%*****************************************************************************80
%
%% CYCLE_FLOYD_TEST03 tests CYCLE_FLOYD for F3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLE_FLOYD_TEST03\n' );
  fprintf ( 1, '  Test CYCLE_FLOYD for F3().\n' );
  fprintf ( 1, '  f3(i) = mod ( 123 * i + 456, 100000 ).\n' );

  x0 = 789;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Starting argument X0 = %d\n', x0 );

  [ lam, mu ] = cycle_floyd ( @f3, x0 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported cycle length is %d\n', lam );
  fprintf ( 1, '  Expected value is 50000\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported distance to first cycle element is %d\n', mu );
  fprintf ( 1, '  Expected value is 0\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bvalues=scd_scheme_bvalue(scheme)
gyro = 42.57; % kHz/mT
scheme(:,8) = gyro*scheme(:,4).*scheme(:,6); % um-1
bvalues=(2*pi*scheme(:,8)).^2.*(scheme(:,5)-scheme(:,6)/3);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [CD CLmax CDv qratio]= wingpolar(CL,AC,tState)
% inputs:
%   AC: see documentation
%   tState: see documentation
%   CL: wing lift coefficient; L = CL*q*Sref_wing
% outputs:
%   CD: wing drag coefficient; D = CD*q*Sref_wing
%   CLmax: WING maximum lift coefficient
%   CDv: drag coefficient of the wing due to the vertical downwash
%   qratio: ratio of downwash q at wing to 1/2*rho*(v_induced)^2

if AC.Wing.S
    e = .9;
    CD0 = .009;
    CD = CD0 + CL.^2./(pi*e*AC.Wing.AR);
else
    CD = .009;
end

if nargout >1
CLmax = 1.2; %this could be improved by adjusting CLmax to the AR
CDv = 1.3; %drag coefficient of the wing to the vertical downwash
qratio = .6; %ratio of downwash q at wing to 1/2*rho*(v_induced)^2
end

author: Sable

repo: mcbench-benchmarks


text:
function LabBookEmail(varargin)
%"LabbookEmail"
%   GUI to handle emailing CERR Labbook image captures.
%
%JRA 5/31/05
%
%Usage:
%   LabBookEmail('init', image, annotation);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

screenSize = get(0,'ScreenSize');
y = 300;
x = 519;
units = 'normalized';

%If no arguments default to init.
if(nargin == 0)
    varargin{1} = 'init';
end

%Command for callback switch statement.
command = varargin{1};

hFig = findobj('tag', 'CERR_LabBookEmail');

%If asking for INIT, and already initialized, return.
if ~isempty(hFig) & strcmpi(command, 'init');
    figure(hFig);
    return;
end

switch upper(command)
    case 'INIT'
        hFig = figure('visible', 'on', 'Name','Email Dialogue', 'units', 'pixels', 'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'CERR_LabBookEmail');
        uicontrol('units',units,'Position',[.01 .03 .54 .94],'style', 'frame');
        uicontrol('units',units,'Position',[.56 .82 .43 .15],'style', 'frame');
        uicontrol('units',units,'Position',[.56 .03 .43 .76],'style', 'frame');

        uicontrol('units',units,'Position',[.03 .9 .5 .05],'String','To:', 'style', 'text', 'HorizontalAlignment', 'left');
        uicontrol('units',units,'Position',[.03 .76 .5 .14],'String','', 'style', 'edit','HorizontalAlignment', 'left', 'Tag', 'LabBookEmailToField');

        uicontrol('units',units,'Position',[.03 .69 .5 .05],'String','Subject:', 'style', 'text', 'HorizontalAlignment', 'left');
        uicontrol('units',units,'Position',[.03 .62 .5 .07],'String','', 'style', 'edit','HorizontalAlignment', 'left', 'Tag', 'LabBookEmailSubjectField');

        uicontrol('units',units,'Position',[.03 .55 .5 .05],'String','Body:', 'style', 'text', 'HorizontalAlignment', 'left');
        uicontrol('units',units,'Position',[.03 .05 .5 .5],'String','', 'style', 'edit','HorizontalAlignment', 'left', 'Tag', 'LabBookEmailBodyField');

        uicontrol('units',units,'Position',[.59 .85 .17 .08],'String','Send', 'callback','LabBookEmail(''SEND'');', 'Tag', 'LabBookEmailSend');
        uicontrol('units',units,'Position',[.8 .85 .17 .08],'String','Cancel', 'callback','LabBookEmail(''CANCEL'');', 'Tag', 'LabBookEmailCancel');

        addresses = getEmailAddresses;
        optEmails = addresses.optionsEmails;
        cachedEmails = addresses.cachedEmails;

        uicontrol('units',units,'Position',[.57 .04 .41 .74],'style', 'listbox', 'string', {optEmails{:}, cachedEmails{:}}, 'callback', 'LabBookEmail(''ADDRESSCLICKED'')');

    case 'MESSAGE'
        saveFile = varargin{2};
        annotation = varargin{3};

        subjectField = findobj(hFig, 'Tag', 'LabBookEmailSubjectField');
        bodyField    = findobj(hFig, 'Tag', 'LabBookEmailBodyField');

        [pathstr,name,ext] = fileparts(saveFile);

        set(subjectField, 'String', [name ext]);
        set(bodyField, 'string', annotation);

        set(hFig, 'userdata', saveFile);

    case 'ADDRESSCLICKED'
        emails = get(gcbo, 'string');
        value  = get(gcbo, 'value');
        emailToAdd = emails{value};

        toBox = findobj(hFig, 'Tag', 'LabBookEmailToField');

        oldEmails = get(toBox, 'string');

        if ~isempty(oldEmails)
            set(toBox, 'string', [oldEmails emailToAdd ';']);
        else
            set(toBox, 'string', [emailToAdd ';']);
        end

    case 'CANCEL'
        delete(hFig);

    case 'SEND'
        optS = opts4Exe([getCERRPath,'CERROptions.json']);
        toField = findobj(hFig, 'Tag', 'LabBookEmailToField');
        subjectField = findobj(hFig, 'Tag', 'LabBookEmailSubjectField');
        bodyField    = findobj(hFig, 'Tag', 'LabBookEmailBodyField');

        subject = get(subjectField, 'string');
        toAddys = get(toField, 'string');
        saveFile = get(hFig, 'userdata');
        body = get(bodyField, 'string');

        toAddys = parseEmailString(toAddys);

        delete(hFig);
        pause(.1);
        sendmail(toAddys, subject, body, {saveFile});
        setEmailAddresses(toAddys);

end

function addresses = getEmailAddresses
%"getEmailAddresses"
%   Returns a list of stored and cached email addresses.

	optS = CERROptions;
	optionsEmails = optS.emailAddresses;

	try
        cachedEmails  = load(fullfile(getCERRPath, 'Labbook', 'cachedEmails'));
        cachedEmails = cachedEmails.cachedEmails;
	catch
        cachedEmails = {};
	end

	cachedEmails = setdiff(cachedEmails, optionsEmails);
	addresses.optionsEmails = optionsEmails;
	addresses.cachedEmails  = cachedEmails;

function setEmailAddresses(newAddresses)
%"setEmailAddresses"
%   Adds the passed addresses to the email cache.

	optS = CERROptions;
	optionsEmails = optS.emailAddresses;

	try
        cachedEmails  = load(fullfile(getCERRPath, 'Labbook', 'cachedEmails'));
	catch
        cachedEmails = [];
	end

	cachedEmails = union(cachedEmails, newAddresses);
	cachedEmails = setdiff(cachedEmails, optionsEmails);
    saveOpt = getSaveInfo;
    try
        if ~isempty(saveOpt);
        	save(fullfile(getCERRPath, 'Labbook', 'cachedEmails'), 'cachedEmails', saveOpt);
        else
        	save(fullfile(getCERRPath, 'Labbook', 'cachedEmails'), 'cachedEmails');
        end
    catch
        warning('Could not save cached email addresses.');
    end

function strings = parseEmailString(emailString)
%"parseEmailString"
%   Convert from a string of emails seperated by ; , and blanks to a cell
%   array of email addresses.

dividers = [0 length(emailString)+1];
dividers = union(dividers, strfind(emailString, ','));
dividers = union(dividers, strfind(emailString, ' '));
dividers = union(dividers, strfind(emailString, ';'));
dividers = unique(dividers);

dividerEnd = dividers(2:end) - 1;
dividers = dividers(1:end-1);

for i=1:length(dividers)
    strings{i} = emailString(dividers(i)+1:dividerEnd(i));
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Reconstruction of image from Laplacian pyramid
%
% Arguments:
%   pyramid 'pyr', as generated by function 'laplacian_pyramid'
%
% tom.mertens@gmail.com, August 2007
%
%
% More information:
%   'The Laplacian Pyramid as a Compact Image Code'
%   Burt, P., and Adelson, E. H.,
%   IEEE Transactions on Communication, COM-31:532-540 (1983).
%

function R = reconstruct_laplacian_pyramid_(pyr)

r = size(pyr{1},1);
c = size(pyr{1},2);
nlev = length(pyr);

% start with low pass residual
R = pyr{nlev};
filter = pyramid_filter;
for l = nlev - 1 : -1 : 1
    % upsample, and add to current level
    odd = 2*size(R) - size(pyr{l});
    R = pyr{l} + upsample_(R,odd,filter);
end

author: mahmoudnafifi


text:
function [W,B] = updateWeight(obj,objhat,N)
% Update the weight vectors

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Tomoaki Takagi

    fmin = min(obj,[],1);
    fmax = max(obj,[],1);
    obj = unique([normalize(obj,'range');objhat],'rows');

    %% Select the representative objective vectors
    LpNormD = pdist2(obj,obj,'minkowski',0.5);
    Choose  = false(1,size(obj,1));
    % Select the extreme objective vectors
    [~,extreme]     = min(pdist2(obj,eye(size(obj,2)),'cosine'),[],1);
    Choose(extreme) = true;
    % Greedy inclusion distance-based subset slection
    while sum(Choose) < N
        Remain   = find(~Choose);
        [~, rho] = max(min(LpNormD(Remain,Choose),[],2));
        Choose(Remain(rho)) = true;
    end
    obj = obj(Choose,:);

    %% Update the weight vectors
    W = obj.*repmat(fmax-fmin,N,1);
    B = pdist2(W,W);
    [~,B] = sort(B,2);
    B = B(:,1:ceil(N/10));
    W = W./vecnorm(W,1,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ n_data, n, x, fx ] = t_polynomial_values ( n_data )

%*****************************************************************************80
%
%% T_POLYNOMIAL_VALUES returns values of the Chebyshev polynomial T(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA
%    from the previous call.
%
%    Output, integer N_DATA, the index of the test data.
%
%    Output, integer N, the order of the function.
%
%    Output, real X, the point where the function is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 13;
  fx_vec = [ ...
     1.0000000000E+00,  0.8000000000E+00,  0.2800000000E+00, ...
    -0.3520000000E+00, -0.8432000000E+00, -0.9971200000E+00, ...
    -0.7521920000E+00, -0.2063872000E+00,  0.4219724800E+00, ...
     0.8815431680E+00,  0.9884965888E+00,  0.7000513741E+00, ...
     0.1315856097E+00 ];
  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10, 11, ...
    12 ];
  x_vec = [ ...
    0.8E+00,  0.8E+00,  0.8E+00, ...
    0.8E+00,  0.8E+00,  0.8E+00, ...
    0.8E+00,  0.8E+00,  0.8E+00, ...
    0.8E+00,  0.8E+00,  0.8E+00, ...
    0.8E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0E+00;
    fx = 0.0E+00;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function out = jaccoeffs(f, n, a, b)
%JACCOEFFS   Jacobi polynomial coefficients of a CHEBFUN.
%   A = JACCOEFFS(F, N, ALPHA, BETA) returns the first N+1 coefficients in the
%   Jacobi series expansion of the CHEBFUN F, so that such that F approximately
%   equals A(1) J_0(x) + ... + A(N+1) J_N(x) where J_N(x) denotes the N-th
%   Jacobi polynomial with parameters ALPHA and BETA. A is a column vector.
%
%   If F is smooth (i.e., numel(f.funs) == 1), then A = JACCOEFFS(F, ALPHA,
%   BETA) will assume that N = length(F).
%
%   JACCOEFFS does not support quasimatrices.
%
% See also CHEBCOEFFS, LEGCOEFFS.

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

if ( (numel(f(1).funs) == 1) && (nargin < 4) )
    b = a;
    a = n;
    n = length(f);
elseif ( isempty(n) )
    n = length(f);
end

if ( numel(f) > 1 )
    out = zeros(n, numel(f));
    f = cheb2cell(f);
    for k = 1:numel(f)
        out(:,k) = jaccoeffs(f{k}, n, a, b);
    end
    return
end

%%
% Special cases:
if ( a == 0 && b == 0 )
    out = legcoeffs(f, n);
    return
elseif ( a == -.5 && b == -.5 )
    out = chebcoeffs(f, n, 'kind', 1);
    nn = 0:(n-2);
    scl = repmat([1 ; cumprod((nn'+.5)./(nn'+1))], 1, size(out, 2));
    out = out./scl;
    return
elseif ( a == .5 && b == .5 )
    out = chebcoeffs(f, n, 'kind', 2);
    nn = 0:(n-2);
    scl = repmat([1 ; cumprod((nn'+1.5)./(nn'+1))], 1, size(out, 2));
    out = (1:n)'.*out./scl;
    return
end

%%
if ( numel(f.funs) == 1 )

    % Compute Jacobi coefficients of the single fun.
    out = jaccoeffs(f.funs{1}, n, a, b);

else

    % Jacobi-Vandermonde matrix:
    Enorm = jacpoly(0:n-1, a, b, f.domain);

    % Compute the projection (with correct scaling):
    out = Enorm \ f;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [x, output] = findExtremePool(fbaModel, obj, printLevel)
% Finds an extreme ray in the left nullspace of the stoichiometric matrix
%
% USAGE:
%
%    [x, output] = findExtremePool(fbaModel, obj, printLevel)
%
% INPUT:
%    fbaModel:       FBA type model
%
% OPTIONAL INPUT:
%    obj:            default = random vector with size depending on `fbaModel.S`
%    printLevel:     argument for `solveCobraLP` function, default = 0
%
% OUTPUTS:
%    x:              `x = output.full`
%    output:         `output = solveCobraLP(LPProblem)`
A    = fbaModel.S';
[n, m] = size(A);

if nargin < 2
    obj = rand(m,1);
end
if nargin < 3
    printLevel = 0;
end

LPProblem.A=sparse([A; ones(1,m)]);
LPProblem.b=[zeros(n,1); 1];
LPProblem.c=obj;
LPProblem.lb=-100*ones(size(LPProblem.A,2),1);
LPProblem.ub= 100*ones(size(LPProblem.A,2),1);
LPProblem.osense=-1;
LPProblem.csense(1:size(LPProblem.A,1),1)='E';
output = solveCobraLP(LPProblem, 'printLevel', printLevel);
x=output.full;

author: opencobra

repo: cobratoolbox


text:
function DPMpcOneBeam9(leak, spectrum_File, planC_File, nhist, OutputError, whichBeam, PBMaxWidth, gradsense, MCsolver, saveIM, sourceModel, doseToWater, fillWater, useWedge, inputPB, Softening, batch)
% written by KZ, JC
% Calculate doseV dose3D directly from DPM, without compressing it.
% Use generateDPMdose.m instead of DPMInfluenceJing.m
% Calculate IM struct, call DPM, get
%
% Usage: IIMCalc = calculateBeamIM(0.018, 'DPM_10beams_planC', 1, 1, 5)
%   leak = 0.018
%   spectrum_File = file name of photon spectrum. eg. photon6MV.spectrum
%   planC_File = file name contained planC, eg. 'DPM_10beams_planC'
%   indexBeam = Beam Index: 1, 2, ...
%   imin = Beamlet index
%   imax = Beamlet index
%
%   eval(['load Beam',num2str(indexBeam),'_w']);
%   where, Beam1_w.mat (Beam2_w.mat) need to be in the current directory. It's the weight of
%   the beamlets for each beam.
%
% JC Feb 27 06
% Call in DPMInfluence.m; use sum(10*clock) as the random seeds input to DPM.
%
% JC. Aug 10, 2005, Need to make planC read from the disk, since in the
% stand-alone mode, no matlab run is expected.
%
% JC Dec 18 06
% Add 'saveIM' input,
%   If it's larger than 0, say 1, then save IM
%   If it's 0, don't save IM
%
% LM: JC Jan 26 2007
% Include non-zero couchAngle
% The desination coordinates should be the patient support
% system.
% The previous assumption is that the "fixed system" is the same as
% the "Patient support system".
% test
%
% LM: JC Mar 03, 2007
% Score dose to water flag. Only valid for VMC++ for now. "doseToWater"
% Add "fillWater" flag. It 1, fill the skin strucutere with all water
% density, i.e. CT value = 1024.
%
% LM: JC Mar 29, 2007
% get rid of 'stateS' arguments to 'generateVMCdose.m', 'generateDPMdose.m'
% get rid of the call for 'IMRTP_plancheck'
% add 'useWedge' flag. default = 0; when useWedge == 1, use Elekta
% universal wedge.
% only implemented for DPM, openField == 0, sourceModel == 0;
%
% open planC_File
% In matlab release 14, "eval" can be compiled.
Error = []; %Initialize Error to NULL.
load(planC_File);
%eval(['load ',planC_File]);
% This planC_File should contains planC and stateS
%leak = 0.018; %leakage in per *100
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ischar(leak)
    leak = str2num(leak)
end

if ischar(nhist)
    nhist = str2num(nhist)
end

if ischar(OutputError)
    OutputError = str2num(OutputError)
end

if ischar(whichBeam)
    whichBeam = str2num(whichBeam)
end

if ischar(PBMaxWidth)
    PBMaxWidth = str2num(PBMaxWidth)
end

if ischar(gradsense)
    gradsense = str2num(gradsense)
end

% If MCsolver == 1, use DPM, default
% If MCsolver == 2, use VMC++
% other, throw out error
if ischar(MCsolver)
    MCsolver = str2num(MCsolver)
end

% If saveIM == 0, calculate doseV, do not calc/save IM
% If saveIM ~= 0, calc/save IM, only valid for DPM for now, ie. MCsolver == 1.
if ischar(saveIM)
    saveIM = str2num(saveIM)
end

% If sourceModel == 0, no source model, only point source
% If sourceModel == 1, source model version 1: point source (Primary) + FF
% + electron contamination.
% Expandable to other source models for the future development.

if ischar(sourceModel)
    sourceModel = str2num(sourceModel)
end

if ischar(doseToWater)
    doseToWater = str2num(doseToWater)
end

if ischar(fillWater)
    fillWater = str2num(fillWater)
end

if ischar(useWedge)
    useWedge = str2num(useWedge)
end

if ischar(inputPB)
    inputPB = str2num(inputPB)           %Flag, inputPB == 1, ask user to input .mat file storing PB info.
                                 % That file can be generated using
                                 % function generateOpenfiledPB.m
end
% Now, only work for DPM, planC{7} not empty.

if ischar(Softening)
    Softening = str2num(Softening)           %Flag,  == 1, use off-axis-softening in dpm calc.

end
% Now, only work for DPM, planC{7} not empty.

if ischar(batch)
    batch = str2num(batch)
end

currentDir = cd;
if (~isempty(planC{7}))
    % When there are multiple beam definition in planC, use the first BEAM
    % data as the default.

    %for indexBeam = 1 : planC{7}(1).FractionGroupSequence.Item_1.NumberOfBeams
    for indexBeam = whichBeam

            bs = planC{7}(1).BeamSequence.(['Item_' num2str(indexBeam)]);

            %Parse BlockData into PB information.
            input = bs.BlockSequence.Item_1.BlockData;
            input = reshape(input, 2, [])';

            [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y]=...
                BlockToPB(input, 10, 10);   %PB size here is in mm.
            w_field = ones(size(xPosV));

            gA = bs.ControlPointSequence.Item_1.GantryAngle;
            iC = bs.ControlPointSequence.Item_1.IsocenterPosition;
            IMBase.beams.gantryAngle = gA;
            % couchAngle is PatientSupportAngle, in DICOM, i.e. IEC format
            % (in deg)
            IMBase.beams.couchAngle = bs.ControlPointSequence.Item_1.PatientSupportAngle;
            % JC Apr. 10 2007,
            % Add Beam Limiting Deviceing Angle, i.e. collimator angle from
            % DICOM, in IED format. (in deg)
            % Not test for "Non-zero angle" yet.
            IMBase.beams.collimatorAngle = bs.ControlPointSequence.Item_1.BeamLimitingDeviceAngle

            try
            if ~isfield(planC{7}(1).PatientSetupSequence,(['Item_' num2str(indexBeam)]))
                position = {planC{7}(1).PatientSetupSequence.(['Item_' num2str(1)]).PatientPosition};
            else
                position = {planC{7}(1).PatientSetupSequence.(['Item_' num2str(indexBeam)]).PatientPosition};
            end
            catch
                position = 'HFS';
                disp('Use default patient posiont HFS')
            end

            if strcmpi(position, 'HFP')
                IMBase.beams(1).isocenter.x = iC(1)/10;
                IMBase.beams(1).isocenter.y = iC(2)/10;
                IMBase.beams(1).isocenter.z = -iC(3)/10;
            else
                IMBase.beams(1).isocenter.x = iC(1)/10;
                IMBase.beams(1).isocenter.y = -iC(2)/10;
                IMBase.beams(1).isocenter.z = -iC(3)/10;
            end

            IMBase.beams(1).isodistance = bs.SourceAxisDistance/10;

            IMBase.beams(1).beamEnergy = bs.ControlPointSequence.Item_1.NominalBeamEnergy;

            xPosV = xPosV/10;
            yPosV = yPosV/10;
            beamlet_delta_x = beamlet_delta_x/10;
            beamlet_delta_y = beamlet_delta_y/10;

            IMBase.goals.PBMargin = 0.5;
            IMBase.goals.structNum = 2;
            IMBase.params.xyDownsampleIndex = 1;
            IMBase.goals.isTarget(1) = 'y';
            IMBase.goals.xySampleRate = 2;
            IMBase.params.numCTSamplePts = 300;

            IMBase.beams(1).zRel = 0;
            IMBase.beams(1).xRel =  IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);
            IMBase.beams(1).yRel =  IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

            % JC In the above calc. for xRel and yRel, the assumption is the couchAngle
            % == 0. i.e. the "fixed system" is the same as the "patient support
            % system".
            % LM: JC Jan 26 2007
            % Include non-zero couchAngle
            % The desination coordinates should be the patient support
            % system.
            % The previous assumption is that the "fixed system" is the same as
            % the "Patient support system".
            % test

            %         patientVectorsM = zeros(size(RTOGVectorsM));
            patientxRel = cosdeg(IMBase.beams.couchAngle) * IMBase.beams(1).xRel- sindeg(IMBase.beams.couchAngle) * IMBase.beams(1).zRel;
            patientyRel = IMBase.beams(1).yRel;
            patientzRel = sindeg(IMBase.beams.couchAngle) * IMBase.beams(1).xRel + cosdeg(IMBase.beams.couchAngle) * IMBase.beams(1).zRel;
            IMBase.beams(1).xRel = patientxRel;
            IMBase.beams(1).yRel = patientyRel;
            IMBase.beams(1).zRel = patientzRel;
            clear patientxRel patientyRel patientzRel

            %IMBase.params.algorithm = 'VMC++';
            %IMCalc = IMRTP_MapCheck(IMBase);

            IMBase.params.algorithm = 'DPM';
            %      load PB_40x40.mat beamlet_delta_x beamlet_delta_y xPosV yPosV w_field
            %      load PB_15x15.mat beamlet_delta_x beamlet_delta_y xPosV yPosV w_field

            if inputPB         % get PB info from a specified .mat file.
                setappdata(0, 'usenativesystemdialogs', false)
                currentDir = cd;

                [FileName,path] = uigetfile('*.mat','Select MAT file containing beamlet_delta_x beamlet_delta_y xPosV yPosV w_field');

                if path == 0
                    errordlg('File Should exist');
                    error('File Should exist');
                end

                cd(path);
                load (FileName, 'beamlet_delta_x', 'beamlet_delta_y', 'xPosV', 'yPosV', 'w_field')
                cd(currentDir);
            end

            disp('Number of pencile beams are:')
            length(xPosV)

            sourceS = IMBase.beams(1);
            [RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
                    getPBRays(xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, sourceS);
            IMCalc = IMBase;
            IMCalc.beams(1).RTOGPBVectorsM_MC = RTOGPBVectorsM_MC;
            IMCalc.beams(1).RTOGPBVectorsM    = RTOGPBVectorsM;
            IMCalc.beams(1).xPBPosV           = xPBPosV;
            IMCalc.beams(1).yPBPosV           = yPBPosV;
            IMCalc.beams(1).rowPBV            = rowPBV;
            IMCalc.beams(1).colPBV            = colPBV;
            % No need of this field
            IMCalc.beams(1).CTTraceS          = struct;
            IMCalc.beams(1).beamletDelta_x    = beamlet_delta_x;
            IMCalc.beams(1).beamletDelta_y    = beamlet_delta_y;

            %RTOG positions of sources
            IMCalc.beams(1).x = IMCalc.beams(1).xRel + IMCalc.beams(1).isocenter.x;
            IMCalc.beams(1).y = IMCalc.beams(1).yRel + IMCalc.beams(1).isocenter.y;
            IMCalc.beams(1).z = IMCalc.beams(1).zRel + IMCalc.beams(1).isocenter.z;
            %IMCalc = IMRTP_MapCheck(IMBase, planC, stateS, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
            % JC. 11 Aug, 05 Added more arguments to IMRTP_MapCheck,

            % JC. Aug. 3, 2005
            % Add fields to IMCalc, necessary for generateDPMInfluence.m
            IMCalc.beams.collimatorAngle = 0;
            IMCalc.params.ScatterMethod = 'exponential';
            % Use 0.002 instead the original default 0.01 as threshold value.
            % Change to 0.004, since the stored IM struct is very large.
            IMCalc.params.Scatter.Threshold = 0.004;
            IMCalc.params.Scatter.RandomStep = 30;
            IMCalc.beamlets = struct;

            %JC Add/porpulate VMC field in IM struct.
            IMCalc.params.VMC.NumParticles = nhist;
            IMCalc.params.VMC.NumBatches = 10;
            IMCalc.params.VMC.scoreDoseToWater = 'No';
            IMCalc.params.VMC.monoEnergy = 0;   %?
            IMCalc.params.VMC.repeatHistory = 0.2510;
            IMCalc.params.VMC.splitPhotons = 'Yes';
            IMCalc.params.VMC.photonSplitFactor = -40;
            IMCalc.params.VMC.base = 2;;
            IMCalc.params.VMC.dimension = 60;    %?
            IMCalc.params.VMC.skip = 0;  %? 1?
            IMCalc.params.VMC.includeError = 'No';
            IMCalc.params.VMC.spectrum = spectrum_File;
            IMCalc.beams.beamModality = 'photons';

            if doseToWater
            IMCalc.params.VMC.scoreDoseToWater = 'Yes';
            disp('VMC++ score dose to Water.');
            end

            % JC Jun 06 2006
            % Still need to fill in the structure of IMCalc, but don't need put
            % it back(?), since it's not to be used to calculate dose.
            % Will also need to use 'DoseScale' inside generateDPMdose.m,
            % since DPM output dose on a per partical base.
            % Need to use w_field as a way of telling DPM how to add beamlets
            % up.

            % JC Aug 30. 2006
            if (MCsolver == 1)

                [IMCalc doseV indV numberParticles] = generateDPMdose7(IMCalc,spectrum_File,nhist, OutputError,planC, indexBeam, w_field, saveIM, sourceModel, fillWater, useWedge, Softening);

            elseif (MCsolver == 2)
                % use VMC++
                [doseV indV] = generateVMCdose(IMCalc, planC, w_field, sourceModel, doseToWater, fillWater,saveIM);
            else
                error('INVALID input: MCsolver == 1, DPM; MCsolver == 2, VMC++')
            end

    end

else %If it's a conventional plan, planC{7} is empty

    indexBeam = whichBeam;
    if (length(planC{6}(indexBeam).file) > 5)
        [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, beamGeometry]=...
            planCToPB(planC, indexBeam, 1, 1);
    else
        % The field is shaped by Jaws, not MLC.
        [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, beamGeometry]=...
            planCToPBJaws(planC, indexBeam, 1, 1);
    end

    IMBase.beams(1).isocenter.x = beamGeometry.isocenter(1);
    IMBase.beams(1).isocenter.y = beamGeometry.isocenter(2);
    IMBase.beams(1).isocenter.z = beamGeometry.isocenter(3);
    % RTOG file has the same couchAngle definition as in IEC1217.
    IMBase.beams(1).couchAngle = beamGeometry.couchAngle;
    % RTOG file has different gantryAngle definition from IEC1217.
    % Use IEC 1217 coordinates.
    IMBase.beams(1).gantryAngle = - beamGeometry.gantryAngle;
    IMBase.beams(1).isodistance = beamGeometry.isoDistance;

    IMBase.goals.PBMargin = 0.5;
    IMBase.goals.structNum = 2;
    IMBase.params.xyDownsampleIndex = 1;
    IMBase.goals.isTarget(1) = 'y';
    IMBase.goals.xySampleRate = 1;
    IMBase.params.numCTSamplePts = 300;
    IMBase.beams(1).beamEnergy = planC{6}(1).beamEnergyMeV;

    IMBase.beams(1).zRel = 0;
    IMBase.beams(1).xRel =  IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);
    IMBase.beams(1).yRel =  IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

    % JC In the above calc. for xRel and yRel, the assumption is the couchAngle
    % == 0. i.e. the "fixed system" is the same as the "patient support
    % system".
    % LM: JC Jan 26 2007
    % Include non-zero couchAngle
    % The desination coordinates should be the patient support
    % system.
    % The previous assumption is that the "fixed system" is the same as
    % the "Patient support system".

    patientxRel = cosdeg(IMBase.beams.couchAngle) * IMBase.beams(1).xRel- sindeg(IMBase.beams.couchAngle) * IMBase.beams(1).zRel;
    patientyRel = IMBase.beams(1).yRel;
    patientzRel = sindeg(IMBase.beams.couchAngle) * IMBase.beams(1).xRel + cosdeg(IMBase.beams.couchAngle) * IMBase.beams(1).zRel;
    IMBase.beams(1).xRel = patientxRel;
    IMBase.beams(1).yRel = patientyRel;
    IMBase.beams(1).zRel = patientzRel;
    clear patientxRel patientyRel patientzRel

    IMBase.beams(1).collimatorAngle = beamGeometry.collimatorAngle;

    disp('Number of pencile beams are:')
            length(xPosV)

    % JC Mar 29, 2007
    % Repace 'IMRTP_plancheck' by the following command
            sourceS = IMBase.beams(1);
            [RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
                    getPBRays(xPosV, - yPosV, beamlet_delta_x, beamlet_delta_y, sourceS);
    % Note: In the above line, use '- yPosV' instead of 'yPosV'. Since IEC and RTOG have different definition about yPosV.
            IMCalc = IMBase;
            IMCalc.beams(1).RTOGPBVectorsM_MC = RTOGPBVectorsM_MC;
            IMCalc.beams(1).RTOGPBVectorsM    = RTOGPBVectorsM;
            IMCalc.beams(1).xPBPosV           = xPBPosV;
            IMCalc.beams(1).yPBPosV           = yPBPosV;
            IMCalc.beams(1).rowPBV            = rowPBV;
            IMCalc.beams(1).colPBV            = colPBV;
            % No need of this field
            IMCalc.beams(1).CTTraceS          = struct;
            IMCalc.beams(1).beamletDelta_x    = beamlet_delta_x;
            IMCalc.beams(1).beamletDelta_y    = beamlet_delta_y;

            %RTOG positions of sources
            IMCalc.beams(1).x = IMCalc.beams(1).xRel + IMCalc.beams(1).isocenter.x;
            IMCalc.beams(1).y = IMCalc.beams(1).yRel + IMCalc.beams(1).isocenter.y;
            IMCalc.beams(1).z = IMCalc.beams(1).zRel + IMCalc.beams(1).isocenter.z;

    % JC Mar 29, 2007
    % Replace 'IMRTP_MapCheck' by above commands.
    % IMCalc = IMRTP_MapCheck(IMBase, planC, stateS, xPosV, - yPosV, beamlet_delta_x, beamlet_delta_y, gA);
    clear IMBase beamlet_delta_x beamlet_delta_y PlanC_File xPosV yPosV

    % Add fields to IMCalc, necessary for generateDPMInfluence.m
    IMBase.params.algorithm = 'DPM';
    IMCalc.params.ScatterMethod = 'exponential';
    IMCalc.params.Scatter.Threshold = 0.004;
    IMCalc.params.Scatter.RandomStep = 30;
    IMCalc.beamlets = struct;

    % For the conventional RT. No change in w_field
    w_field = ones(length(IMCalc.beams.beamletDelta_x),1);
    %[doseV indV] = generateDPMdose(IMCalc,spectrum_File,nhist, OutputError,planC,stateS, indexBeam, w_field);
    % Dec 18, 2006
    % Use the new modified DPM, with Source Model.

    %JC Add/porpulate VMC field in IM struct.
    IMCalc.params.VMC.NumParticles = nhist;
    IMCalc.params.VMC.NumBatches = 10;
    IMCalc.params.VMC.scoreDoseToWater = 'No';
    IMCalc.params.VMC.monoEnergy = 0;   %?
    IMCalc.params.VMC.repeatHistory = 0.2510;
    IMCalc.params.VMC.splitPhotons = 'Yes';
    IMCalc.params.VMC.photonSplitFactor = -40;
    IMCalc.params.VMC.base = 2;;
    IMCalc.params.VMC.dimension = 60;    %?
    IMCalc.params.VMC.skip = 0;  %? 1?
    IMCalc.params.VMC.includeError = 'No';
    IMCalc.params.VMC.spectrum = spectrum_File;
    IMCalc.beams.beamModality = 'photons';

    if doseToWater
            IMCalc.params.VMC.scoreDoseToWater = 'Yes';
            disp('VMC++ score dose to Water.');
     end

    if (MCsolver == 1)

        [IMCalc doseV indV numberParticles] = generateDPMdose7(IMCalc,spectrum_File,nhist, OutputError,planC, indexBeam, w_field, saveIM, sourceModel, fillWater, useWedge, Softening);

    elseif (MCsolver == 2)
        % use VMC++
        [doseV indV] = generateVMCdose(IMCalc, planC, w_field, sourceModel, doseToWater, fillWater,saveIM);
    else
        error('INVALID input: MCsolver == 1, DPM; MCsolver == 2, VMC++')
    end

    % ?? Right place to put end???
end

%%%%%%%%%%%%%%%
    siz = getUniformizedSize(planC);
    dose3D = zeros(siz);

    if (OutputError == 1)            % OutputError only valid for DPM,ie., MCsolver == 1
        doseV = doseV(:,1);
        Error_3D = zeros(siz);
        Error_3D(indV') = doseV(:,2);
    else
        disp('Do not output error')
        dose3D(indV') = doseV;
    end

    %        filename = ['doseV_',num2str(indexBeam),'_', num2str(nhist), '_', num2str(batch)];
    %        save(filename, 'doseV', 'indV');

    filename = ['dose3D_',num2str(indexBeam),'_', num2str(nhist), '_', num2str(batch)];
    save(filename, 'dose3D');

    % Only Output w_field for the first batch
    if(batch == 1)
        filename = ['w_field',num2str(indexBeam)];
        save(filename, 'w_field');%dose3D = (dose3D/max(dose3D(:)));
    end

if (saveIM == 1)
%saveIM  %Only valid for DPM, i.e. IM
    filename = ['IMCalc_Beam',num2str(indexBeam),'_', num2str(nhist), '_', num2str(batch)];
    save(filename, 'IMCalc');
end
clear IMCalc dose3D w_field

end


text:
function [retPath,retDist,cycleNodes]=BellmanFordAlg(adjMat,sourceIdx,destIdx)
%%BELLMANFORDALG Use the Bellman-Ford algorithm to find the shortest path
%                from a given source node through a graph represented using
%                an adjacency matrix. Negative edge costs are allowed and
%                if a negative cycle reachable from the source exists in
%                the graph, then that cycle will be returned (and no
%                minimum distance path will exist).
%
%INPUTS: adjMat  An adjacency matrix full of costs for the graph.
%                adjMat(i,j) is the cost of going from vertex i to vertex
%                j. If there is no connection between the vertices in the
%                graph, then set that element to infinity. Cycles
%                accessible from the source involving negative edge weights
%                are detected by the algorithm and mean that no solution
%                exists.
%     sourceIdx  The index of the node (vertex) from which the shortest
%                paths to other vertices is desired.
%        destIdx The index of the destination path. If this parameter is
%                provided, then the path of nodes from the source to the
%                destination is returned along with the path length.
%                Otherwise, information allowing one to reconstruct the
%                whole path is returned if this parameter is omtted.
%
%OUTPUTS: retPath If destIdx is given, this is the minimum cost sequence of
%                 nodes to get from sourceIdx to destIdx, including the end
%                 nodes; if there is no path to the node, an empty matrix
%                 is returned. Otherwise, this is prevNodes, such that one
%                 can recreate all shortest paths. Starting at index idx,
%                 the next node on the path in the reverse direction from
%                 destimation to source is retPath(idx). By looping, one
%                 can recreate the whole path. If there is no path to a
%                 node, then retPath(idx)=0. If a negative cycle exists
%                 in the path, then an empty matrix is returned.
%         retDist If destIdx is given, this is the shortest distance from
%                 the source to the destination. If destIdx is not given,
%                 then this is an array where each index holds the shortest
%                 distance from the source to that node. If a negative
%                 cycle exists in the graph, then an empty matrix is
%                 returned.
%      cycleNodes If no negative cost cycles exist in the graph, then this
%                 is an empty matrix. If any negative cycles exist, then
%                 this will be the nodes in order of the first negative
%                 cycle found. The last node in cycleNodes is connected
%                 back to the first node, finishing the loop. A cycle is a
%                 loop in the graph; a negative cost cycle is one where the
%                 sum of the costs in the cycle is negative.
%
%The algorithm implemented is that of Chapter 24.1 of [1], where the
%extraction of a negative cost cycle accessible from the source is added
%when a negative cost cycle is detected.
%
%Ad an example, consider the adjacency matric with two negative cost
%cycles:
% adjMat=[Inf,    -20,    Inf;
%         Inf,    Inf,    10,
%         5,      -15,    Inf];
% [retPath,retDist,cycleNodes]=BellmanFordAlg(adjMat,1,3);
%Looking for a path from 1 to 3, the algorithm would encounter the cycle.
%Thus, the return values would be [] and [] for retPath and retDist and
%cycleNodes=[3;2], because that is the first cycle found.
%
%REFERENCES:
%[1] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein,
%    Introduction to Algorithms, 2nd ed. Cambridge, MA: The MIT Press,
%    2001.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The number of vertices.
numNodes=length(adjMat);

%The shortest distance yet found from the source to each of the nodes.
dist=inf(numNodes,1);
%The distance from the source to itself is zero. This is the first node to
%be visited.
dist(sourceIdx)=0;
%The previous node in the optimal path from the given node to the source.
prevNodes=zeros(numNodes,1);

for i=1:(numNodes-1)
    for curNode=1:numNodes
        %We have to go over all edges from the current node. If
        %adjMat(curNode,destNode) is infinite, then no edge exists to destNode.
        for destNode=1:numNodes
            w=adjMat(curNode,destNode);%The edge weight
            if(~isfinite(w))
                continue;
            end

            %Perform the relaxation step.
            if(dist(curNode)+w<=dist(destNode))
                dist(destNode)=dist(curNode)+w;
                prevNodes(destNode)=curNode;
            end
        end
    end
end

%Now, we check for a negative weight cycle (reachable from the source). If
%a negative weight cycle is found, then return it. We have to look at all
%edges in the graph.
cycleStart=[];
for curNode=1:numNodes
    for destNode=1:numNodes
        w=adjMat(curNode,destNode);%The edge weight
        if(~isfinite(w))
            continue;
        end

        if(dist(curNode)+w<dist(destNode))
            cycleStart=destNode;
            prevNodes(destNode)=curNode;
            break;
        end
    end
end

%If a negative weight cycle is present, then trace back along it to get the
%cycle to return. However, a problem is that from the starting node, we
%might get lost in a cycle that does not return to the starting node. Thus,
%we have to keep track of when we visited each edge and declare a cycle
%once we revisit any edge, not necessarily the starting edge.

if(~isempty(cycleStart))
    %Allocate the maximum space possible for the cycle.
    cycleNodes=zeros(numNodes,1);
    nodeVisited=zeros(numNodes,1);

    curIdx=1;
    cycleNodes(curIdx)=cycleStart;
    nodeVisited(cycleStart)=curIdx;
    while(1)
        nextNode=prevNodes(cycleNodes(curIdx));
        %if we have seen this node before
        if(nodeVisited(nextNode)~=0)
            break;
        end
        cycleNodes(curIdx+1)=prevNodes(cycleNodes(curIdx));
        nodeVisited(nextNode)=curIdx+1;

        curIdx=curIdx+1;
    end

    %Shrink the list of cycle nodes down to the size of the number of
    %actual nodes in the cycle and reverse the order so that the nodes are
    %in the correct order.
    cycleNodes=flipud(cycleNodes(nodeVisited(nextNode):curIdx));

    %It cannot solve for the optimal path, because of the negative cycle.
    retPath=[];
    retDist=[];
    return;
else
    cycleNodes=[];
end

%If a destination index is given, the recreate the path, if there are no
%cycles.
if(nargin>2)
    %If there is no path to the desired node.
    if(~isfinite(dist(destIdx)))
        retPath=[];
        retDist=Inf;
        return;
    end

    path=zeros(numNodes,1);
    path(1)=destIdx;
    pathStep=1;
    while(path(pathStep)~=sourceIdx)
        path(pathStep+1)=prevNodes(path(pathStep));
        pathStep=pathStep+1;
    end
    retPath=flipud(path(1:pathStep));
    retDist=dist(destIdx);
else
    retPath=prevNodes;
    retDist=dist;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function coef=spreadfun(T)
%SPREADFUN  Spreading function of a matrix
%   Usage:  c=spreadfun(T);
%
%   `spreadfun(T)` computes the spreading function of the operator *T*,
%   represented as a matrix. The spreading function represent the operator *T*
%   as a weighted sum of time-frequency shifts. See the help text for
%   |spreadop| for the exact definition.
%
%   See also:  spreadop, tconv, spreadinv, spreadadj

complainif_argnonotinrange(nargin,1,1,mfilename);

if ndims(T)>2 || size(T,1)~=size(T,2)
    error('Input symbol T must be a square matrix.');
end;

L=size(T,1);

% The 'full' appearing on the next line is to guard the mex file.
coef=comp_col2diag(full(T));

coef=fft(coef)/L;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [P] = unmixP_NCM(X,E,Sigma,Parameters)
%% Input:
%        X - d-by-N HSI data, matrix
%        E - d-by-M endmember mean set, matrix
%        Sigma - d-by-d-by-M endmember covariance set, matrix
%  Output:
%        P - M-by-N proportion matrix
% Author: Alina Zare et al. Rewriten by Sheng Zou
% Department of Electrical and Computer Engineering, University of Florida
% 09/07/2017

%% Initialization of P
[d,N] = size(X);
M = size(E,2);
P = 1/M.*ones(M,N); % intialize proportion values to be 1/M

%Initialize all Likelihood and Prior Values
LogLikelihoodOld          = ComputeLogLikelihoodAll(X, E, P, Sigma, N);

for iteration = 2:Parameters.NumberIterations+1
    Y = randg(1, N, M) ;
    v = sum(Y,2);
    samples = (Y./repmat(v,[1,size(Y,2)]))';
    LogLikelihoodNew      = ComputeLogLikelihoodAll(X, E, samples, Sigma, N);
    Ratios                = exp(LogLikelihoodNew - LogLikelihoodOld);
    rands                 = rand(1,N);
    Vals                  = rands < Ratios;
    Vrep                  = repmat(Vals,M,1);
    P                     = samples.*Vrep  + P.*(1-Vrep);

    LogLikelihoodOld = (1-Vals).*LogLikelihoodOld + (Vals).*LogLikelihoodNew;
    if mod(iteration, round(Parameters.NumberIterations/10)) == 0
        disp(strcat('iteration = ',num2str(iteration)));
        disp(strcat('loglikelihood = ',num2str(sum(LogLikelihoodOld))));
    end
end

end

function [LogLikelihoodAll] = ComputeLogLikelihoodAll(X, E, P, cov, N)
% compute log likelihood of all points

term1 = zeros(1,N);
term2 = zeros(1,N);

for s = 1:size(E,2)
    statement(s)=isdiag(squeeze(cov(:,:,s))) && length(unique(diag(cov(:,:,s))))==1; % check if all the covariance matrices are diagonal and isotropic
end

if mean(statement) ==1 % all diagonal and isotropic
    for z = 1:size(P,1)
        a(z) = unique(diag(cov(:,:,z)));
    end
    term1 = -.5*size(X,1)*log(a*(P.^2));
    term2 = -.5*(sum((X-E*P).^2)./(a*(P.^2)));
else

    parfor t = 1:N
        P3D = zeros(1,1,size(E,2));
        P3D(:,:,1:size(E,2)) = P(:,t).^2;
        P3D_full = repmat(P3D,size(X,1),size(X,1));
        term1(t) = -.5*logdet(sum(P3D_full.*cov,3));
        term2(t) = -.5*(X(:,t)-E*P(:,t))'*(sum(P3D_full.*cov,3)\(X(:,t)-E*P(:,t)));
    end

end

LogLikelihoodAll = term1 + term2;
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [DialedNum] = FrameSequeneRules(A);

%Computes the most probable dialed number based on design rules
%Input Variables
%           A    = Input Frame Decode Matrix
%
% Output Variables
%    DialedNum   = Number determined after analyzing a block of frames

% Checking for minimum number of arguments
if nargin < 1
    error('Not enough input arguments');
end

PossibleNums=cell(1,length(A));
RepeatSequence=[];
t=1;
k=1;

% Check to see what is repeating most number of times other than a pause
% signal
for i=1:(length(A)-1)
    if  strcmp(A(i),A(i+1)) && A(i) ~= 'p'
        k = k+1;
    else
        t=t+1;
        k = 1;
    end

       PossibleNums{t}=A(i);
       RepeatSequence(t)=k;
end

[R,C]=max(RepeatSequence);
DialedNum = PossibleNums{C};


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REAL-TIME BCILAB-LSL DEMONSTRATION SCRIPT
%
% This script demonstrates how to set up a default BCILAB-LSL pipeline for
% real-time data extraction and processing from an EEG device or file
%
% This script requires the following tools to be installed:
% BCILAB: https://www.sccn.ucsd.edu/wiki/BCILAB
% LSL:    https://code.google.com/p/labstreaminglayer/
%
% Optional (if using SourceLocalization)
% MobiLab: http://sccn.ucsd.edu/wiki/MoBILAB
%
% Author: Tim Mullen,  October 14, 2013, SCCN/INC/UCSD
% Email:  tmullen@ucsd.edu
%
% .........................................................................
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
% .........................................................................
%
% If you use this script as a basis for your published work, kindly
% consider acknowledging us in your manuscript. A possible example:
% "Thanks to Tim Mullen and Christian Kothe (SCCN, UCSD) for providing
%  scripting support"
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Force EEGLAB to use double precision
pop_editoptions( 'option_single', false);

%% Define config options (USER MUST DEFINE THIS)
% -------------------------------------------------------------------------
% If runlsl = true, then stream 'online' from device
% If runlsl = false then stream data from opts.PlaybackDataFile (below)
opts.runlsl = false;

% Run SIFT on the pipeline (following BCILAB filter pipeline)
opts.runSIFT = true;

% Set up the name of the LSL stream
opts.lsl.StreamName          = 'EEGDATA';
opts.lsl.SelectionProperty   = 'type';  % can also be 'name'
opts.lsl.SelectionValue      = 'EEG';

% Determine a fixed 'sliding' window length (seconds) for pulling data from LSL
% If winlen=0, pull the most recent (variable length) chunk of unprocessed data
opts.winlen = 1;

% suppress console output in processing loop?
% if false, this can reusult in a LOT of text printed to console
opts.silence = true;

% Establish file paths
% ..........................................................................
% User Tip:
%
% All paths are relative to opts.datapath. This is platform-independent and
% itself can be relative to bcilab root folder (i.e. data:/ is the userdata
% folder in the bcilab root dir)
% ..........................................................................
opts.datapath         = 'data:/';
% (required) calibration data set2
opts.TrainingDataFile = 'Simulation_EEGLAB_XVII.set'; %'MyCalibrationFile.set';
% (optional) file to "playback" if opts.runlsl = false
opts.PlaybackDataFile = 'Simulation_EEGLAB_XVII.set'; %'MyPlaybackFile.set';
% (optional) pipeline config file
opts.BCILAB_PipelineConfigFile = ''; 'MyFilteringPipeline.mat';
% (option) SIFT pipeline config file
opts.SIFT_PipelineConfigFile = ''; 'MySiftFilteringPipeline.mat';
% (optional) name of file containing MobiLab head model
opts.HeadModelDataFile= ''; %'NameOfHeadModelFile.mat';   % Look in <SIFT_RootFolder>\resources\headmodels\standard-Colin27-385ch.mat

%% Load head model object
% -------------------------------------------------------------------------
if ~isempty(opts.HeadModelDataFile)
    hmObj = headModel.loadFromFile(env_translatepath([opts.datapath opts.HeadModelDataFile]));
end

%% load calibration recording
% -------------------------------------------------------------------------
calibData = exp_eval_optimized(io_loadset([opts.datapath opts.TrainingDataFile], ...
                         'markerchannel',{'remove_eventchns',false}));
flt_pipeline('update');

%% Start LSL streaming
% -------------------------------------------------------------------------
if opts.runlsl == true
    run_readlsl('MatlabStream',opts.lsl.StreamName,                 ...
                'SelectionProperty',opts.lsl.SelectionProperty,     ...
                'SelectionValue',opts.lsl.SelectionValue);
else
    % ... OR read from datafile
    run_readdataset('MatlabStream',opts.lsl.StreamName,                         ...
                    'Dataset',io_loadset([opts.datapath opts.PlaybackDataFile], ...
                    'markerchannel',{'remove_eventchns',false}));
end

%% Set up the pipeline
% -------------------------------------------------------------------------
% load existing config file(s) (if it exists)
try    fltPipCfg = exp_eval(io_load([opts.datapath opts.BCILAB_PipelineConfigFile]));
catch, disp('-- no existing pipeline --'); fltPipCfg = {}; end
try    siftPipCfg = exp_eval(io_load([opts.datapath opts.SIFT_PipelineConfigFile]));
catch, disp('-- no existing pipeline --'); siftPipCfg = {}; end

% enforce usage of user-defined head model
if ~isempty(opts.HeadModelDataFile)
    fltPipCfg.psourceLocalize.hmObj = [opts.datapath opts.HeadModelDataFile]; end

% display GUI to allow user to configure the BCILAB pipeline...
fltPipCfg = arg_guipanel('Function',@flt_pipeline, ...
                         'Parameters',[{'signal',calibData} fltPipCfg], ...
                         'PanelOnly',false);

% display GUI to allow user to configure the SIFT pipeline...
if opts.runSIFT
    if isfield(calibData,'CAT')
        calibData = rmfield(calibData,'CAT');
    end
    siftPipCfg = arg_guipanel('Function',@onl_siftpipeline, ...
                             'Parameters',[{'EEG',calibData} siftPipCfg], ...
                             'PanelOnly',false);
    % force the modeling window length to agree with chunk length
    siftPipCfg.modeling.winlen = opts.winlen;
end
% ..........................................................................
% User Tip:
%
% You might want to save the pipeline config (fltPipCfg) at this point for
% later reload (from 'opts.BCILAB_PipelineConfigFile')
% save(env_translatepath([opts.datapath opts.SIFT_PipelineConfigFile]),'-struct','siftPipCfg');
% save(env_translatepath([opts.datapath opts.BCILAB_PipelineConfigFile]),'-struct','fltPipCfg');
% ..........................................................................

% ...apply the pipeline to calibration data
disp('-- Calibrating pipeline on training data (please wait) --');
cleaned_data = exp_eval(flt_pipeline('signal',calibData,fltPipCfg));

% ..........................................................................
% User Tip:
%
% Instead of using the GUI, you can also easily define pipelines in the
% script itself. Here are some examples:
%
% % remove flatline channels:
% noflats      = flt_clean_flatlines(flt_seltypes(calibData,'EEG'));
% % apply a high-pass filter:
% highpassed   = flt_fir(noflats,[0.5 1],'highpass');
% % discard noisy channels:
% goodchannels = flt_clean_channels('signal',highpassed,'min_corr',0.35,'ignored_quantile',0.3);
%
% % you always need to make sure to apply the pipeline
% cleaned_data = exp_eval(goodchannels);
% ..........................................................................

% ..........................................................................
% User Tip:
%
% After this step, cleaned_data will contain the results of having applied
% the entire pipeline to your calibration data. Thus, the steps up to this
% point essentially define an offline data processing pipeline which you
% can use for your publications. It is also helpful at this stage to
% examine the contents of cleaned_data for deficiencies, which can help you
% fine-tune your pipeline for later online use.
% ..........................................................................

% initialize the pipeline for streaming data
pipeline     = onl_newpipeline(cleaned_data,opts.lsl.StreamName);

% render a panel for viewing streams
gui_vis_filtered;

%% Main loop
% -------------------------------------------------------------------------
disp('-- Running pipeline on playback data --');
chunk_len = round(opts.winlen*cleaned_data.srate);
ConnData = [];
while true
    % grab a chunk of data from the stream and send through the filter pipeline
    [eeg_chunk,pipeline] = onl_filtered(pipeline, chunk_len, opts.silence);

    % process data chunk with SIFT
    if opts.runSIFT
        [eeg_chunk cfg] = hlp_scope({'is_online',1},@onl_siftpipeline,'EEG',eeg_chunk,siftPipCfg);
    end

    % .....................................................................
    % User Tip:
    %
    % The structure 'eeg_chunk' now contains the results of applying the
    % pipeline to the chunk of raw data acquired from the LSL stream.
    % This is an augmented EEGLAB data structure.
    %
    % - You can add code here to apply your own additional processing steps
    % - You'll also want to add code to store the results you'd like keep
    %
    % Hint:
    %  eeg_chunk.data contains processed/cleaned data
    %
    %  (if using source localization)
    %  eeg_chunk.srcpot contains current density for ROIs
    %
    %  (if using source localization and keepFullCsd is enabled)
    %  eeg_chunk.srcpot_all contains current density for all mesh vertices
    %
    %  (if using ICA)
    %  eeg_chunk.icaact contains ica activations
    %
    % (if using SIFT)
    % eeg_chunk.CAT  contains a SIFT object/datastructure
    %  CAT
    %   .Conn       Connectivity object
    %       .winCenterTimes
    %       .freqs
    %       .dims   dimensions for connectivity matrix
    %               e.g. (chans_to, chans_from, freqs, time)
    %       .<measure_name>    Connectivity matrix (conforms to dims above)
    %   .MODEL      Model structure
    %
    % .....................................................................

    % (optionally) give Matlab a little break...
    pause(0.0001);
end

repo: REST


function [net,net_asdn,poss,hardnegs] = mdnet_finetune_hnm_update(net,net_asdn,pos_data,neg_data,varargin)
%
%
global gpu;

opts.useGpu = gpu;
opts.conserveMemory = true ;
opts.sync = true ;

opts.maxiter = 30;
opts.learningRate = 0.001;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;

opts.batchSize_hnm = 256;
opts.batchAcc_hnm = 4;

opts.batchSize = 128;
opts.batch_pos = 32;
opts.batch_neg = 96;

opts = vl_argparse(opts, varargin) ;
% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

for i=1:numel(net.layers)
    if strcmp(net.layers{i}.type,'conv')
        net.layers{i}.filtersMomentum = zeros(size(net.layers{i}.filters), ...
            class(net.layers{i}.filters)) ;
        net.layers{i}.biasesMomentum = zeros(size(net.layers{i}.biases), ...
            class(net.layers{i}.biases)) ; %#ok<*ZEROLIKE>

        if opts.useGpu
            net.layers{i}.filtersMomentum = gpuArray(net.layers{i}.filtersMomentum);
            net.layers{i}.biasesMomentum = gpuArray(net.layers{i}.biasesMomentum);
        end
    end
end

for i=1:numel(net_asdn.layers)
    if strcmp(net_asdn.layers{i}.type,'conv')
        net_asdn.layers{i}.filtersMomentum = zeros(size(net_asdn.layers{i}.filters), ...
            class(net_asdn.layers{i}.filters)) ;
        net_asdn.layers{i}.biasesMomentum = zeros(size(net_asdn.layers{i}.biases), ...
            class(net_asdn.layers{i}.biases)) ; %#ok<*ZEROLIKE>

        if opts.useGpu
            net_asdn.layers{i}.filtersMomentum = gpuArray(net_asdn.layers{i}.filtersMomentum);
            net_asdn.layers{i}.biasesMomentum = gpuArray(net_asdn.layers{i}.biasesMomentum);
        end
    end
end

%% initilizing
if opts.useGpu
    one = gpuArray(single(1)) ;
else
    one = single(1) ;
end
res = [] ;
res_asdn=[];
res_test=[];

n_pos = size(pos_data,4);
n_neg = size(neg_data,4);
train_pos_cnt = 0;
train_neg_cnt = 0;

% extract positive batches
train_pos = [];
remain = opts.batch_pos*opts.maxiter;
while(remain>0)
    if(train_pos_cnt==0)
        train_pos_list = randperm(n_pos)';
    end
    train_pos = cat(1,train_pos,...
        train_pos_list(train_pos_cnt+1:min(end,train_pos_cnt+remain)));
    train_pos_cnt = min(length(train_pos_list),train_pos_cnt+remain);
    train_pos_cnt = mod(train_pos_cnt,length(train_pos_list));
    remain = opts.batch_pos*opts.maxiter-length(train_pos);
end

% extract negative batches
train_neg = [];
remain = opts.batchSize_hnm*opts.batchAcc_hnm*opts.maxiter;
while(remain>0)
    if(train_neg_cnt==0)
        train_neg_list = randperm(n_neg)';
    end
    train_neg = cat(1,train_neg,...
        train_neg_list(train_neg_cnt+1:min(end,train_neg_cnt+remain)));
    train_neg_cnt = min(length(train_neg_list),train_neg_cnt+remain);
    train_neg_cnt = mod(train_neg_cnt,length(train_neg_list));
    remain = opts.batchSize_hnm*opts.batchAcc_hnm*opts.maxiter-length(train_neg);
end

% learning rate
lr = opts.learningRate ;
lr_asdn = opts.learningRate*5;

% for saving positives
poss = [];

% for saving hard negatives
hardnegs = [];

% objective fuction
objective = zeros(1,opts.maxiter);

%% training on training set
% fprintf('\n');
for t=1:opts.maxiter
%     fprintf('\ttraining batch %3d of %3d ... ', t, opts.maxiter) ;
    iter_time = tic ;

    % ----------------------------------------------------------------------
    % hard negative mining
    % ----------------------------------------------------------------------
    score_hneg = zeros(opts.batchSize_hnm*opts.batchAcc_hnm,1);
    hneg_start = opts.batchSize_hnm*opts.batchAcc_hnm*(t-1);
    for h=1:opts.batchAcc_hnm
        batch = neg_data(:,:,:,...
            train_neg(hneg_start+(h-1)*opts.batchSize_hnm+1:hneg_start+h*opts.batchSize_hnm));
        if opts.useGpu
            batch = gpuArray(batch) ;
        end

        % backprop
        net.layers{end}.class = ones(opts.batchSize_hnm,1,'single') ;
        res = vl_simplenn(net, batch, [], res, ...
            'disableDropout', true, ...
            'conserveMemory', opts.conserveMemory, ...
            'sync', opts.sync) ;

        score_hneg((h-1)*opts.batchSize_hnm+1:h*opts.batchSize_hnm) = ...
            squeeze(gather(res(end-1).x(1,1,2,:)));
    end
    [~,ord] = sort(score_hneg,'descend');
    hnegs = train_neg(hneg_start+ord(1:opts.batch_neg));
    im_hneg = neg_data(:,:,:,hnegs);
%     fprintf('hnm: %d/%d, ', opts.batch_neg, opts.batchSize_hnm*opts.batchAcc_hnm) ;
    hardnegs = [hardnegs; hnegs];

%----------------------------yb added---------------------------
    batch_asdn=pos_data(:,:,:,train_pos((t-1)*opts.batch_pos+1:t*opts.batch_pos));
    if opts.useGpu
        batch_asdn = gpuArray(batch_asdn) ;
    end

    res_asdn=vl_simplenn(net_asdn,batch_asdn,[],res_asdn,...
            'disableDropout', true, ...
            'conserveMemory', opts.conserveMemory, ...
            'sync', opts.sync);
    feat_asdn=squeeze(gather(res_asdn(end-1).x(:,:,1,:)));

    num=size(feat_asdn,3);
    mask_asdn=ones(3,3,512,num,'single');
    for i=1:num
        feat_=feat_asdn(:,:,i);
        featlist=reshape(feat_,9,1);
        [~,idlist]=sort(featlist,'ascend');
        idxlist=idlist(1:3);

        for k=1:length(idxlist)
            idx=idxlist(k);
            row=floor((idx-1)/3)+1;
            col=mod((idx-1),3)+1;
            mask_asdn(col,row,:,i)=0;
        end
    end
    batch_asdn=batch_asdn.*mask_asdn;
%---------------------------------------------------------------

    % ----------------------------------------------------------------------
    % get next image batch and labels
    % ----------------------------------------------------------------------
    poss = [poss; train_pos((t-1)*opts.batch_pos+1:t*opts.batch_pos)];

    %batch = cat(4,pos_data(:,:,:,train_pos((t-1)*opts.batch_pos+1:t*opts.batch_pos)),...
        %im_hneg);
    batch = cat(4,batch_asdn,im_hneg);
    labels = [2*ones(opts.batch_pos,1,'single');ones(opts.batch_neg,1,'single')];
    if opts.useGpu
        batch = gpuArray(batch) ;
    end

    % backprop
    net.layers{end}.class = labels ;
    res = vl_simplenn(net, batch, one, res, ...
        'conserveMemory', opts.conserveMemory, ...
        'sync', opts.sync) ;

    % gradient step
    for l=1:numel(net.layers)
        if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

        net.layers{l}.filtersMomentum = ...
            opts.momentum * net.layers{l}.filtersMomentum ...
            - (lr * net.layers{l}.filtersLearningRate) * ...
            (opts.weightDecay * net.layers{l}.filtersWeightDecay) * net.layers{l}.filters ...
            - (lr * net.layers{l}.filtersLearningRate) / opts.batchSize * res(l).dzdw{1} ;

        net.layers{l}.biasesMomentum = ...
            opts.momentum * net.layers{l}.biasesMomentum ...
            - (lr * net.layers{l}.biasesLearningRate) * ....
            (opts.weightDecay * net.layers{l}.biasesWeightDecay) * net.layers{l}.biases ...
            - (lr * net.layers{l}.biasesLearningRate) / opts.batchSize * res(l).dzdw{2} ;

        net.layers{l}.filters = net.layers{l}.filters + net.layers{l}.filtersMomentum ;
        net.layers{l}.biases = net.layers{l}.biases + net.layers{l}.biasesMomentum ;
    end

    % print information
    objective(t) = gather(res(end).x)/opts.batchSize ;
    iter_time = toc(iter_time);
    fprintf('objective %.3f, %.2f s\n', mean(objective(1:t)), iter_time) ;

%%-------------------yb added----------------------
    iter_time = tic;
    net_fc=net;
    net_fc.layers=net_fc.layers(1:end-1);
    prob_k=zeros(9,1);
    for i=1:9
        row=floor((i-1)/3)+1;
        col=mod((i-1),3)+1;
        batch=pos_data(:,:,:,train_pos((t-1)*opts.batch_pos+1:t*opts.batch_pos));
        batch(col,row,:,:)=0;

        if opts.useGpu
            batch = gpuArray(batch) ;
        end

        res_test = vl_simplenn(net_fc, batch, [], res_test, ...
            'disableDropout', true, ...
            'conserveMemory', true, ...
            'sync', true) ;

        feat = gather(res_test(end).x) ;

        X=feat;
        E = exp(bsxfun(@minus, X, max(X,[],3))) ;
        L = sum(E,3) ;
        Y = bsxfun(@rdivide, E, L) ;
        prob_k(i)=sum(Y(1,1,1,:));
    end
    [~,idx]=min(prob_k);
    row=floor((idx-1)/3)+1;
    col=mod((idx-1),3)+1;

    batch=pos_data(:,:,:,train_pos((t-1)*opts.batch_pos+1:t*opts.batch_pos));
    labels = ones(3,3,1,opts.batch_pos,'single');
    labels(col,row,:)=0;

    if opts.useGpu
        batch = gpuArray(batch) ;
    end

    net_asdn.layers{end}.class = labels ;
    res_asdn = vl_simplenn(net_asdn, batch, one, res_asdn, ...
        'conserveMemory', opts.conserveMemory, ...
        'sync', opts.sync) ;

    for l=1:numel(net_asdn.layers)
        if ~strcmp(net_asdn.layers{l}.type, 'conv'), continue ; end

        net_asdn.layers{l}.filtersMomentum = ...
            opts.momentum * net_asdn.layers{l}.filtersMomentum ...
            - (lr_asdn * net_asdn.layers{l}.filtersLearningRate) * ...
            (opts.weightDecay * net_asdn.layers{l}.filtersWeightDecay) * net_asdn.layers{l}.filters ...
            - (lr_asdn * net_asdn.layers{l}.filtersLearningRate) / opts.batchSize * res_asdn(l).dzdw{1} ;

        net_asdn.layers{l}.biasesMomentum = ...
            opts.momentum * net_asdn.layers{l}.biasesMomentum ...
            - (lr_asdn * net_asdn.layers{l}.biasesLearningRate) * ...
            (opts.weightDecay * net_asdn.layers{l}.biasesWeightDecay) * net_asdn.layers{l}.biases ...
            - (lr_asdn * net_asdn.layers{l}.biasesLearningRate) / opts.batchSize * res_asdn(l).dzdw{2} ;

        net_asdn.layers{l}.filters = net_asdn.layers{l}.filters + net_asdn.layers{l}.filtersMomentum ;
        net_asdn.layers{l}.biases = net_asdn.layers{l}.biases + net_asdn.layers{l}.biasesMomentum ;
    end
    objective(t) = gather(res_asdn(end).x)/opts.batchSize ;
    iter_time = toc(iter_time);
    fprintf('asdn objective %.3f, %.2f s\n', mean(objective(1:t)), iter_time) ;

%------------------------------------------------------------
end % next batch


text:
function calpak_test014 ( )

%*****************************************************************************80
%
%% CALPAK_TEST014 tests JED_TO_YMDF_BAHAI and YMDF_TO_JED_BAHAI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST014\n' );
  fprintf ( 1, '  For the Bahai calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_BAHAI: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_BAHAI: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF	       JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_bahai ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_bahai ( jed1 );

      s2 = ymdf_to_s_numeric ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_bahai ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
function obj = get_metric(obj)

global Metric Nk Nt Nr Ns n;

V_equal = zeros(Nt,Ns);
W_equal = zeros(Nr,Ns,Nk);

if obj.V_RF == 0
    V_equal = obj.V_B;
    W_equal = obj.W_B;

else
    for k = 1:Nk
        V_equal = obj.V_RF*obj.V_B;
        W_equal(:,:,k) = obj.W_RF*obj.W_B(:,:,k);
    end
end

if (Metric.rate)
    [obj.rate(n),obj.A(:,n)] = get_scrate(V_equal,W_equal);
end

if (Metric.mse)
    obj.mse(n)= get_scmse(V_equal,W_equal);
end

if (Metric.ber)
    [obj.ber(n),obj.suberror(:,n)] = get_scber(V_equal,W_equal);
end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
% test harness for jacob_dot

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

if exist('p560') == 0
    mdl_puma560
end

e = 1e-4;
col = 1;
dq=[0 0 0 0 0 0];
dq(col) = 1;

Jd = (p560.jacob0(qn+dq*e)-p560.jacob0(qn))/e;
%Jd(:,col)
Jd
jacob_dot(p560, qn, dq)

author: petercorke

repo: robotics-toolbox-matlab


text:
function r8poly_values_horner_test ( )

%*****************************************************************************80
%
%% R8POLY_VALUES_HORNER_TEST tests R8POLY_VALUES_HORNER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8POLY_VALUES_HORNER_TEST\n' );
  fprintf ( 1, '  R8POLY_VALUES_HORNER evaluates a polynomial at a\n' );
  fprintf ( 1, '  point, using Horner''s method.\n' );

  m = 4;
  c = [ 24.0, -50.0, +35.0, -10.0, 1.0 ];
  r8poly_print ( m, c, '  The polynomial:' );

  n = 16;
  x_lo = 0.0;
  x_hi = 5.0;
  x = r8vec_linspace ( n, x_lo, x_hi );

  p = r8poly_values_horner ( m, c, n, x );

  r8vec2_print ( n, x, p, '  X, P(X):' );

  return
end


text:
function score = scoreModel_helper(metric, files, imnum, map, FixationDBCleaned, FIXATIONMAPS, STIMULI, HASH)

score = nan;

switch metric
        case 'ROC'
            % find the ROC performance
            [~, ~, fixations] = showAllEyeDataForImage2(files(imnum).name, FixationDBCleaned, STIMULI, HASH);
            score = AUC_Judd(map, fixations);

        case 'S'
            % find the similarity
            fixMap = im2double(imread(fullfile(FIXATIONMAPS, files(imnum).name)));
            score= similarity(map, fixMap);

        case 'EMD'
            % find the EMD performance
            fixMap = im2double(imread(fullfile(FIXATIONMAPS, files(imnum).name)));
            score = EMD(map, fixMap);

       case 'ROC_borji'
            % find Borji's ROC performance
            [~, ~, fixations] = showAllEyeDataForImage2(files(imnum).name, FixationDBCleaned, STIMULI, HASH);
            score = AUC_Borji(map, fixations);

        case 'sROC_borji'
            % find the shuffled ROC performance
            [~, ~, fixations] = showAllEyeDataForImage2(files(imnum).name, FixationDBCleaned, STIMULI, HASH);

            otherMap = zeros(size(fixations));

            % choose nOther other images to sample fixations from
            nOther = 10; % number of images to get fixations from
            otherImInd = [1:(imnum-1),(imnum+1):length(files)]; % indices of images other than the current one
            whichOther = otherImInd(randperm(length(otherImInd))); % randomize choice of images

            % accumulate fixations from nOther other images into one
            % fixation map: otherMap
            for ii = 1:nOther
                [~, ~, fixations2] = showAllEyeDataForImage2(files(whichOther(ii)).name, FixationDBCleaned, STIMULI, HASH);
                otherMap = otherMap + imresize(fixations2,size(otherMap));
            end

            % find the shuffled AUC score
            score = AUC_shuffled(map, fixations, otherMap);

    case 'CC'
            % find the cross correlation performance
            fixMap = im2double(imread(fullfile(FIXATIONMAPS, files(imnum).name)));
            score = CC(map,fixMap);

    case 'NSS'
            % find the normalized scanpath saliency
            [~, ~, fixations] = showAllEyeDataForImage2(files(imnum).name, FixationDBCleaned, STIMULI, HASH);
            score = NSS(map,fixations);
 end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function blockwrite(f)
%BLOCKWRITE  Append block to an existing file
%   Usage: blockwrite(f);
%
%   Input parameters:
%      f    : Block stream input.
%
%   Function appends *f* to a existing file. The file must have been
%   explicitly defined as the 'outfile' parameter of |block| prior
%   calling this function. If not, the function does nothing.
%
%   The function expect exactly the same format of *f* as is returned by
%   |blockread|.
%
%   See also: block

% Authors: Bjoern Ohl, Zdenek Prusa

complainif_notenoughargs(nargin,1,'BLOCKWRITE');

filestruct = block_interface('getOutFile');

if isempty(filestruct)
   % Do nothing if the file was not setup in block.
   return;
   %error('%s: Output file was not specified in block function.',...
   %      upper(mfilename));
end

filename = filestruct.filename;

% Reformat f if necessary
f = comp_sigreshape_pre(f,'BLOCKPLAY',0);

[L, W] = size(f);

Wread = filestruct.Nchan;

if Wread ~= W
    error(['%s: %s was initialized to work with %i channels but', ...
           ' only %i provided. '],upper(mfilename),filename,Wread,W);
end

if W>2
    error('%s: Cannot work with more than 2 channels.',upper(mfilename));
end

% prepare data depending on mono/stereo:
if W == 2           % stereo
  f = f.';
  f = f(:);
end

% flength = dlength + 36;
% We need to read one field from the header and
% update two.
fid = fopen(filename,'r+');

fseek(fid,40,-1);
dataLenInBytes = fread(fid,1,'uint32');
dataLenInBytes = dataLenInBytes + filestruct.alignment*L;
fileLenInBytes = dataLenInBytes + 36;

try
    if fseek(fid,4,-1) ~= 0
        error('d');
    end
    if fwrite(fid,fileLenInBytes,'uint32')<=0
        error('d');
    end

    if fseek(fid,40,-1)~=0
         error('d');
    end
    if fwrite(fid,dataLenInBytes,'uint32') <=0
         error('d');
    end
catch
    % We have to check whether the header was modified properly.
    error(['%s: An error has ocurred when modifying header of the ',...
          ' wav file. The file might be unreadable. Consider',...
          ' starting over.'],upper(mfilename));

end
fclose(fid);

% And now we can append the actual data
fid = fopen(filename,'a');
%write data into file (amplified by 2^15 to suit int16-range (from -2^15 to +2^15):

maxval = 1-(1/2^16);
minval = -1;

% Allow clipping since it is hard to do some sensible normalization.
f(f >= maxval)  = maxval;
f(f <= minval)  = minval;

% clipping check:
%if (max(tempvec) >= maxval) || (min(tempvec) <= minval)
    %We have no way how to find out how to properly normalize in blocks
    %warning('Clipping! Audio data limited to [-1, +1)');
%end

fwrite(fid, f*2^15, 'int16');
fclose(fid);    %close file


text:
i = imread('lab.pgm');

%Make image greyscale
if length(size(i)) == 3
	im =  double(i(:,:,2));
else
	im = double(i);
end

cs = fast_corner_detect_9(im, 30);
c = fast_nonmax(im, 30, cs);

image(im/4)
axis image
colormap(gray)
hold on
plot(cs(:,1), cs(:,2), 'r.')
plot(c(:,1), c(:,2), 'g.')
legend('9 point FAST corners', 'nonmax-suppressed corners')
title('9 point FAST corner detection on an image')


text:
function areaBB = computeArea(bb)
%computes area of the bb=[xmin ymin xmax ymax]

if ((bb(1) > bb(3)) || (bb(2) > bb(4)))
    areaBB = 0;
else
    areaBB = (bb(3) - bb(1) + 1) * (bb(4) - bb(2) + 1);
end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function parfor_sparseGV( qname, dbname, params )

[~, dbbasename, ~] = fileparts(dbname);
this_sparsegv_matname = fullfile(params.output.gv_sparse.dir, qname, [dbbasename, params.output.gv_sparse.matformat]);

if exist(this_sparsegv_matname, 'file') ~= 2
    %load features
    qfmatname = fullfile(params.input.feature.dir, params.data.q.dir, [qname, params.input.feature.q_sps_matformat]);
    if exist(qfmatname, 'file') ~= 2
        Iqname = fullfile(params.data.dir, params.data.q.dir, qname);
        [ f, d ] = Features_WUSTL( im2single(rgb2gray(imread(Iqname))) );
        d = relja_rootsift(d);
        [qfdir, ~, ~] = fileparts(qfmatname);
        if exist(qfdir, 'dir') ~= 7
            mkdir(qfdir);
        end
        save('-v6', qfmatname, 'f', 'd');
    end
    features_q = load(qfmatname);

    dbfmatname = fullfile(params.input.feature.dir, params.data.db.cutout.dir, [dbname, params.input.feature.db_sps_matformat]);
    if exist(dbfmatname, 'file') ~= 2
        Idbname = fullfile(params.data.dir, params.data.db.cutout.dir, dbname);
        [ f, d ] = Features_WUSTL( im2single(rgb2gray(imread(Idbname))) );
        d = relja_rootsift(d);
        [dbfdir, ~, ~] = fileparts(dbfmatname);
        if exist(dbfdir, 'dir') ~= 7
            mkdir(dbfdir);
        end
        save('-v6', dbfmatname, 'f', 'd');
    end
    features_db = load(dbfmatname);

    %geometric verification
    if size(features_db.d, 2) < 6
        H = nan(3, 3);
        inls_qidx = [];
        inls_dbidx = [];
        inliernum = 0;
        matches = [];
        inliers = [];
    else

        %geometric verification (homography lo-ransac)
        [matches, inliers, H, ~] = at_sparseransac(features_q.f,features_q.d,features_db.f,features_db.d,3,10);
        inliernum = length(inliers);
        inls_qidx = inliers(1, :); inls_dbidx = inliers(2, :);
    end

    %save
    if exist(fullfile(params.output.gv_sparse.dir, qname), 'dir') ~= 7
        mkdir(fullfile(params.output.gv_sparse.dir, qname));
    end
    save('-v6', this_sparsegv_matname, 'H', 'inliernum', 'inls_qidx', 'inls_dbidx', 'matches', 'inliers');

%     %debug
%     Iq = imread(fullfile(params.data.dir, params.data.q.dir, qname));
%     Idb = imread(fullfile(params.data.dir, params.data.db.cutout.dir, dbname));
%     figure();
%     ultimateSubplot ( 2, 1, 1, 1, 0.01, 0.05 );
%     imshow(rgb2gray(Iq));hold on;
%     plot(features_q.f(1, inls_qidx), features_q.f(2, inls_qidx),'g.');
%     ultimateSubplot ( 2, 1, 2, 1, 0.01, 0.05 );
%     imshow(rgb2gray(Idb));hold on;
%     plot(features_db.f(1, inls_dbidx), features_db.f(2, inls_dbidx),'g.');
%
%     keyboard;

end

end

author: HajimeTaira


text:
%{
B_LoopDeLoopCessna.m

This simply shows how the imported cessna patch object in A can be plugged directly into
the loop-de-loop example with no particular complications.

Matt Sheen, mws262@cornell.edu
%}

close all; clear all;
fig = figure;

%% Make the patch object:
patchData = stlread('cessna.stl');

p1 = patch(patchData,'EdgeColor','none','FaceColor','[0.8 0.8 1.0]','FaceAlpha',0.9);

axis equal
material('metal')
camlight('headlight');

p1.FaceAlpha = 0.5;

OrigVerts = p1.Vertices;

fig.Children.Clipping = 'off'; %Still see the plane when it leaves the axis area.
axis(2*[-30 30 -30 30 -10 50])
view(3)

%% Do the loop example.
vel = 5*[-1 0 0];
rot = angle2dcm(0,-0.04,0);
dt = 0.5;
center = [0 0 0];

totalRot = eye(3,3);

for i = 1:1000
    totalRot = totalRot*rot;
    center = dt*(totalRot*vel')' + center; %transform the velocity, keep track of the center position.
    p1.Vertices = (totalRot*OrigVerts')' + repmat(center,[size(OrigVerts,1),1]); %Rotate the patch object about its center and add new center point.

    pause(0.01);
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


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

%%
clear; close all; clc;

%% Syntax
% |[cMap]=magma(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function pref = determineDiscretization(N, lengthDom, pref)
%DETERMINEDISCRETIZATION    Determine discretization for a CHEBOP object.
%
%   PREFOUT = DETERMINEDISCRETIZATION(N, LENGTHDOM, PREFIN) choses the
%   correct discretization PREFOUT.DISCRETIZATION to be used when solving
%   problems (BVP/EIGS/EXPM), modeled by a CHEBOP N.
%
% See also CHEBOP/CLEARPERIODICBC, CHEBOP/SOLVEBVP, CHEBOP/EIGS.

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

% Case 1. A string was passed:
if ( isa(pref.discretization, 'char') )

    % Case 1.1. Determine the discretization if the user wants a discretization
    % using values:
    if ( strcmpi(pref.discretization, 'values') )
        % The default for the periodic case if TRIGCOLLOC. But, since TRIGCOLLOC
        % does not support breakpoints, it will be used only if there are no
        % breakpoints.
        if ( isa(N.bc, 'char') && strcmpi(N.bc, 'periodic') && lengthDom < 3 )
            pref.discretization = @trigcolloc;
        % Otherwise (i.e. periodic + breakpoints or other boundary
        % conditions), use CHEBCOLLOC2:
        else
            pref.discretization = @chebcolloc2;
        end

    % Case 1.2. Determine the discretization if the user wants a discretization
    % using coefficients:
    elseif ( strcmpi(pref.discretization, 'coeffs') )
        % Same here with TRIGSPEC and breakpoints:
        if ( isa(N.bc, 'char') && strcmpi(N.bc, 'periodic') && lengthDom < 3 )
            pref.discretization = @trigspec;
        else
            pref.discretization = @ultraS;
        end
    % Error:
    else
        error('CHEBFUN:CHEBOP:solvebvp:determineDiscretization', ...
        'PREF.DISCRETIZATION should be VALUES or COEFFS.\n')
    end

% Case 2. A function handle referring to a class was passed:
else
    % If TRIGCOLLOC or TRIGSPEC were chosen and there are breakpoints, we need
    % to throw an error:
    if ( ( isequal(pref.discretization, @trigcolloc) || ...
            isequal(pref.discretization, @trigspec) ) && lengthDom > 2 )
        error('CHEBFUN:CHEBOP:solvebvp:breakpointsInDomain', ...
        ['Problems with periodic boundary conditions where breakpoints \n', ...
        'are present cannot be solved with the TRIGCOLLOC/TRIGSPEC class.\n' ...
        'Please change the discretization to CHEBCOLLOC1/2 or ULTRAS.'])
    end
end

end