text:
function planC = insertDCMStruct(planC,dcmFileName)
%
%Insert structures from dcm file into the currently open planC.
%
% 06/05/2015
%
% 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/>.
%
% AI 2/13/20 Updated to account for pt orientation

global stateS

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

if ~exist('dcmFileName','var')
    [fileName, dirName] = uigetfile('*.dcm','Pick Structure file');
    dcmFileName = fullfile(dirName,fileName);
end

[dcmObj, isDcm] = scanfile_mldcm(dcmFileName);
dcmdirS = [];
if isDcm
    dcmdirS = dcmdir_add(dcmFileName, dcmObj, dcmdirS);
    dcmObj.clear;
else
    error('Invalid DICOM data')
end

% Read CERROptions.json
pathStr = getCERRPath;
optName = [pathStr 'CERROptions.json'];
optS = opts4Exe(optName);

scanOriS = planC{indexS.scan};
for scanNum = 1:length(scanOriS)
    imageOrientationPatient = scanOriS(scanNum).scanInfo(1).imageOrientationPatient;
    scanOriS(scanNum).imageOrientationPatient = imageOrientationPatient;
end

dataS = populate_planC_field('structures', dcmdirS.PATIENT, optS, scanOriS);

% Tolerance to determine oblique scan (think about passing it as a
% parameter in future)
numScans = length(planC{indexS.scan});
obliqTol = 1e-3;
isObliqScanV = zeros(1,numScans);

% Check scan to associate the strucutres
scanUIDc = {planC{indexS.scan}.scanUID};
scanTypesC = {};
for i = 1 : numScans
    scanTypesC{i} = [num2str(i) '.  ' planC{indexS.scan}(i).scanType];

    ImageOrientationPatientV = planC{indexS.scan}(i).scanInfo(1).imageOrientationPatient;

    % Check for obliqueness
    %     if ~isempty(ImageOrientationPatientV) && max(abs((abs(ImageOrientationPatientV) - [1 0 0 0 1 0]'))) <= obliqTol
    %         isObliqScanV(i) = 0;
    %     end
    if ~isempty(ImageOrientationPatientV)
        if max(abs((ImageOrientationPatientV(:) - [1 0 0 0 1 0]'))) < 1e-3
            isObliqScanV(i) = 0;
        elseif max(abs((ImageOrientationPatientV(:) - [-1 0 0 0 1 0]'))) < 1e-3
            isObliqScanV(i) = 0;
        elseif max(abs((ImageOrientationPatientV(:) - [-1 0 0 0 -1 0]'))) < 1e-3
            isObliqScanV(i) = 0;
        elseif max(abs((ImageOrientationPatientV(:) - [1 0 0 0 -1 0]'))) < 1e-3
            isObliqScanV(i) = 0;
        else
            isObliqScanV(i) = 1;
        end
    else
        if ~isempty(planC{indexS.scan}(i).scanInfo(1).imageOrientationPatient)
            ImageOrientationPatientV = planC{indexS.scan}(i).scanInfo(1).imageOrientationPatient;
        else
            ImageOrientationPatientV = [];
            isObliqScanV(i) = 1;
        end

    end

end

% Return if not a valid RTStruct file
if isempty(dataS)
    CERRStatusString('File not valid.')
    return
end

if ~ismember(dataS(1).assocScanUID,scanUIDc)
    if ~exist('dcmFileName','var')
        scanInd = listdlg('PromptString','Select Scan to associate structures',...
            'SelectionMode','single',...
            'ListString',scanTypesC);
    else
        scanInd = 1;
    end
    if ~isempty(scanInd)
        [dataS.assocScanUID] = deal(scanUIDc{scanInd});
    else
        return
    end
end
numStructs = length(planC{indexS.structures});

for i=1:length(dataS)
    dataS(i) = sortStructures(dataS(i),isObliqScanV,optS,planC);
    colorNum = numStructs + i;
    if isempty(dataS(i).structureColor)
        color = stateS.optS.colorOrder( mod(colorNum-1, size(stateS.optS.colorOrder,1))+1,:);
        dataS(i).structureColor = color;
    end
end
%Find any structures in dataS not already in planC.
toDeleteV = [];
for i=length(dataS):-1:1
    %     newStructName = temp_planC{temp_planC{end}.structures}(i).structureName;
    newStructData = dataS(i).contour;
    for j=1:length(planC{planC{end}.structures})
        oldStructName = planC{planC{end}.structures}(j).structureName;
        oldStructData = planC{planC{end}.structures}(j).contour;
        if isequal(newStructData, oldStructData)
            disp(['Structure "',oldStructName,'" already exists. Duplicate structure will not be inserted.']);
            toDeleteV = [toDeleteV, i];
            break;
        end
    end
end
dataS(toDeleteV) = [];
if ~isempty(planC{indexS.structures})
    for strNum = 1:length(dataS)
        planC{indexS.structures} = dissimilarInsert(planC{indexS.structures},dataS(strNum),length(planC{indexS.structures})+1);
    end
else
    planC{indexS.structures} = dataS;
end

editStructNumV = (numStructs+1):length(planC{indexS.structures});
planC = getRasterSegs(planC,editStructNumV);
%planC = setUniformizedData(planC);
for strNum = editStructNumV
    planC = updateStructureMatrices(planC, strNum);
end

if isfield(stateS, 'CERRFile')
    stateS.structsChanged = 1;
    CERRRefresh
end

CERRStatusString('Done inserting structures.')

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================

function [a,b] = fitRegressionLine(x, y)
    A11 = x' * x;
    A12 = sum(x);
    A21 = A12;
    A22 = size(x,1);

    B11 = y(:,1)' * x;
    B21 = sum(y(:,1));

    a = (A12*B21 - A22*B11) / (A12*A21 - A11*A22);
    b = (A11*B21 - A21*B11) / (A11*A22 - A12*A21);
end

author: hhyyti

repo: dcm-imu

save_path: github-repos/MATLAB/hhyyti-dcm-imu


text:
function truncated_normal_a_cdf_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_A_CDF_TEST tests TRUNCATED_NORMAL_A_CDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_A_CDF_TEST\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_A_CDF evaluates the lower Truncated Normal CDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean = mu\n' );
  fprintf ( 1, '    standard deviation = sigma\n' );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval [a,+oo)\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                                                  Stored         Computed\n' );
  fprintf ( 1, '       X        Mu         S         A             CDF             CDF\n' );
  fprintf ( 1, '\n');

  n_data = 0;

  while ( 1 )

    [ n_data, mu, sigma, a, x, cdf1 ] ...
      = truncated_normal_a_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    cdf2 = truncated_normal_a_cdf ( x, mu, sigma, a );

    fprintf( 1, '  %8.1f  %8.1f  %8.1f  %8.1f  %14g  %14g\n', x, mu, sigma, a, cdf1, cdf2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Stiffness Tensor
%
%% PLEASE HELP AND ADD CONTENT HERE
%

author: mtex-toolbox

repo: mtex

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


text:
function [m] = yd2m(yd)
% Convert length from yards to meters.
% Chad Greene 2012
m = yd*0.9144;

author: Sable

repo: mcbench-benchmarks


text:
function [dat] = read_edf(filename, hdr, begsample, endsample, chanindx)

% READ_EDF reads specified samples from an EDF datafile. It neglects all trial or
% data block boundaries as if the data was acquired in non-continuous mode.
%
% Note that since FieldTrip only accommodates a single sampling rate in a given
% dataset, whereas EDF allows specification of a sampling rate for each channel.  If
% there are heterogenous sampling rates then this function will automatically choose
% a subset.  If the last such channel is different from the rest, the assumption will
% be made that it is the annotation channel and the rest will be selected.  If that
% is not the case, then the largest subset of channels with a consistent sampling
% rate will be chosen.  To avoid this automatic selection process, the user may
% specify their own choice of channels using chanindx.  In this case, the automatic
% selection will only occur if the user selected channels still have heterogenous
% sampling rates.  In this case the automatic selection will occur amongst the user
% specified channels.  While reading the header the resulting channel selection
% decision will be stored in hdr.orig.chansel and the contents of this field will
% override chanindx during data reading.
%
% Use as
%   [hdr] = read_edf(filename)
% where
%    filename        name of the datafile, including the .edf extension
% This returns a header structure with the following elements
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel
%   hdr.orig         detailled EDF header information
%
% Or use as
%   [dat] = read_edf(filename, hdr, begsample, endsample, chanindx)
% where
%    filename        name of the datafile, including the .edf extension
%    hdr             header structure, see above
%    begsample       index of the first sample to read
%    endsample       index of the last sample to read
%    chanindx        index of channels to read (optional, default is all)
% This returns a Nchans X Nsamples data matrix
%
% Or use as
%   [evt] = read_edf(filename, hdr)
% where
%    filename        name of the datafile, including the .edf extension
%    hdr             header structure, see above
% This returns an Nsamples data vector of just the annotation channel

% Copyright (C) 2006, Robert Oostenveld and others
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

switch nargin
  case 1
    chanindx=[];
  case 2
    chanindx=[];
  case 3
    chanindx=begsample;
  case 4
end

needhdr = (nargin==1)||(nargin==3);
needevt = (nargin==2);
needdat = (nargin==5);

if needhdr
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read the header, this code is from EEGLAB's openbdf
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  FILENAME = filename;

  % defines Seperator for Subdirectories
  SLASH='/';
  BSLASH=char(92);

  cname=computer;
  if cname(1:2)=='PC' SLASH=BSLASH; end

  fid=fopen_or_error(FILENAME,'r','ieee-le');

  EDF.FILE.FID=fid;
  EDF.FILE.OPEN = 1;
  EDF.FileName = FILENAME;

  PPos=min([max(find(FILENAME=='.')) length(FILENAME)+1]);
  SPos=max([0 find((FILENAME=='/') | (FILENAME==BSLASH))]);
  EDF.FILE.Ext = FILENAME(PPos+1:length(FILENAME));
  EDF.FILE.Name = FILENAME(SPos+1:PPos-1);
  if SPos==0
    EDF.FILE.Path = pwd;
  else
    EDF.FILE.Path = FILENAME(1:SPos-1);
  end
  EDF.FileName = [EDF.FILE.Path SLASH EDF.FILE.Name '.' EDF.FILE.Ext];

  H1=char(fread(EDF.FILE.FID,256,'uint8=>char')');
  EDF.VERSION=H1(1:8);                          % 8 Byte  Versionsnummer
  %if 0 fprintf(2,'LOADEDF: WARNING  Version EDF Format %i',ver); end
  EDF.PID = deblank(H1(9:88));                  % 80 Byte local patient identification
  EDF.RID = deblank(H1(89:168));                % 80 Byte local recording identification
  %EDF.H.StartDate = H1(169:176);               % 8 Byte
  %EDF.H.StartTime = H1(177:184);               % 8 Byte
  EDF.T0=[str2num(H1(168+[7 8])) str2num(H1(168+[4 5])) str2num(H1(168+[1 2])) str2num(H1(168+[9 10])) str2num(H1(168+[12 13])) str2num(H1(168+[15 16])) ];

  % Y2K compatibility until year 2090
  if EDF.VERSION(1)=='0'
    if EDF.T0(1) < 91
      EDF.T0(1)=2000+EDF.T0(1);
    else
      EDF.T0(1)=1900+EDF.T0(1);
    end
  else
    % in a future version, this is hopefully not needed
  end

  EDF.HeadLen = str2num(H1(185:192));  % 8 Byte  Length of Header
  % reserved = H1(193:236);            % 44 Byte
  EDF.NRec = str2num(H1(237:244));     % 8 Byte  # of data records
  EDF.Dur  = str2num(H1(245:252));     % 8 Byte  # duration of data record in sec
  EDF.NS   = str2num(H1(253:256));     % 8 Byte  # of signals

  EDF.Label      = char(fread(EDF.FILE.FID,[16,EDF.NS],'uint8=>char')');
  EDF.Transducer = char(fread(EDF.FILE.FID,[80,EDF.NS],'uint8=>char')');
  EDF.PhysDim    = char(fread(EDF.FILE.FID,[ 8,EDF.NS],'uint8=>char')');

  EDF.PhysMin= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'uint8=>char')'));
  EDF.PhysMax= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'uint8=>char')'));
  EDF.DigMin = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'uint8=>char')'));
  EDF.DigMax = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'uint8=>char')'));

  % check validity of DigMin and DigMax
  if (length(EDF.DigMin) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Digital Minimum\n');
    EDF.DigMin = -(2^15)*ones(EDF.NS,1);
  end
  if (length(EDF.DigMax) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Digital Maximum\n');
    EDF.DigMax = (2^15-1)*ones(EDF.NS,1);
  end
  if (any(EDF.DigMin >= EDF.DigMax))
    fprintf(2,'Warning OPENEDF: Digital Minimum larger than Maximum\n');
  end
  % check validity of PhysMin and PhysMax
  if (length(EDF.PhysMin) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Physical Minimum, taking Digital Minimum instead\n');
    EDF.PhysMin = EDF.DigMin;
  end
  if (length(EDF.PhysMax) ~= EDF.NS)
    fprintf(2,'Warning OPENEDF: Failing Physical Maximum, taking Digital Maximum instead\n');
    EDF.PhysMax = EDF.DigMax;
  end
  idx_PhysMin_ge_PhysMax = EDF.PhysMin >= EDF.PhysMax;
  if (any(idx_PhysMin_ge_PhysMax))
    tmplabel = cellfun(@(x) [x ' '], cellstr(EDF.Label(idx_PhysMin_ge_PhysMax,:)),'UniformOutput',false)';
    fprintf(2,['Warning OPENEDF: Physical Minimum larger than Maximum.\nPLEASE recheck if the scaling and polarity in the following channels are still correct if used:\n' tmplabel{:} '\n']);
    %EDF.PhysMin = EDF.DigMin;
    %EDF.PhysMax = EDF.DigMax;
  end
  EDF.PreFilt= char(fread(EDF.FILE.FID,[80,EDF.NS],'uint8=>char')');
  EDF.SPR = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'uint8=>char')'));  % samples per data record

  fseek(EDF.FILE.FID,32*EDF.NS,0);

  EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./(EDF.DigMax-EDF.DigMin);
  EDF.Cal(isnan(EDF.Cal)) = 0;
  EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;

  %tmp = find(EDF.Cal < 0);
  %EDF.Cal(tmp) = ones(size(tmp));
  %EDF.Off(tmp) = zeros(size(tmp));

  EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
  %EDF.Calib=sparse(diag([1; EDF.Cal]));
  %EDF.Calib(1,2:EDF.NS+1)=EDF.Off';

  EDF.SampleRate = EDF.SPR / EDF.Dur;

  EDF.FILE.POS = ftell(EDF.FILE.FID);
  if EDF.NRec == -1                            % unknown record size, determine correct NRec
    fseek(EDF.FILE.FID, 0, 'eof');
    endpos = ftell(EDF.FILE.FID);
    EDF.NRec = floor((endpos - EDF.FILE.POS) / (sum(EDF.SPR) * 2));
    fseek(EDF.FILE.FID, EDF.FILE.POS, 'bof');
    H1(237:244)=sprintf('%-8i',EDF.NRec);      % write number of records
  end

  EDF.Chan_Select=(EDF.SPR==max(EDF.SPR));
  for k=1:EDF.NS
    if EDF.Chan_Select(k)
      EDF.ChanTyp(k)='N';
    else
      EDF.ChanTyp(k)=' ';
    end
    if contains(upper(EDF.Label(k,:)),'ECG')
      EDF.ChanTyp(k)='C';
    elseif contains(upper(EDF.Label(k,:)),'EKG')
      EDF.ChanTyp(k)='C';
    elseif contains(upper(EDF.Label(k,:)),'EEG')
      EDF.ChanTyp(k)='E';
    elseif contains(upper(EDF.Label(k,:)),'EOG')
      EDF.ChanTyp(k)='O';
    elseif contains(upper(EDF.Label(k,:)),'EMG')
      EDF.ChanTyp(k)='M';
    end
  end

  if isempty(chanindx)
    chanindx=1:EDF.NS;
  else
    if any(chanindx>EDF.NS)
      error('The selected channels are not present in the data');
    end
  end

  EDF.AS.spb = sum(EDF.SPR);    % Samples per Block
  bi=[0;cumsum(EDF.SPR)];

  idx=[];idx2=[];
  for k=1:EDF.NS
    idx2=[idx2, (k-1)*max(EDF.SPR)+(1:EDF.SPR(k))];
  end
  maxspr=max(EDF.SPR);
  idx3=zeros(EDF.NS*maxspr,1);
  for k=1:EDF.NS, idx3(maxspr*(k-1)+(1:maxspr))=bi(k)+ceil((1:maxspr)'/maxspr*EDF.SPR(k));end

  %EDF.AS.bi=bi;
  EDF.AS.IDX2=idx2;
  %EDF.AS.IDX3=idx3;

  % close the file
  fclose(EDF.FILE.FID);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % convert the header to Fieldtrip-style
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if all(EDF.SampleRate(chanindx)==EDF.SampleRate(chanindx(1)))
    chansel=chanindx;
    hdr.Fs           = EDF.SampleRate(chanindx(1));
    hdr.nChans       = length(chansel);
    hdr.label        = cellstr(EDF.Label(chansel,:));
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.NRec * EDF.SPR(chansel(1));
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = 1;
    hdr.chanunit     = cellstr(EDF.PhysDim(chansel,:));
    hdr.chantype     = repmat({'unknown'}, size(hdr.chanunit));  % start with unknown
    hdr.chantype(strcmp(hdr.chanunit, 'uV')) = {'eeg'};          % it might also be EOG, ECG, EMG, etc
    hdr.chantype(strcmp(hdr.chanunit, 'Boolean')) = {'trigger'};
    hdr.orig         = EDF;
    % this will be used on subsequent reading of data
    if length(chansel) ~= EDF.NS
      hdr.orig.chansel = chansel;
    else
      hdr.orig.chansel = 1:hdr.nChans;
    end
    hdr.orig.annotation = find(strcmp(cellstr(hdr.orig.Label), 'EDF Annotations'));

  elseif all(EDF.SampleRate(1:end-1)==EDF.SampleRate(1))
    % only the last channel has a deviant sampling frequency
    % this is the case for EGI recorded datasets that have been converted
    % to EDF+, in which case the annotation channel is the last
    chansel = find(EDF.SampleRate==EDF.SampleRate(1));
    % continue with the subset of channels that has a consistent sampling frequency
    hdr.Fs           = EDF.SampleRate(chansel(1));
    hdr.nChans       = length(chansel);
    ft_warning('Skipping "%s" as continuous data channel because of inconsistent sampling frequency (%g Hz)', deblank(EDF.Label(end,:)), EDF.SampleRate(end));
    hdr.label        = cellstr(EDF.Label(chansel,:));
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.NRec * EDF.SPR(chansel(1));
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = 1;
    hdr.orig         = EDF;
    % this will be used on subsequent reading of data
    hdr.orig.chansel    = chansel;
    hdr.orig.annotation = find(strcmp(cellstr(hdr.orig.Label), 'EDF Annotations'));

  else
    % select the sampling rate that results in the most channels
    [a, b, c] = unique(EDF.SampleRate);
    chancount = nan(size(a));
    for i=1:length(a)
      chancount(i) = sum(c==i);
    end
    [dum, indx] = max(chancount);
    chansel = find(EDF.SampleRate == a(indx));

    % continue with the subset of channels that has a consistent sampling frequency
    hdr.Fs           = EDF.SampleRate(chansel(1));
    hdr.nChans       = length(chansel);
    hdr.label        = cellstr(EDF.Label);
    hdr.label        = hdr.label(chansel);
    % it is continuous data, therefore append all records in one trial
    hdr.nSamples     = EDF.NRec * EDF.SPR(chansel(1));
    hdr.nSamplesPre  = 0;
    hdr.nTrials      = 1;
    hdr.orig         = EDF;
    % this will be used on subsequent reading of data
    hdr.orig.chansel    = chansel;
    hdr.orig.annotation = find(strcmp(cellstr(hdr.orig.Label), 'EDF Annotations'));

    ft_warning('channels with different sampling rate not supported, selecting subset of %d channels at %f Hz', length(hdr.label), hdr.Fs);
  end

  % return the header
  dat = hdr;

elseif needdat || needevt
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read the data
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % retrieve the original header
  EDF = hdr.orig;

  % There can be an optional chansel field containing a list of predefined channels.
  % These channels are in that case also the only ones represented in the FieldTrip
  % header, which means that the other channels are simply not visible to the naive
  % user. This field can be present because the user specified an explicit channel
  % selection in FT_READ_HEADER or because the read_edf function had to automatically
  % choose a subset to cope with heterogenous sampling rates or even both.  In any
  % case, at this point in the file reading process the contents of the chansel field
  % has the proper specification for channel selection, taking into account both the
  % user channel selection as well as any correction that might have been made due to
  % heterogenous sampling rates.

  if     ~isempty(chanindx) && ~isfield(EDF, 'chansel')
    % a subset of channels should been selected from the full list of channels in the file
    chanindx = chanindx; % keep as it is
    useChanindx = true;
  elseif ~isempty(chanindx) &&  isfield(EDF, 'chansel')
    % a subset of channels should been selected from the predefined list
    chanindx = EDF.chansel(chanindx);
    useChanindx = true;
  elseif  isempty(chanindx) &&  isfield(EDF, 'chansel')
    % all channels from the predefined list should be selected
    chanindx = EDF.chansel(chanindx);
    useChanindx = true;
  elseif  isempty(chanindx) && ~isfield(EDF, 'chansel')
    %  simply select all channels that are present in the file
    chanindx = 1:EDF.NS;
    useChanindx = false;
  end

  if needevt
    % read the annotation channel, not the data channels
    chanindx = EDF.annotation;
    begsample = 1;
    endsample = EDF.SPR(EDF.annotation)*EDF.NRec;
  end

  if useChanindx
    epochlength = EDF.SPR(chanindx(1));   % in samples for the selected channel
    blocksize   = sum(EDF.SPR);           % in samples for all channels
    chanoffset  = EDF.SPR;
    chanoffset  = round(cumsum([0; chanoffset(1:end-1)]));
    nchans      = length(chanindx);       % get the selection from the subset of channels
  else
    epochlength = EDF.SPR(1);             % in samples for a single channel
    blocksize   = sum(EDF.SPR);           % in samples for all channels
    nchans      = EDF.NS;                 % use all channels
  end

  % determine the trial containing the begin and end sample
  begepoch    = floor((begsample-1)/epochlength) + 1;
  endepoch    = floor((endsample-1)/epochlength) + 1;
  nepochs     = endepoch - begepoch + 1;

  % allocate memory to hold the data
  dat = zeros(length(chanindx),nepochs*epochlength);

  % read and concatenate all required data epochs
  for i=begepoch:endepoch
    if useChanindx
      % only a subset of channels with consistent sampling frequency is read
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      % read the complete data block
      buf = readLowLevel(filename, offset, blocksize); % see below in subfunction
      for j=1:length(chanindx)
        % cut out the part that corresponds with a single channel
        dat(j,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf((1:epochlength) + chanoffset(chanindx(j)));
      end

    elseif length(chanindx)==1
      % this is more efficient if only one channel has to be read, e.g. the status channel
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      offset = offset + (chanindx-1)*epochlength*2;
      % read the data for a single channel
      buf = readLowLevel(filename, offset, epochlength); % see below in subfunction
      dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf;

    else
      % read the data from all channels, subsequently select the desired channels
      offset = EDF.HeadLen + (i-1)*blocksize*2; % in bytes
      % read the complete data block
      buf = readLowLevel(filename, offset, blocksize); % see below in subfunction
      buf = reshape(buf, epochlength, nchans);
      dat(:,((i-begepoch)*epochlength+1):((i-begepoch+1)*epochlength)) = buf(:,chanindx)';
    end
  end

  % select the desired samples
  begsample = begsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
  endsample = endsample - (begepoch-1)*epochlength;  % correct for the number of bytes that were skipped
  dat = dat(:, begsample:endsample);

  % Calibrate the data
  if useChanindx
    calib = diag(EDF.Cal(chanindx));
  end
  if length(chanindx)>1
    % using a sparse matrix speeds up the multiplication
    dat = sparse(calib) * dat;
  else
    % in case of one channel the sparse multiplication would result in a sparse array
    dat = calib * dat;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for reading the 16 bit values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function buf = readLowLevel(filename, offset, numwords)
is_below_2GB = offset < 2*1024^3;
read_16bit_success = true;
if is_below_2GB
  % use the external mex file, only works for <2GB
  try
    buf = read_16bit(filename, offset, numwords);
  catch e
    read_16bit_success = false;
  end
end
if ~is_below_2GB || ~read_16bit_success
  % use plain matlab, thanks to Philip van der Broek
  fp = fopen(filename,'r','ieee-le');
  status = fseek(fp, offset, 'bof');
  if status
    ft_error(['failed seeking ' filename]);
  end
  [buf,num] = fread(fp,numwords,'bit16=>double');
  fclose(fp);
  if (num<numwords)
    ft_error(['failed reading ' filename]);
    return
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%kputanimate 'Non-Interactive Animation of Image Sequence'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros putanimate.pane file
%
% Parameters:
% InputFile: i 'Input File', required: 'File containing image sequence'
% InputFile: cmap 'Input Colormap ', optional: 'File containing alternate colormap filename'
% InputFile: clip 'Input Clip Gate', optional: 'File containing clip mask'
% Integer: x 'X Placement', default: -1: 'X location for GUI autoplacement'
% Integer: width 'Window Width ', default: 256: 'Window width'
% Integer: y 'Y Placement', default: -1: 'Y location for GUI autoplacement'
% Integer: height 'Window Height', default: 256: 'Window height'
% Double: update 'Update time', default: 2: 'How often to check input file for modification'
%
% Example: kputanimate({i, cmap, clip}, {'i','';'cmap','';'clip','';'x',-1;'width',256;'y',-1;'height',256;'update',2})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% putanimate - Non-Interactive Animation of Image Sequence
%
%  DESCRIPTION
% Putanimate is a visualization program that displays animated
% sequences of data.  It is expected that the input data has multiple
% frames organized in the Khoros Polymorphic Data Model along either
% depth, time, or elements.
%
% The input file containing the data to be animated is monitored by
% default; the file is checked every 2 seconds for change, and if it has
% been modified, the visualization display is updated accordingly.  The
% interval at which the input file is checked for modification may be
% specified using the [-update] argument.
%
% A clip mask may be used to dictate the portion of the data that is
% displayed in the animation using the [-clip] argument.
%
% An alternate input color map may be used to change the current color map
% using the [-cmap] argument.
%
% An animate direction value may be provided which specifies the direction
% in which images are sequenced using the [-dir] argument.  There are five
% different direction values which are as follows:
%
%  "" 5
% 1  -  "<<" Causes the animation to sequence in a backward direction.
% 2  -  "<"  Advances a single frame in a backwards direction (frame N to N-1).
% 3  -  "Stop"  Causes the animation to stop.
% 4  -  ">"  Advances a single frame in a forward direction (frame N to N+1)
% 5  -  ">>"  Causes the animation to sequence in a forward direction.
%
% An animate control value may be provided which controls how animation
% sequencing is performed when the animation is in motion using the [-control]
% argument.  Animation is considered "in motion" when the animation direction is
% sequencing forward or backward.  The animation may be performed once
% completely, repeated indefinitely, or performed in one direction then reversed
% and repeated in the other direction.  There are three different control values
% which are as follows:
%
%  "" 5
% 1  -  "Single" Does a single complete sequence of the animation in the
% current direction, and then stops.
% 2  -  "Loop"  Does a full sequence through the animation in the current
% animation direction.  As soon as the sequence is finished, it is started
% again, so that the animation is put into a loop.  This procedure will
% repeat until the animation is stopped.
% 3  -  "Autoreverse"  Does a single complete sequence of the animation, in the
% current animation direction, then reverses the direction and sequences back.
% This procedure will repeat until the animation is stopped.
%
% An animation speed value may be provided which controls how fast the frames
% in the image are to be sequenced using the [-speed] argument.  The speed is a
% value in seconds ranging between 0 and 5, where 0.0 is the fastest possible
% on a given machine, and 5.0 is very slow.  Note that 1/60 sec (0.016 sec)
% corresponds to the normal frame rate for television.  It is "NOT"
% recommended that the frame speed be set to 0.0 for most situations since this
% is an intensive operation that may make it difficult to re-gain control of
% the graphical user interface.
%
% By default, the visualization display will share its colormap will all
% other applications running at the same time; that is, it does not not use a
% .I private colormap,
% but rather makes use of the
% .I default colormap.
% It is possible to specify that the visualization display allocate its
% own private colormap or "grab" all available colors for its own use.
% This is done using the [-priv] argument.  When [-priv] is specified,
% moving the mouse pointer into the display window will cause the display
% to have its private colormap installed; moving the pointer out of the
% display window will cause the private colormap to be de-installed.  This
% results in the "technoflashing" phenomenon characteristic of private
% colormap installation.
%
% Color allocation can also be controlled and set to either
% .I  Read Only
% or
% .I Read Write
% using the [-alloc] argument.  When set to
% .I read only,
% once a color cell has been allocated, it can have its color set only once;
% from then on, the color cell can be shared by multiple applications, but
% not changed.  If the visualization display requires the color displayed
% to change, it must re-allocate the color cell, forcing a re-display of
% data.  This can be an expensive procedure.  In contrast, after a
% .I read/write
% color cell is allocated, it can have its color changed at any time
% without re-allocation; the data being displayed does not need to be
% redisplayed, and the color update process is much more efficient.
% However, the colors used in the visualization display cannot be shared by
% other applications.
%
% On creation, the animation display window may be placed manually (the
% default method), or placed automatically.  For automatic placement of the
% animation display window, specify the desired location in device (screen)
% coordinates, using the [-x] and [-y] arguments.
%
% While the animation display window should be created with a default size
% that is appropriate to display the data, a width and height for the window
% can be specified explicitly using the [-width] and [-height] arguments.
% Note that interactive resizing of the animate window using the window manager
% is currently NOT supported.
%
%
%
%  EXAMPLES
% % putanimate -i sequence:baby
% % putanimate -i sequence:bushes
%
%  "SEE ALSO"
% putdata(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [deltaTNew,xNew,tNew,k,dxdtCur,exitCode]=performOneAdaptiveRKStep(xCur,tCur,f,deltaT,deltaTMinMag,deltaTMaxMag,dxdtCur,order,solutionChoice,AbsTol,RelTol)
%%PERFORMONEADAPTIVERKSTEP Perform a single adaptive Runge-Kutta step,
%                   returning the new adjusted stepsize. This function is a
%                   subroutine of the function RKAdaptiveOverRange. Though
%                   most folks will just directly use the function
%                   RKAdaptiveOverRange to perform adaptive Runge-Kutta
%                   integration over a particular timespan, this subroutine
%                   has been broken out separately as one might want to
%                   program Runge-Kutta integration routines that
%                   adaptively integrate until a certain criterion is
%                   satisfied. As this function is meant as an efficient
%                   subroutine, none of the inputs provide default
%                   parameters if omitted.
%
%INPUTS: xCur The NX1 state vector at time tCur.
%        tCur The scalar current time of integration.
%           f The function handle for f(x,t)=dxdt over which integration is
%             to be performed. The output is NX1-dimensional.
%      deltaT The current stepsize to be taken in t. This can be positive
%             or negative.
% deltaTMinMag The minimum allowable magnitude of the step size.
% deltaTMaxMag The maximum allowable magnitude of the step size.
%      dxdtCur The value f(xCur,tCur). This is requested so that methods
%             that are FSAL (See comments to the function RungeKStep) can
%             avoid additional computations.
% order,solutionChoice  A pair of optional parameters that specify the
%              highest order of the embedded Runge-Kutta pair to use as
%              well as the specific algorithm to use. Details are given in
%              the comments to the RungeKStep function. If omitted or empty
%              matrices are passed, the default order of 5 is used and the
%              default solutionChoice of 0 is used.
%       RelTol The maximum relative error tolerance allowed, a positive
%              scalar (its use is explained in more detail below).
%       AbsTol The absolute error tolerance allowed, a positive scalar, of
%              the same for all components of x, or a positive NX1 vector.
%
%OUTPUTS: deltaTNew The value of deltaT that can be used for the next
%               adaptive step (i.e. pass it to this function). If the
%               function fails (exitCode~=0), then this will be set to the
%               last step size used.
%          xNew The updated NX1 state vector.
%          tNew The time of the updated state vector.
%             k The values of the derivatives f evaluated at various
%               points as determined by the algorithm used for the
%               step. This can be passed to functions like RKInterpPolys to
%               perform interpolation.
%      exitCode A code indicating whether the step could be successfully
%               performed. Possible values are
%               0: Integration was successful.
%               1: Unable to get a small enough step size.
%               3: Non-finite number encountered.
%
%The algorithm used for the step is described in the comments to the
%function RKAdaptiveOverRange.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %Determine the orders of the main and subsidiary embedded Runge-Kutta
    %formulae that were chosen. The order of convergence for the step size
    %testing is the smallest order taken. The isFSAL flag indicates whether
    %k(:,end) is the value of f evaluated at the next step.
    [orders,isFSAL]=RungeKStep(order,solutionChoice);
    RKOrder=min(orders);

    deltaTMag=abs(deltaT);
    deltaTSign=sign(deltaT);

    %The first time the choice in step size fails, it is adjusted the
    %"optimal" way in the Runge-Kutta-Fehlberg method. Additional times,
    %the step size is just halved in the hope that it will reach an
    %accepted value more quickly.
    failedReducingStepSize=false;
    moveOnToNextStep=false;
    while(moveOnToNextStep==false)
        %The returned value dxValdt in k(:,1) is for xCur and tCur, which
        %is curStep-1.
        [xPredMain,xPredSubsid,k]=RungeKStep(xCur,tCur,f,deltaT,dxdtCur,order,solutionChoice);

        %Integration can only be over finite functions.
        if(any(~isfinite(xPredMain))||any(~isfinite(xPredSubsid)))
            xNew=[];
            tNew=[];
            exitCode=3;
            deltaTNew=deltaTMag;
            return;
        end

        %The local error estimate. This must be transformed into a
        %combination relative/ absolute error term to determine whether
        %the step should be rejected.
        normFactor=max(max(abs(xPredMain),abs(xCur)),AbsTol/RelTol);
        theError=max(abs((xPredMain-xPredSubsid)./normFactor));

        if(theError>RelTol)
            if(deltaTMag<deltaTMinMag)
                %If the step size got too small, then return.
                xNew=[];
                tNew=[];
                exitCode=1;
                deltaTNew=[];
                return;
            end

            if(failedReducingStepSize==false)
                failedReducingStepSize=true;

                %The Fehlberg step reduction (using the relative error).
                deltaTMag=max(deltaTMinMag,deltaTMag*max(0.1, 0.8*(RelTol/theError)^(1/RKOrder)));
            else
                %Just halve the step size.
                deltaTMag=deltaTMag/2;
            end
            deltaT=deltaTSign*deltaTMag;
        else
            %The step was successful; leave the loop.
            moveOnToNextStep=true;
        end
    end

    %The step is successful, save the results from the step, including
    %information so that interpolation can be perfromed, if necessary.
    xNew=xPredMain;
    tNew=tCur+deltaT;

    %Save the current value to be reused on the next step, if the method
    %is an FSAL function, so that an evaluation of f can be avoided.
    if(isFSAL)
        dxdtCur=k(:,end);
    else
        dxdtCur=f(xNew,tNew);
    end

    %If a step is successful, then increase the step size for the
    %next step in the standard manner used with Runge-Kutta-
    %Fehlberg methods, but limit the maximum size of the increase
    %to a scale factor of 4. This avoid huge step sizes when the
    %predicted error is very small.
    deltaTMag=min(deltaTMaxMag,deltaTMag*min(4,0.8*(RelTol/theError)^(1/RKOrder)));

    %Since the minimum step size changes every loop, this makes sure that
    %deltaT does not go beneath it just because the loop changed.
    deltaTMag=max(deltaTMag,deltaTMinMag);
    deltaTNew=deltaTMag*deltaTSign;

    exitCode=0;
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:
% Author Deepak.
% Just provide the path of the image to add the noise or to restore the
% image.
clc;
f = imread('Place the path of your image file');
figure
imshow(f),title('Original Image')
[M N] = size(f);
% Any type of noise can b added to the image provided. Just see the
% imnoise2 file to see the various noise effects. Here an example is shown
% where we are adding the pepper noise to the image.
r = imnoise2('salt & pepper',M,N,0,0.1);
figure
imshow(r),title('Noise to be added');
c = find(r == 1);
gp = f;
gp(c) = 255;
figure
imshow(gp),title('Image after adding the Noise');
% The image distorted by adding the noise can be restored by using the
% function imrest. Here an example is shown by using contraharmonic filter.
% you can use various type of other filter to restore the image. For more
% detail just see the funcion imrest.
fp = imrest(gp,'chmean',3,3,-5.5);
figure
imshow(fp),title('Image after restoration.')


text:
function legendre_polynomial_test08 ( )

%*****************************************************************************80
%
%% LEGENDRE_POLYNOMIAL_TEST08 tests PMN_POLYNOMIAL_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2012
%
%  Author:
%
%    John Burkardt
%
  mm = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_POLYNOMIAL_TEST08:\n' );
  fprintf ( 1, '  PMN_POLYNOMIAL_VALUES stores values of\n' );
  fprintf ( 1, '  the Legendre polynomial Pmn(n,m,x).\n' );
  fprintf ( 1, '  PMN_POLYNOMIAL_VALUE evaluates the polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                                Tabulated                 Computed\n' );
  fprintf ( 1, '     N     M        X           Pmn(N,M,X)                Pmn(N,M,X)             Error\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, m, x, fx1 ] = pmn_polynomial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    v = pmn_polynomial_value ( mm, n, m, x );
    fx2 = v(1,n+1);

    e = fx1 - fx2;

    fprintf ( 1, '  %4d  %4d  %12g  %24g  %24g  %8g\n', n, m, x, fx1, fx2, e );

  end

  return
end


text:
function [GraphObj]=createMetIntrcNetwork(model,metAbbr,varargin)
% Create directed metabolite-metabolite interaction network using given metabolites and the model.
% The produced network consists of given metabolites and its first neighbours.
% Colour and width of edges will be adjusted based on flux values If fluxes of the model are given.
% Another metabolite-metabolite interaction network will be generated if the node is clicked;
%
% Left click :Generate sub-metabolite-metabolite interaction network from the created figure.
%            This functionality was added for better looking at the created network
%            and showing flux values on edges lines.
%
% Right click :Generate metabolite-metabolite interaction network from model. This
%             functionality was added for creating metabolite-metabolite network using
%             clicked metabolite and model.
% The produced figure by right and left click on the main figure also
% has same property with the main figure, and it is clickable too.
%
% Note : If the fluxes are not given, irreversible COBRA model structure should use for
%        representing reversible reactions. See convertToIrreversible function at github page of cobratoolbox
%        https://opencobra.github.io/cobratoolbox/stable/modules/reconstruction/refinement/
% USAGE:
%
%    [GraphObj]=createMetIntrcNetwork(model,metAbbr,varargin)
%
% INPUTS:
%
%    model:             COBRA model structure
%    metAbbr:           List of metabolite abbreviation as a cell array
%
% OPTIONAL INPUT:
%    varargin:      Optional Inputs provided as 'ParameterName', Value
%                   pairs. the following parameternames are available:
%
%                   * fluxes:             flux vector
%                   * Graphtitle:         Title for a figure as a string
%                   * excludedMets:       As a cell array, Metabolite abbreviations that desired for exclude
%                   * nodeColour:         Colour for nodes (RGB Triplet) e.g. [0 1 0]
%                                         see https://www.mathworks.com/help/matlab/ref/colorspec.html
%                   * Hnodecolour:         Colour for nodes that will be highlighted e.g. [0 1 1]
%                   * scaleMin            Minimum value for scaling colorbar
%                   * scaleMax            Maximum value for scaling colorbar
%                   * nodeSize            Size of Nodes
%                   * HnodeSize           Highlighted node size
%                   * arrowSize           Arrow Size
%                   * threshold           Treshold for edges, the edges that has the flux value
%                                         below the treshold will be deleted, and if the node don't
%                                         have any relationship with any other nodes will be deleted.
%                   * excNodesWithDeg     It deletes nodes that has not degree in a given range.
%                                         It should be given as [minimum degree maximum degree].
%                                         Example, [1 50]. it deletes nodes that has lower than 1 degree
%                                         and higher than 50 degree.
%
% OUTPUT:
%    GraphObj:              Matlab digraph object
%
%
% EXAMPLES:
%    1) create network with succ[c] and akg[c]
%       [GraphObj]=createMetIntrcNetwork(model,{'succ[c]','akg[c]'})
%    2) create network with flux values
%       [GraphObj]=createMetIntrcNetwork(model,{'succ[c]','akg[c]'},'fluxes',fluxvector)
%    3) create network with flux values, and exculude exclude very employed metabolites
%       [GraphObj]=createMetIntrcNetwork(model,{'succ[c]','akg[c]'},'fluxes',fluxvector,'excludedMets',{'atp[c]','adp[c]'})
%    4) create network with scaling colorbar, [0 1]
%       [GraphObj]=createMetIntrcNetwork(model,{'succ[c]','akg[c]'},'fluxes',fluxvector,'excludedMets',{'atp[c]','adp[c]'},'scaleMax',1,'scaleMin',0)
%
% .. Author: - Kadir KOCABAS November/2020

defaultGraphtitle ='Metabolite-Metabolite Interaction Network';
defaultexcludedMets ={};
defaultFluxes=ones(size(model.rxns,1),1);
defaultnodeColour=[0.9 0.7 1];
defaultHnodecolour=[0 1 0];
defaultscaleMax=1e-6;
defaultscaleMin=0;
defaultNodeSize=15;
defaultHNodeSize=25;
defaultarrowSize=9;
defaultThreshold=1e-7;
defaultexcNodesWithDeg=[0 1e+6];

p = inputParser;
addRequired(p,'model',@isstruct);
addRequired(p,'metAbbr',@iscell);
addParameter(p,'fluxes',defaultFluxes,@isvector);
addParameter(p,'Graphtitle',defaultGraphtitle,@(x) isstring(x) || ischar(x));
addParameter(p,'excludedMets',defaultexcludedMets,@iscell);
addParameter(p,'nodeColour',defaultnodeColour,@isvector);
addParameter(p,'Hnodecolour',defaultHnodecolour,@isvector);
addParameter(p,'nodeSize',defaultNodeSize,@isnumeric);
addParameter(p,'HnodeSize',defaultHNodeSize,@isnumeric);
addParameter(p,'arrowSize',defaultarrowSize,@isnumeric);
addParameter(p,'scaleMax',defaultscaleMax,@isnumeric);
addParameter(p,'scaleMin',defaultscaleMin,@isnumeric);
addParameter(p,'threshold',defaultThreshold,@isnumeric);
addParameter(p,'excNodesWithDeg',defaultexcNodesWithDeg,@isvector);

parse(p,model,metAbbr,varargin{:});
p.KeepUnmatched = true;
model=p.Results.model;
metAbbr=p.Results.metAbbr;
fluxes=p.Results.fluxes;
Graphtitle=p.Results.Graphtitle;
excludedMets=p.Results.excludedMets;
nodeColour=p.Results.nodeColour;
Hnodecolour=p.Results.Hnodecolour;
scaleMin=p.Results.scaleMin;
scaleMax=p.Results.scaleMax;
nodeSize=p.Results.nodeSize;
HnodeSize=p.Results.HnodeSize;
arrowSize=p.Results.arrowSize;
threshold=p.Results.threshold;
excNodesWithDeg=p.Results.excNodesWithDeg;

%Controlling scale values
if ~(scaleMax>scaleMin)
    error('scaleMax should bigger than scaleMin');
end

%controlling degre values
if (excNodesWithDeg(1)>excNodesWithDeg(2))
    error('maximum degree should bigger or equal to minimum degree');
end

%controlling if flux vakues given
if (sum(fluxes)==size(model.rxns,1)&all(ismember(fluxes,1)))
    fcont=1;

else
    fcont=0;
end

if ((fcont==1) & ~(scaleMax==1e-6 && scaleMin==0) )
    warning('Scales are given without giving fluxes')
end

if ~all(ismember(metAbbr,model.mets))
    str = strjoin( metAbbr(~ismember(metAbbr,model.mets)), ' ' );
    warning(['These metabolites are not found in model : ' str])
end

%if user wants to delete edges with 0 fluxes treeshold assign 1e-5. Because
%weight of edges with zero fluxes shown as 1e-6 since weight of edges cannot
%be 0 in the matlab graph objects

if threshold==0
    threshold=1e-5;
end

%Create table, creating table without fluxrates if fluxes are not supplied
IndexesofMets=find(ismember(model.mets,metAbbr));
metMatrix=~ismember(model.S(IndexesofMets,:),0);
Rxns=model.rxns(any(metMatrix,1));
FluxRes=fluxes(ismember(model.rxns,Rxns));
%
%Get sub S matrix from S matrix
RxnsMatrx=~ismember(model.S(:,any(metMatrix,1)),0);
InvolvedMets=model.mets(any(RxnsMatrx,2));
adjunc_mat=model.S(any(RxnsMatrx,2),any(metMatrix,1));
%
%Multiply sub S matrix with flux values to keep flux value information
FluxMatrix  =FluxRes+1e-6;
FluxMatrix = repmat(FluxMatrix',size(adjunc_mat,1),1);
adjunc_mat=adjunc_mat.*FluxMatrix;

%Create left and right matrix( left matrix : comsumed metabolites, righ matrix:
%produced metabolites)
leftMatrix=adjunc_mat;
leftMatrix(leftMatrix > 0)=0;
rightMatrix=adjunc_mat;
rightMatrix(rightMatrix < 0)=0;
rightMatrix(rightMatrix > 0)=1;

%create adjacency matrix
A=(leftMatrix*rightMatrix')*(-1);

%Delete metabolites that is not first neighbour of metAbbr
idx=ismember(InvolvedMets,metAbbr);
newadjMtrx=zeros(size(A));
newadjMtrx(idx,:)=A(idx,:);
newadjMtrx(:,idx)=A(:,idx);
rows=any(newadjMtrx,1);
cols=any(newadjMtrx,2);
indofMets=rows'|cols;
metnames=InvolvedMets(indofMets);
newadjMtrx=newadjMtrx(indofMets,indofMets);
A=newadjMtrx;

%create graph object.
G = digraph(A,metnames);

GraphObj=G;
excludedMets=excludedMets(~ismember(excludedMets,metAbbr));

%Exclude given metabolites
G=rmnode(G,excludedMets);

if threshold~=1e-7
    edesBelowTresholdIdx=find(G.Edges.Weight<threshold);
    G=rmedge(G,edesBelowTresholdIdx);
end

%controlling if there is a node with given range
nodesIndegree=indegree(G);
nodesoutdegree=outdegree(G);
totalDegree=nodesIndegree+nodesoutdegree;
minDegree=excNodesWithDeg(1);
maxDegree=excNodesWithDeg(2);

if (max(totalDegree) < minDegree | min(totalDegree) > maxDegree)
    error('There is no node with degree in the given range')
end

G=rmnode(G,find(totalDegree<minDegree | totalDegree>maxDegree ));

figure;
hold on

%chech if there is a edge in the network
if ~isempty(G.Edges.Weight)
    G.Edges.LWidths = (G.Edges.Weight-min(G.Edges.Weight))/(max(G.Edges.Weight)-min(G.Edges.Weight))+0.00001;
    edgeCont=true;
else
    warning('There is no edges in the Network');
    edgeCont=false;
end

%create plot
h=plot(G,'MarkerSize',nodeSize,'NodeColor',nodeColour,'ArrowSize',arrowSize,'NodeLabelMode','auto','Layout','force','UseGravity',true);
metAbbr=metAbbr(ismember(metAbbr,G.Nodes.Name));
highlight(h,metAbbr,'NodeColor',Hnodecolour,'MarkerSize',HnodeSize);

nl = h.NodeLabel;
h.NodeLabel = '';
xd = get(h, 'XData');
yd = get(h, 'YData');
title(Graphtitle,'Interpreter', 'none');
txt=text(xd, yd, nl, 'FontSize',8, 'FontWeight','bold', 'HorizontalAlignment','center', 'VerticalAlignment','middle');
set(txt,'Interpreter', 'none');
%Define function for texts
set(txt,'ButtonDownFcn',{@hit_node,A,metnames,model,fluxes,excludedMets,Graphtitle,nodeColour,Hnodecolour,fcont,scaleMin,scaleMax,nodeSize,HnodeSize,arrowSize,threshold,excNodesWithDeg});

%define color and width of edges if fluxes are given
if edgeCont % control if there is edge in the network
    if  fcont==1
        h.LineWidth=1;
    elseif ~any(isnan(G.Edges.LWidths))
        h.LineWidth=G.Edges.LWidths*2.5;
        colormap jet(10)
        h.EdgeCData=G.Edges.Weight;
        hcb=colorbar;
        colorTitleHandle = get(hcb,'Title');
    else
        h.LineWidth=1;
        colormap jet(10)
        h.EdgeCData=G.Edges.Weight;
        hcb=colorbar;
        colorTitleHandle = get(hcb,'Title');

    end

    if fcont==0
        if ~(scaleMax==1e-6 && scaleMin==0)
            caxis([scaleMin scaleMax])
            titleBar = ['Scaled Fluxes',' ', '(', num2str(scaleMin),' - ',num2str(scaleMax),')'];
            set(colorTitleHandle ,'String',titleBar,'FontWeight','Bold');
        else
            set(colorTitleHandle ,'String','Fluxes','FontWeight','Bold');
        end
    end
end

set(gca,'XTickLabel',{' '});
set(gca,'YTickLabel',{' '});
set(gca,'YTick',[]);
set(gca,'XTick',[]);
set(gca,'XColor', 'none','YColor','none');
end

author: opencobra

repo: cobratoolbox


text:
filename='Cantilever_tetrahedra';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.15;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef LDClassifier < Algorithm

    properties (Access = public)
        shouldTrain = false;
    end

    properties (Access = private)
        classifier;
    end

    methods (Access = public)

        function obj = LDClassifier()
            obj.name = 'LinearDiscriminant';
            obj.inputPort = DataType.kTable;
            obj.outputPort = DataType.kTable;
        end

        function dataOut = compute(obj,data)
            if obj.shouldTrain
                obj.train(data);
                dataOut = [];
            else
                dataOut = obj.test(data);
            end
        end

        function train(obj,table)

            obj.classifier = fitcdiscr(...
                table.features, ...
                table.label, ...
                'DiscrimType', 'linear', ...
                'Gamma', 0, ...
                'FillCoeffs', 'off');
        end

        function metrics = computeMetrics(obj,table)
            flops = timeit(@()obj.test(table)) / Constants.kReferenceComputingTime;
            memory = Helper.ComputeObjectSize(obj.classifier);
            outputSize = table.height * Constants.kClassificationResultBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function labels = test(obj,table)
            labels = predict(obj.classifier,table.features);
        end

        function str = toString(obj)
            str = sprintf('%s',obj.name);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%wrap - Set radian angles in range [0,2pi] or [-pi,pi].
%
%  USAGE
%
%    y = wrap(x,range)
%
%    x              angles in radians
%    range          optional:  1 for [-pi,pi] (default)
%                              2 for [0,2pi]
%
%  SEE ALSO
%
%    See also isradians.
%

% Copyright (C) 2010-2011 by Michaël Zugaro
%
% 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.

function y = wrap(x,range)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help wrap'' for details).');
end

if nargin < 2,
	range = 1;
end

if ~isa(x,'double'), y = []; return; end

% Determine angle in [0,2*pi]
y = mod(x,2*pi);

% Change range if necessary
if range == 1,
	change = y > pi;
	y(change) = y(change)-2*pi;
end


text:
function C = cs_droptol (A, tol)                                            %#ok
%CS_DROPTOL remove small entries from a sparse matrix.
%   C = cs_droptol(A,tol) removes entries from A of magnitude less than or
%   equal to tol.  Same as A = A .* (abs (A) >= tol).
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ;
%       cspy (abs (A) >= 1e-10) ;
%       C = cs_droptol (A, 1e-10) ;
%       cspy (C) ;
%
%   See also: RELOP, ABS

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_droptol mexFunction not found') ;


text:
function [sshcmd,bsubcmd,info] = wrapjobcmd_bsub_ssh(jobcmd,varargin)

[info.jobname,info.resfile,info.reserrfile,info.ncores] = myparse(varargin,...
  'jobname',sprintf('job%s',datestr(now,'yyyymmddTHHMMSS')),...
  'resfile',[tempname,'.out'],...
  'reserrfile',[tempname,'.err'],...
  'ncores',1);

bsubcmd = sprintf('bsub -n %d -J %s -o %s -e %s "%s"',info.ncores,info.jobname,info.resfile,info.reserrfile,strrep(jobcmd,'"','\"'));
sshcmd = sprintf('ssh login1 "%s"',strrep(strrep(bsubcmd,'\','\\'),'"','\"'));

author: kristinbranson

repo: JAABA


text:
%{
helix_example.m
Example of how to use rebin_helix.m for helical cone-beam CT reconstruction.
Copyright 2010-7-21, Gregory Handy and Jeff Fessler, University of Michigan
%}

if 0 % test a particular geometry (case 2795)
	f.down = 4; % down sample a lot to save time
	cg = ct_geom('ge2', 'nt', 32, 'na', 3625, ...
...%		'dfs', inf, ...
		'dfs', 0, ...
		'down', f.down, ...
		'pitch', .53125, 'source_z0', -50.5494160, ...
		'orbit', 3625/984*360, ... % 1326.21951219512195,
		'orbit_start', 109.12139);

	ig = image_geom('nx', 320, 'ny', 320, 'nz', 3*64, ... % 61
		'down', f.down, ...
		'dx', 2.191162, 'dz', 0.625, 'offset_z', 56.0);
%	cg.plot(ig), return
end

if ~isvar('cg'), printm 'cg: cone-beam CT geometry'
	f.down = 4; % down sample a lot to save time
	f.nturn = 12;
	cg = ct_geom('ge2', 'pitch', 0.5, 'source_z0', -100, ...
		'na', 984*f.nturn, 'orbit', 360*f.nturn, 'orbit_start', 17, ...
		'down', f.down);
end

if ~isvar('ig'), printm 'ig: image geometry'
	ig = image_geom('nx', 256, 'ny', 256, 'nz', 160, ...
			'dx', 2, 'dz', 0.625, 'down', f.down);
	im clf, cg.plot(ig);
prompt
end

if ~isvar('ell'), printm 'ell: ellipsoid object'
	z0 = ig.offset_z * ig.dz;
	ell = [ ...
		[0 0 -z0	[0.3 0.1]*ig.fov 0.6*ig.zfov	0 0 1000];
		[80 10 -z0	30 30 10	0 0 1000];
		[-10 -40 75	40 40 40	0 0 1000];
		[-10 80 -20	30 30 30	0 0 1000];
		[0 0 -36	20 20 5	0 0 1000];
		[0 0 -12	20 20 5	0 0 1000];
		[0 0 12	20 20 5	0 0 1000];
		[0 0 36	20 20 5	0 0 1000];
	];
	clear z0
end

if ~isvar('xtrue'), printm 'xtrue: true image volume'
	xtrue = ellipsoid_im(ig, ell, 'oversample', 2);

	clim = [0 2000];
	im plc 2 3
%	im(1, ig.x, ig.y, xtrue, clim), cbar
	im(1, 'mid3', xtrue, clim), cbar
	titlef('x true, z=%g to %g', ig.z(1), ig.z(end))
prompt
end

if ~isvar('proj'), printm 'proj: analytical ellipsoid projection views'
	proj = ellipsoid_proj(cg, ell);

	im(4, 'row', 1, permute(proj, [1 3 2]), 'true helix sinograms'), cbar
prompt
end

if ~isvar('sino'), printm 'sino: rebin helix to fan-beam'
%	f.short = 1;
	f.short = 0; % for rebinned projection views below
	f.itype = 'linear';
%	f.itype = 'nearest';
	[sino, orbits, used] = rebin_helix(cg, ig, proj, ...
		'type', f.itype, 'short', f.short, 'collapse', 0);

	if exist('helix_rebin_mex') == 3 % matlab vs mex
		[sino0, orbits0, used] = rebin_helix(cg, ig, proj, ...
			'use_mex', 0, ...
			'type', f.itype, 'short', f.short, 'collapse', 0);

		equivs(orbits, orbits0)
		equivs(sino0, sino)
		clear sino0 orbits0
	end

	im(5, sino, 'SSRB sinos'), cbar, axis normal
	im(6, used), cbar, axis normal
prompt
end

%% examine the projection views
if false && ~f.short % not possible for short scans!
%	pview = permute(sino, [1 3 2]); % projection views
%	jf_slicer(pview) % awful, understandably because orbits vary!

	% rebin each fan-beam sinogram onto a common orbit_start parallel-beam
	sf = sino_geom('fan', 'dsd', cg.dsd, 'dso', cg.dso, ...
		'ns', cg.ns, 'ds', cg.ds, 'offset_s', cg.offset_s, ...
		'orbit', orbits(1,1), ...
		'na', size(sino, 2)); % # views in fan-beam sinogram, not helix cg.na
	sp = sino_geom('par', 'nb', cg.ns, 'dr', cg.ds/2, ...
		'orbit', 360, 'na', size(sino, 2)); % parallel-beam

	sino2 = zeros(size(sino));
	for iz = 1:ig.nz % each slice
		sf.orbit_start = orbits(iz,2); % sino_geom for this slice!
		sino2(:,:,iz) = rebin_sino(sino(:,:,iz), sf, sp);
	end
%	jf_slicer(sino2)

	pview = permute(sino2, [1 3 2]); % projection views
	jf_slicer(pview) % animate views
return
end

% helical cone-beam reconstruction based on SSRB
if ~isvar('xssrb'), printm 'fbp from ssrb'
	[xssrb, tmp] = fbp_helix_stack(cg, ig, sino, orbits, ...
		'short', f.short, 'window', 'hanning,1.0');
	xssrb = xssrb .* (ig.circ(cg.rmax) > 0); % mask

%	im(2, xssrb, 'SSRB recon', clim), cbar
	im(2, 'mid3', xssrb, 'SSRB recon', clim), cbar
	im(6, tmp, 'after weighting'), cbar, axis normal
	im(3, xssrb - xtrue, 'error', [-500 500]), cbar
prompt
end

if im % profile, for checking HU accuracy
	im subplot 6
	iz = 1:ig.nz;
	ix = ig.nx/2+1; iy = ig.ny/2+1;
	plot(	ig.z, squeeze(xtrue(ix,iy,iz)), 'o-', ...
		ig.z, squeeze(xssrb(ix,iy,iz)), '.--')
	axis([minmax(ig.z)' -200 2100])
	titlef('profile at (ix,iy)=(%g,%g)', ix,iy), xlabel z
end

if 0 % toggle
	fun = @(x) x;
	fun = @(x) jf_mip3(x, 'type', 'mid');
	figure(2)
	im_toggle(fun(xssrb .* ig.circ), fun(xtrue))
return
end

return % below here is comparisons with GH original method

if ~isvar('xssrb_gh'), printm 'fbp from ssrb gh'
	[xssrb_gh, sino_gh] = fbp_helix_gh(cg, ig, proj, 'short', f.short);
	im(3, xssrb_gh, 'SSRB GH recon', clim), cbar
	im(4, sino, 'SSRB GH sinos'), cbar, axis normal
prompt
end

if 0
	figure(2)
%	im_toggle(fun(xssrb .* ig.circ), fun(xtrue), fun(xssrb_gh .* ig.circ), clim)
	im_toggle(fun(xssrb .* ig.circ), fun(xssrb_gh .* ig.circ), clim)
end

if im % profile
	iz = 21;
	iz = 1:ig.nz;
	im(3, xssrb_gh(:,:,iz) - xtrue(:,:,iz), 'GH Error'), cbar

	im subplot 6
	ix = 1:ig.nx; iy = ceil(ig.ny/2); iz = ceil(ig.nz/2);
	iz = round(25/40 * ig.nz);
	plot(ix, xtrue(ix,iy,iz), '-', ix, xssrb_gh(ix,iy,iz), '--')
	axis([1 ig.nx -200 2100])
	titlef('slice %d', iz), xlabel 'ix'
end

% xfdk = easyhelix(cg, ig, li_hat, 'use_mex', has_mex_jf, 'w1cyl', 0);
% xfdk = myFeldkamp(cg, ig, li_hat, 'use_mex', 0);
% im(4, xfdk(:,:,25), 'FDK recon Other'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Stores the result of a classifier
classdef ClassificationResult < Data
    properties (Access = public)
        predictedClasses;
        table;
    end

    properties (Dependent)
        classNames;
        truthClasses;
    end

    methods
        function classNames = get.classNames(obj)
            classNames = obj.table.classNames;
        end

        function set.classNames(obj, classNames)
            obj.table.classNames = classNames;
        end

        function truthClasses = get.truthClasses(obj)
            truthClasses = obj.table.label;
        end

        function set.truthClasses(obj, truthClasses)
            obj.table.label = truthClasses;
        end
    end

    methods (Access = public)
        function obj = ClassificationResult(predictedClasses,table)
            if nargin > 0
                obj.predictedClasses = predictedClasses;
                obj.table = table;
            end
            obj.type = DataType.kClassificationResult;
        end
    end
end


text:
%--- help for generic/load_parameters ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
function S = randsym(n, N)
% Generates random symmetric matrices with normal entries.
%
% function S = randsym(n)
% function S = randsym(n, N)
%
% S is an n-by-n-by-N array where each slice S(:, :, i) for i = 1..N is a
% random symmetric matrix with upper triangular entries distributed
% independently following a normal distribution (Gaussian, zero mean, unit
% variance).
%
% By default, N = 1.
%
% See also: randrot randskew randherm randskewh

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Oct. 23, 2018.
% Contributors:
% Change log:
%       Oct. 23, 2018 (NB):
%           This is not technically necessary for the rotations factory,
%           but it is counter-intuitive to have access to a function called
%           randskew yet not have one for randsym.
%       June 19, 2019 (NB):
%           Now handles the case n = 1 properly.

    if nargin < 2
        N = 1;
    end

    if n == 1
        S = randn(1, 1, N);
        return;
    end

    % Subindices of the (strictly) upper triangular entries of an n-by-n
    % matrix.
    [I, J] = find(triu(ones(n), 1));

    K = repmat(1:N, n*(n-1)/2, 1);

    % Indices of the strictly upper triangular entries of all N slices of
    % an n-by-n-by-N array.
    L = sub2ind([n n N], repmat(I, N, 1), repmat(J, N, 1), K(:));

    % Allocate memory for N random symmetric matrices of size n-by-n and
    % populate each upper triangular entry with a random number following a
    % normal distribution and copy them on the corresponding lower
    % triangular side.
    S = zeros(n, n, N);
    S(L) = randn(size(L));
    S = S + multitransp(S);

    % Now populate the diagonal entries:

    % Subindices of the diagonal entries of an n-by-n matrix.
    [I, J] = find(eye(n));

    K = repmat(1:N, n, 1);

    % Indices of the diagonal entries of all N slices of an n-by-n-by-N
    % array.
    L = sub2ind([n n N], repmat(I, N, 1), repmat(J, N, 1), K(:));

    S(L) = randn(size(L));

end

author: NicolasBoumal

repo: manopt


text:
function [out, param] = pid_loop(y_c, y, dy, kp, ki, kd, limit, Ts, tau, param)

integrator      = param.int;
differentiator  = param.diff;
error_prev      = param.error_prev;
y_c_prev        = param.y_c_prev;

% Update error
error           = y_c - y;

% Update integrator
if isfinite(dy)
    integrator  = integrator + Ts*(error + dy/2);
else
    integrator  = integrator + (Ts/2)*(error + error_prev);
end

% Update differentiator
if isfinite(1/kd)
    if isfinite(dy)
        dy_c = y_c - y_c_prev;          % Compute command diff
        derror = dy_c - dy;             % Compute error diff
        y_c_prev = y_c;                 % Update the command for next step
    else
        derror = error - error_prev;    % Compute error diff
        error_prev = error;             % Update the error for next step
    end
    differentiator = (2*tau-Ts)/(2*tau+Ts)*differentiator...
                        + 2/(2*tau+Ts)*(derror);
end

% Update output before considering saturation
out_unsat = kp*error + ki*integrator + kd*differentiator;
% Check saturation
out = saturate(out_unsat, limit);

% Implement integrator anti-windup
    if ki~=0
        integrator = integrator + Ts/ki * (out - out_unsat);
    end

param.int           = integrator;
param.diff          = differentiator;
param.error_prev    = error_prev;
param.y_c_prev      = y_c_prev;

end

function out = saturate(out_unsat, limit)

    if out_unsat > limit
        out = limit;
    elseif out_unsat < -limit
        out = -limit;
    else
        out = out_unsat;
    end

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function sb_write_materials(filename,values,labels,resolution)
%
% SB_WRITE_MATERIALS writes a .mat file (read by the simbio-vgrid mesher)
% in which each tissue compartment is assigned to an integer value
%
% Use as
%   sb_write_materials(values,labels,weights)
%
% where:
%   values are the integer values assigned to the different compartments
%   labels is a cell with the labels descibing the comartments (default: {'scalp' 'skull' 'brain'})
%   weights are the weights assigned to the different compartments (default: [1.0 1.0 1.0])
%
% $Id$

% Copyright (C) 2011, Cristiano Micheli

if nargin<3 & length(values)==3
  weights = [1. 1. 1.];
  labels  = {'skin', 'skull', 'brain'}; % it assumes this order
elseif nargin<3 & length(values)~=3
  error('You should specify the labels for models with more (or less) than 3 compartments')
end

if isempty(resolution)
  resolution = 1;
end

% open the file and write the header
try
  fid = fopen(filename, 'w');
  fprintf(fid,'%s\n','material bg 1 0 0   2   1.0');
  for i=1:numel(values)
    fprintf(fid,'%s%s %d %d %d %d %0.1f\n','material ',labels{i},resolution,values(i),values(i),1,1);
  end
  fclose(fid);
catch
  disp('Error in writing the file')
  rethrow(ME)
end


text:
classdef prtKernelRbfNeighborhoodScaled < prtKernel
    % prtKernelRbfNeighborhoodScaled  Radial basis function kernel where
    % each entry is scaled according to the distance to neighbors in the
    % training set.
    %
    %  KERNOBJ = prtKernelRbfNeighborhoodScaled Generates a kernel object
    %  implementing a radial basis function.  Kernel objects are widely
    %  used in several prt classifiers, such as prtClassRvm and
    %  prtClassSvm.  RBF kernels implement the following function for 1 x N
    %  vectors x1 and x2:
    %
    %   k(x1,x2) = exp(-sum((x1-x2).^2)./sigma1.^2);
    %
    %  sigma1 is learned based on the neighboord of the data in
    %  feature space.
    %
    %  KERNOBJ = prtKernelRbfNeighborhoodScaled(PROPERTY1, VALUE1, ...) constructs a
    %  prtKernelRbfNeighborhoodScaled object KERNOBJ with properties as specified by
    %  PROPERTY/VALUE pairs. prtKernelRbfNeighborhoodScaled objects have the following
    %  user-settable properties:
    %
    %   neighborhoodPercentile - Quantile distance that is used to define
    %       the neighboorhood. THis is a value between 1 and 100. Smaller
    %       values will make more local kernels. (Default value is 5)
    %
    %   prtKernelRbf objects inherit the TRAIN and RUN methods from prtKernel.
    %
    %   % Example
    %   ds = prtDataGenMoon;                     % Generate a dataset
    %   k1 = prtKernelRbfNeighborhoodScaled;    % Create a prtKernel object with
    %                                            % default value of neighborhoodPercentile
    %   k2 = prtKernelRbfNeighborhoodScaled('neighborhoodPercentile',2); % Create a prtKernel object with
    %                                                                     % the specified value of neighborhoodPercentile
    %
    %   k1 = k1.train(ds); % Train
    %   g1 = k1.run(ds); % Evaluate
    %
    %   k2 = k2.train(ds); % Train
    %   g2 = k2.run(ds); % Evaluate
    %
    %   subplot(2,1,1); imagesc(g1.getObservations);  %Plot the results
    %   subplot(2,1,2); imagesc(g2.getObservations);
    %
    %   See also: prtKernel,prtKernelSet, prtKernelDc, prtKernelRbf, prtKernelDirect,
    %   prtKernelHyperbolicTangent, prtKernelPolynomial,
    %   prtKernelRbfNdimensionScale,

    properties (SetAccess = private)
        name = 'RBF Kernel Scaled'; % RBF Kernel Scaled
        nameAbbreviation = 'RBF'; % RBF
    end

    properties
        neighborhoodPercentile = 5;
        sigmas = []; % The inverse kernel width for each dimension
    end

    methods (Access = protected, Hidden = true)
        function Obj = trainAction(Obj,ds)
            Obj.internalDataSet = ds;

            D = prtDistanceEuclidean(ds,ds);
            D = sort(D,'ascend');

            nPointsAway = min(max(round(size(D,1)*Obj.neighborhoodPercentile/100),1),size(D,1));

            Obj.sigmas = sqrt(D(nPointsAway,:))';
            Obj.sigmas(Obj.sigmas<=0) = 1; % Put weird values to 1

            Obj.isTrained = true;
        end

        function dsOut = runAction(Obj,ds)
            if ~Obj.isTrained
                error('prtKernelRbfNeighboorhoodScaled:run','Attempt to run an untrained kernel; use kernel.train(ds) to train');
            end
            if Obj.internalDataSet.nObservations == 0
                dsOut = prtDataSetClass;
            else
                gram = prtKernelRbfNeighborhoodScaled.kernelFn(ds.getObservations,Obj.internalDataSet.getObservations,Obj.sigmas);
                dsOut = ds.setObservations(gram);
            end
        end
    end

    methods
        function Obj = prtKernelRbfNeighborhoodScaled(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
    end

    methods(Hidden = true)
        function Obj = retainKernelDimensions(Obj,keepLogical)
            Obj.sigmas = Obj.sigmas(keepLogical);
            Obj = retainKernelDimensions@prtKernel(Obj,keepLogical);
        end

        function varargout = plot(obj)
            x = obj.internalDataSet.getObservations;

            if size(x,2) <= 3
                if size(x,2) == 1 && obj.internalDataSet.isLabeled
                    xy = cat(2,x,obj.internalDataSet.getTargets);
                    h = prtPlotUtilScatter(xy, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                else
                    h = prtPlotUtilScatter(x, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                end
            else
                h = nan;
            end

            varargout = {};
            if nargout
                varargout = {h};
            end
        end
    end

    methods (Static, Hidden = true)
        function gram = kernelFn(x,y,sigmas)
            [n1, d] = size(x);
            [n2, nin] = size(y);
            if d ~= nin
                error('size(x,2) must equal size(y,2)');
            end

            %dist2 = prtDistanceLNorm(x,y,2);
            dist2 = repmat(sum((x.^2), 2), [1 n2]) + repmat(sum((y.^2),2), [1 n1]).' - 2*x*(y.');

            if numel(sigmas) == 1
                gram = exp(-dist2/(sigmas.^2));
            else
                gram = exp(-bsxfun(@rdivide,dist2,(sigmas.^2)'));
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function initexcel(pref)

resultsFile = fullfile(pref.results_path, [pref.results_sub '.xlsx'] );
if exist(resultsFile, 'file') == 2
    delete(resultsFile);
end

% then copy the blank excel file from the files to the data folder
sourcefile = [fileparts(mfilename('fullpath')) filesep 'results.xlsx'];
destinationfile = fullfile(pref.results_path, [pref.results_sub '.xlsx']);
if exist(pref.results_path, 'dir') == 0
    mkdir(pref.results_path)
end
copyfile(sourcefile,destinationfile);

author: european-central-bank

repo: BEAR-toolbox


text:
function [resp] = responseCox(X,coeff)
% -------------------------------------------------------------------------
% function [resp] = responseLR(X,coeff)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the multivariable response of an input set of
% features from a given set of logistic regression coefficients.
% -------------------------------------------------------------------------
% INPUTS:
% - X: Matrix of size [nInst X nFeat], specifying the numerical data of the
%      features of the input training data, where 'nInst' refers to the
%      number of instances in X, and 'nFeat' to the number of features in
%      X. Each column is a different feature.
% - coeff: Column vector of size [nFeat+1 X 1] representing the set of
%          logistic regression coefficients computed from
%          drxlr_apply_logistic_regression.m. One coefficient is present
%          for each feature in 'X', in addition to one offset coefficient.
% -------------------------------------------------------------------------
% OUTPUTS:
% - resp: Multivariable response vector of size [nInst X 1].
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

nInst = size(X,1);
nFeat = size(X,2);
resp = zeros(nInst,1);
for j = 1:nFeat
    resp(:) = resp(:) + X(:,j).*coeff(j);
end

end

author: mvallieres

repo: radiomics