text:
function aij = r8cc_get ( m, n, nz_num, colptr, rowind, a, i, j )

%*****************************************************************************80
%
%% R8CC_GET gets a value of a R8CC matrix.
%
%  Discussion:
%
%    It is legal to request entries of the matrix for which no storage
%    was set aside.  In that case, a zero value will be returned.
%
%    The R8CC format is the double precision sparse compressed column
%    format.  Associated with this format, we have an M by N matrix
%    with NZ_NUM nonzero entries.  We construct the column pointer
%    vector COL of length N+1, such that entries of column J will be
%    stored in positions COL(J) through COL(J+1)-1.  This indexing
%    refers to both the ROW and A vectors, which store the row indices
%    and the values of the nonzero entries.  The entries of the
%    ROW vector corresponding to each column are assumed to be
%    ascending sorted.
%
%    The R8CC format is equivalent to the MATLAB "sparse" format,
%    and the Harwell Boeing "real unsymmetric assembled" (RUA) format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero entries.
%
%    Input, integer COLPTR(N+1), indicate where each column's data begins.
%
%    Input, integer ROWIND(NZ_NUM), the row indices.
%
%    Input, real A(NZ_NUM), the nonzero entries.
%
%    Input, integer I, J, the indices of the value to retrieve.
%
%    Output, real AIJ, the value of A(I,J).
%

%
%  Seek sparse index K corresponding to full index (I,J).
%
  k = r8cc_ijk ( m, n, nz_num, colptr, rowind, i, j );
%
%  If no K was found, then be merciful, and simply return 0.
%
  if ( k == -1 )
    aij = 0.0;
  else
    aij = a(k);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function judge = NearZero(near)
% *** BASIC HELPER FUNCTIONS ***
% Takes a scalar.
% Checks if the scalar is small enough to be neglected.
% Example Input:
%
% clear; clc;
% near = -1e-7;
% judge = NearZero(near)
%
% Output:
% judge =
%     1

judge = norm(near) < 1e-6;
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
% MATLAB computation of pulse transformer model - DS Method
% File:  c:\M_files\shortcuts\xfrmrds.m
% 9/19/02; 4/17/04; 2/15/07
%
tic;clc;clear;
K=1e3;pF=1e-12;mH=1e-3;uH=1e-6;ns=1e-9;ps=1e-12; % unit suffixes
%
% Components
%
R1=10;R2=1.5;R3=20*K;R4=1.5;R5=1*K;R6=0.5;R7=1;
C1=20*pF;C2=5*pF;C3=20*pF;L1=1*uH;L2=2*mH;L3=1*uH;
%
% Get A, B, D, & E arrays; this function called only once.
%
Nom=[R1 R2 R3 R4 R5 R6 R7 C1 C2 C3 L1 L2 L3];
[A,B,D,E,I]=tfrmr2(Nom);
%
% * * * * * * * * * * * * Frequency response * * * * * * * * * * * *
%
Ein=10; % Change Ein from 1V to 10V.
%
BF=2;ND=6;PD=50;NP=ND*PD+1;L=linspace(BF,BF+ND,NP);
%
% Since the output is vC3, we dont need the D and E arrays.
% The cv output below is [vC1 vC2 vC3 iL1 iL2 iL3]'
% (a column vector).  Hence we need vC3 or cv(3).
%
for i=1:NP
   F=10^L(i);s=2*pi*F*j;
   cx=(s*I-A)\B*Ein;
%   cy=D*cx+E*Ein; % cy not used
   Vo=abs(cx(3)); % vC3 = Vo
   Vf(i)=20*log10(Vo);
end
%
% * * * * * * * * * * * * * Transient response * * * * * * * * * * *
%
Tx=1/max(max(abs(A)));
disp('Shortest circuit time constant');Tx
%Per=input('Sweep time? (sec)');
% set Sweep time to 200ns = 200e-9 to match Spice run.
Per=200*ns;
kmax=1e5; % kmax increased due to fast time constant Tx
dt = 2*ps
N=6;
%dt=Per/kmax;N=6;
t1=linspace(0,Per,kmax);IV=zeros(N,kmax);
%
% input ramp parameters
%
p=Ein/(5*dt);b=6*dt;pw=5e4*dt;c=pw+6*dt;d=pw+11*dt;
Ea1=ramp1(p,t1(1),dt)-ramp1(p,t1(1),b)-ramp1(p,t1(1),c)+ramp1(p,t1(1),d);
% initialize k = 1
IV(:,1)=B*Ea1*dt;
%
% iterate for k = 2,3,...kmax
%
for k=2:kmax
   Eak=ramp1(p,t1(k),dt)-ramp1(p,t1(k),b)-ramp1(p,t1(k),c)+ramp1(p,t1(k),d);
   IV(:,k)=A*IV(:,k-1)*dt+B*Eak*dt+IV(:,k-1);
end
%
% Plot frequency response
%
subplot(2,1,1)
h=plot(L,Vf,'k');
set(h,'LineWidth',2);
grid on;
axis([BF BF+ND -40 30]);
XT=linspace(BF,BF+ND,7);
set(gca,'xtick',XT);
ylabel('dBV');title('AC Output Vc3');
xlabel('Log Freq(Hz)');
%
% Plot time response
%
subplot(2,1,2)
h=plot(t1/ns,IV(1,:),'k',t1/ns,IV(3,:),'r');
set(h,'LineWidth',2);
axis auto
grid on;ylabel('Volts');title('Transient response, Vc1 & Vc3');
xlabel('nsec');
legend('Vc1','Vc3');

figure(1) % display plot on screen.
%
disp(' ');disp('Execution time in seconds');
ET=toc

author: Sable

repo: mcbench-benchmarks


text:
function planC = getFFDCND(planC,leak,sigma,kernelSize)
%
% 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/>.

%leak = 0.035; %0.018; %leakage in per *100

%sigma = 0.14; %0.02

%kernelSize = 73; %5

indexS = planC{end};

doseNum = 1;

i = find(strcmpi({planC{indexS.structures}.structureName}, 'skin'));

if isempty(i)
    i = find(strcmpi({planC{indexS.structures}.structureName}, 'Body'));
end

if isempty(i)
    i = find(strcmpi({planC{indexS.structures}.structureName}, 'External'));
end

if isempty(i)
    error('Skin structure must be defined to generate a uniform CT.');
else
    i = i(1);
end

structNum = i;

scanSet = 1;

%currentDir = cd;

% [FileName,path] = uigetfile('*.*','Select DICOM Dose File, RD');
%
% if path == 0
%     errordlg('DICOM Dose File Should exist');
%     error('DICOM Dose File Should exist');
% end
%
% cd(path);
%
% fileList = dir;
%
% filesToRun = {};
%
% matches = strmatch('RD', {fileList.name});

% cd(currentDir);
%
% for i = 1:length(matches)
%     [maxDose, dA] = getDICOMMaxDose([path fileList(matches(i)).name]);
%     dM(i) = maxDose;
% end

for i = 1 : planC{7}.FractionGroupSequence.Item_1.NumberOfBeams

    BeamMeterset = planC{7}.FractionGroupSequence.Item_1.ReferencedBeamSequence.(['Item_' num2str(i)]).BeamMeterset;

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

    LS = getDICOMLeafPositions(bs);

    if ~isfield(LS,'manufacturer') && isfield(planC{7},'Manufacturer')
        LS.manufacturer = planC{7}.Manufacturer;
    elseif ~isfield(LS,'manufacturer') && ~isfield(planC{7},'Manufacturer')
        errordlg('Manufacturer is not specified. Using Varian as default.');
        LS.manufacturer = 'Varian';
    end

    [inflMap, xV, yV, colDividerXCoord, rowDividerYCoord, rowLeafPositions,MLCopening] = getLSInfluenceMapFactor(LS,leak,bs.BeamNumber);

    gA = bs.ControlPointSequence.Item_1.GantryAngle;
    iC = bs.ControlPointSequence.Item_1.IsocenterPosition;

    if ~isfield(planC{7},'PatientSetupSequence')
        disp('Patient Setup Missing. Setting to defalt HSF');
        position = 'HSF';
    else
        if ~isfield(planC{7}.PatientSetupSequence,(['Item_' num2str(i)]))
            position = {planC{7}.PatientSetupSequence.(['Item_' num2str(1)]).PatientPosition};
        else
            position = {planC{7}.PatientSetupSequence.(['Item_' num2str(i)]).PatientPosition};
        end
    end

    isocenter = [];

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

    isodistance = bs.SourceAxisDistance/10;

    isocenter = [isocenter.x,isocenter.y, isocenter.z];

    [DoseF, rDepth, xMin, xMax, yMin, yMax, PBSizeX, PBSizeY, maxFlu] = getDose_UnitD(sigma,kernelSize,bs,inflMap, xV, yV,MLCopening);

    clear inflMap

    [xD, yD, zD] = getDoseXYZVals(planC{indexS.dose}(doseNum));

    [xM, yM, zM] = meshgrid(xD, yD, zD);

    coll3V = scan2Collimator([xM(:) yM(:) zM(:)], (gA/360)*2*pi, 0, 0, isocenter, isodistance);

%    distsquared = sepsq(coll3V', [0 0 0]');

    coll3V = (coll3V./ repmat(coll3V(:,3), [1 3])) * -isodistance;

    toKeep = coll3V(:,1) <= xMax + PBSizeX*(kernelSize+1)/2  & coll3V(:,1) >= xMin - PBSizeX*(kernelSize+1)/2 & coll3V(:,2) >= yMin - PBSizeY*(kernelSize+1)/2  & coll3V(:,2) <= yMax + PBSizeY*(kernelSize+1)/2 ;

    pointsToInterpolate = coll3V(toKeep,:);

    clear coll3V distsquared

    radDepthV = getRadiologicalDepth(xM(toKeep), yM(toKeep), zM(toKeep), (gA/360)*2*pi, isocenter, isodistance, structNum, scanSet, planC);

    dV = finterp3(pointsToInterpolate(:,1), pointsToInterpolate(:,2), radDepthV, DoseF, [(xMin - PBSizeX*(kernelSize+1)/2) PBSizeX (xMax + PBSizeX*(kernelSize+1)/2)], [(yMin - PBSizeY*(kernelSize+1)/2)  PBSizeY (yMax + PBSizeY*(kernelSize+1)/2)], rDepth, 0);

    %Should be checked for now turn off
    %dV = dV./distsquared(toKeep);

    dose3D = zeros(size(xM));

    dose3D(toKeep) = dV;

    %dose3D = dose3D/max(dose3D(:));

    %Normalize to max dose
    %dose3D = dM(i)*(dose3D/max(dose3D(:)));

    %Normalize to dose at isocentr
%     [maxDose, dA] = getDICOMMaxDose([path fileList(matches(i)).name]);
%
%     planC{indexS.dose}(end + 1) = planC{indexS.dose}(1);
%     planC{indexS.dose}(end).doseArray = dA;
%
%     Dose_Cl = getDoseAt(length(planC{indexS.dose}),isocenter(1),isocenter(2),isocenter(3),planC);
%
%     planC{indexS.dose}(end) = [];
%
%     planC{indexS.dose}(end + 1) = planC{indexS.dose}(1);
%     planC{indexS.dose}(end).doseArray = dose3D;
%
%     Dose_N = getDoseAt(length(planC{indexS.dose}),isocenter(1),isocenter(2),isocenter(3),planC);
%
%     planC{indexS.dose}(end) = [];

 %   dose3D = (Dose_Cl/Dose_N)*dose3D;

 dose3D = BeamMeterset*dose3D;

    if ~exist('dose')
        dose = dose3D;
    else
        dose = dose + dose3D;
    end

    clear xM yM zM dose3D;

end

mask = find(planC{8}(1).doseArray);

dose1 = zeros(size(planC{8}(1).doseArray));

dose1(mask) = dose(mask);

planC{8}(end+1) = planC{8}(1);

planC{8}(end).doseArray = dose1;

planC{8}(end).fractionGroupID = 'planCheck_Dose';

%m = max(planC{8}(1).doseArray(:));

mM = mean(planC{8}(1).doseArray(:));

%planC{8}(end).doseArray = m*planC{8}(end).doseArray./max(planC{8}(end).doseArray(:));

planC{8}(end).doseArray = mM*planC{8}(end).doseArray./mean(planC{8}(end).doseArray(:));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function z_callBack_virtualTeachPendant_J2_minus()
%% Soft realtime joints control Joint X button callback
% Copyright: Mohammad SAFEEA, 27th-July-2018
global virtual_smartPad_w_joints_command;
disp('Moving J2 in negative direction');
virtual_smartPad_w_joints_command(2)=-1;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% Read NHDR/NRRD files
%
% headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
%
% Inputs:
% * nrrdFileName (char): path to NHDR/NRRD file, either a detached header,
% a detached header pointing to detached data files or a NRRD standalone
% file with header and data included
% * bReadData (bool): set to true to read the data and store it in
% headerInfo.data, set to false to just import the header without the data
%
% Outputs:
% * headerInfo (struct): contains all the fields specified in the file
% header. If data was read it is contained in the 'data' field. That
% structure can be fed to the nhdr_nrrd_write module as is and produce
% valid NHDR/NRRD files.
%
% Format definition:
% http://teem.sourceforge.net/nrrd/format.html
%
% A few supported NRRD features:
% - detached headers with all variants of 'data file:'
% - raw, txt/text/ascii, gz/gzip encodings
% - definition of space and orientation
% - handling of diffusion-weighted MRI data with '<key>:=<value>' lines
%   'modality:=DWMRI', 'DWMRI_b-value:=' and 'DWMRI_gradient_0000:=',
%   'DWMRI_gradient_0001:=', 'DWMRI_gradient_0002:=', etc.
%   (see https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format)
%
%
% Other features:
% - exits cleanly upon error (no file should be left open)
% - informative error messages whenever possible
% - warnings printed to console
% - NHDR/NRRD writer fills out missing fields by inspecting the data
%   whenever possible
%
%
% Unsupported features:
%
% In general, any header field that we were unable to parse is reported
% in a message printed to the console. Specific examples of
% unsupported features include:
%
% - reading data along more than 1 dimension or along a dimension other
%   than the slowest (last) axis specified by the optional <subdim>, as in
%   'data file: <format> <min> <max> <step> [<subdim>]' or 'data file: LIST
%   [<subdim>]'
% - storing all the comments found in headers
% - hex encoding
% - bz2/bzip2 encoding
% - byte skip; can only accept -1 with raw encoding, 0 for all other
%   encodings
% - line skip: can only accept 0
% - checking that field specifications of the form '<field>: <desc>'
%   appear no more than once in the NRRD header (unlike '<key>:=<value>'
%   lines)
%
% Date: July 2018
% Author: Gaetan Rensonnet
%
%
%
% Notes :
%
% This is intended as a more comprehensive and accurate implementation of
% the NRRD format specification than most of the Matlab scripts that have
% been proposed so far. We try to fail graciously when an unsupported
% feature of the NRRD format is encountered. One of our main contributions
% is to propose a writer module which is compatible with the read module,
% in that the output of one can be given as an argument to the other to
% read or produce equivalent NHDR/NRRD files. This is still a version with
% much room for improvement.
%
% The following contributions inspired parts of our Matlab read/write
% modules:
%
% 1. The body of the writer module was pretty much written from scratch but
% the general structure of the reader's main body is based on the Matlab
% functions maReadNrrdHeader.m and maReadNrrdData.m by marc@bwh.harvard.edu
% and kquintus@bwh.harvard.edu (unpublished). Many additions were made and
% a few bugs fixed.
%
% 2. Jeff Mather's NRRD reader
% (http://nl.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader)
% and
% http://jeffmatherphotography.com/dispatches/2012/02/writing-a-file-reader-in-matlab/)
% provided the auxiliary functions:
% - adjustEndian
% - getDatatype, which we renamed nrrd_getMatlabDataType and now throws
% a gracious error if it encounters 'block'-type data,
% - readData was adapted to include a cross-platform fix to delete
% temporary files when using gzip encoding. David Feng's fix used a
% Windows-specific command to delete temporary files
% (https://nl.mathworks.com/matlabcentral/fileexchange/50830-nrrd-format-file-reader).
%
% 3. mdcacio's nrrdWriter
% (https://nl.mathworks.com/matlabcentral/fileexchange/48621-nrrdwriter-filename--matrix--pixelspacing--origin--encoding-)
% provided the auxiliary functions:
% - writeData(): we got rid of the 'unexpected end of input stream when
% attempting to gunzip the file' error when using gzip encoding, which we
% later found had been fixed by Quan Chen independenlty and in a similar
% manner.
% - setDatatype(), which we renamed get_nrrd_datatype() and is the
% reciprocal of getDatatype
%

function headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)

[mainFpath,mainFileName,mainFext] = fileparts(nrrdFileName);

headerInfo = struct();

% default header:
headerInfo.content = mainFileName;        % default value, overwritten if content field is set
headerInfo.data = [];

fidr = fopen(nrrdFileName, 'r');

if (fidr == -1)
    error('ABORT: %s does not exist.\n', nrrdFileName);
end

try
    if ~(strcmpi(mainFext, '.nhdr') || strcmpi(mainFext, '.nrrd'))
        warning('%s looks like a %s file, not a nhdr or nrrd file.\n', nrrdFileName, mainFext );
    end

    % Magic line
    cs  = fgetl(fidr);
    assert(numel(cs) >= 8 && isequal(cs(1:4), 'NRRD'),...
        'Bad signature. First line should be magic line of type "NRRD000X" with X an integer between 1 and 5.'); % FIXME should throw an error if a bad integer is provided
    nrrd_version = sscanf(cs(5:end), '%d');
    if nrrd_version > 5
        error('This reader only supports versions of the NRRD file format up to 5. Detected %d.', nrrd_version)
    end

    % Always optional: defining orientation
    define_orientation = 0;         % internal-use variable

    % Parse header
    while ~feof(fidr)

        cs = fgetl(fidr);         % content string

        if isempty(cs)
            % End of header
            break;
        end

        if foundKeyword( 'CONTENT:', cs )

            headerInfo.content = strtrim( cs( length('CONTENT:')+1:end ) );

        elseif foundKeyword('TYPE:', cs )

            headerInfo.type = strtrim( cs( length('TYPE:')+1:end ) );

        elseif foundKeyword('ENDIAN:', cs )

            headerInfo.endian = strtrim( cs( length('ENDIAN:')+1:end ) );

        elseif foundKeyword('ENCODING:', cs )

            headerInfo.encoding = strtrim( cs( length('ENCODING:')+1:end ) );

        elseif foundKeyword('DIMENSION:', cs )

            headerInfo.dimension = sscanf( cs( length('DIMENSION:')+1:end ), '%i' );

        elseif foundKeyword('SIZES:', cs )

            iSizes = sscanf( cs(length('SIZES:')+1:end), '%i' );
            headerInfo.sizes = iSizes(:)';      % store as row vector

        elseif foundKeyword('KINDS:', cs )

            headerInfo.kinds = extractStringList( cs(length('KINDS:')+1:end) ); % bug fixed with extractStringList where 2 entries are present
            % FIXME: check that axis sizes match each kind according to nrrd standard

        elseif foundKeyword('SPACE:', cs )
            % Starts defining orientation (either space or space dimension,
            % not both)
            define_orientation = 1;

            if isfield(headerInfo, 'spacedimension')
                fprintf(['WARNING nhdr_nrrd_read %s:\n ''space'' field specifier will ' ...
                    'be checked for consistency but will be ignored afterwards ' ...
                    'because ''space dimension'' was specified before.\n'], fopen(fidr));
            end

            tmp_space = strtrim( cs( length('SPACE:')+1:end ) );
            tmp_spacedimension = nrrd_getSpaceDimensions(tmp_space);
            if tmp_spacedimension <= 0
                error('%s: unrecognized ''space'' descriptor ''%s''.', fopen(fidr), tmp_space)
            end

            if isfield(headerInfo, 'spacedimension')
                % internal_spacedimension already set
                if internal_spacedimension ~= tmp_spacedimension
                    error(['%s: ''space'' field specifier implies a spatial ' ...
                            '(world) dimension equal to %d, which differs from ' ...
                            'the ''space dimension'' field specifier set to %d.'],...
                            fopen(fidr), tmp_spacedimension, internal_spacedimension)
                end
                % if no inconsistencies found, just ignore space field
            else
                % Set space info for the first time:
                headerInfo.space = tmp_space;
                internal_spacedimension = tmp_spacedimension;  % for internal use only
            end

        elseif foundKeyword('SPACE DIMENSION:', cs)
            % Starts defining orientation (either space or space dimension,
            % not both)
            define_orientation = 1;
            if isfield(headerInfo, 'space')
                fprintf(['WARNING nhdr_nrrd_read %s:\n ''space dimension'' field specifier ' ...
                    ' will be checked for consistency but will be ignored afterwards ' ...
                    'because ''space'' was specified before.\n'],...
                    fopen(fidr));
            end

            tmp_spacedimension = sscanf( cs( length('SPACE DIMENSION:')+1:end), '%i' );
            if numel(tmp_spacedimension) ~= 1
                error(['%s: ''space dimension'' should be specified as one' ...
                    ' integer number. Found %d element(s) instead.'],...
                      fopen(fidr), numel(tmp_spacedimension))
            end
            if tmp_spacedimension <= 0
                error('%s: ''space dimension'' should be specified as a strictly positive integer (found %d).',...
                                fopen(fidr), tmp_spacedimension)
            end

            if isfield(headerInfo, 'space')
                if tmp_spacedimension ~= internal_spacedimension
                    error(['%s: ''space dimension'' field specifier set to %d, ' ...
                        'which differs from the space (world) dimension implied by the ' ...
                        '''space'' field specifier which is equal to %d.'],...
                            fopen(fidr), tmp_spacedimension, internal_spacedimension)
                end
                % if no inconsistencies found, just ignore space dimension
                % field
            else
                % Set space info for the first time:
                headerInfo.spacedimension = tmp_spacedimension;
                internal_spacedimension = tmp_spacedimension;  % for internal use only
            end

        elseif foundKeyword('SPACE DIRECTIONS:', cs )
            % Required if orientation defined but must come after space or
            % space dimension
            % space directions: <vector[0]> <vector[1]> ... <vector[dim-1]>
            % The format of the <vector> is as follows. The vector is
            % delimited by "(" and ")", and the individual components are
            % comma-separated.
            if ~define_orientation
               error('%s: field specifier ''space directions'' cannot be set before ''space'' or ''space dimension''.',fopen(fidr))
            end

            space_dir_tmp = strtrim(cs(length('SPACE DIRECTIONS:')+1:end));     % remove leading and trailing spaces

            spacedir_vecs = strsplit(space_dir_tmp); %  cell array of strings after split at {' ','\f','\n','\r','\t','\v'}
            SD_data = zeros(internal_spacedimension, internal_spacedimension);

            % Check each vector: either none or (f1,...,f_spaceDim) with fi
            % a floating-point number
            cnt_space_vectors = 0;
            for i = 1:length(spacedir_vecs)
                none_start_index = strfind(lower(spacedir_vecs{i}), 'none');
                if ~isempty(none_start_index)
                    % Axis-specific entry contains substring none
                    if ~strcmpi(spacedir_vecs{i}, 'none')
                        fprintf(['WARNING nhdr_nrrd_read: detected %s instead of ' ...
                            'expected none for axis %d of the per-axis field specifications "space directions:".\n' ...
                            ' There should be no quotation marks, parentheses or any other characters, just plain none.\n'],...
                            spacedir_vecs{i}, i);
                        % Clean none vector specification
                        spacedir_vecs{i} = 'none';
                    end
                else
                    % Axis-specific entry is a numerical vector
                    cnt_space_vectors = cnt_space_vectors + 1;
                    if cnt_space_vectors > internal_spacedimension
                        error(['%s:\n ''space directions'' field specifier: ' ...
                            'number of space vectors detected exceeds space (world)' ...
                            ' dimension, which is equal to %d.'],...
                            fopen(fidr), internal_spacedimension)
                    end
                    btw_parentheses = spacedir_vecs{i}(1) == '(' ...
                                        && spacedir_vecs{i}(end) == ')';
                    axis_vector = regexprep(spacedir_vecs{i}, '[()]', ''); % stripped off all parens
                    if ~btw_parentheses
                        fprintf(['WARNING nhdr_nrrd_read: vector should be delimited ' ...
                                'by parentheses for axis %d of the per-axis field ' ...
                                'specifications "space directions:".\n' ...
                                ' At least one missing parenthesis in ''%s''.\n'],...
                                i, spacedir_vecs{i})
                    end
                    % Clean up by forcing single enclosing parentheses:
                    spacedir_vecs{i} = ['(',  axis_vector,  ')'];
                    % Check vector and extract numerical data
                    vector_entries = strsplit(axis_vector, ',');
                    if length(vector_entries) ~= internal_spacedimension
                        error(['%s:\n vector for data axis %d (space axis %d) of the ' ...
                                'per-axis field specifications "space directions:" should ' ...
                                'contain %d entries corresponding to the space (or world) dimension ' ...
                                'specified in the "space" or "space dimension" field specification.' ...
                                ' Found %d here.'],...
                                fopen(fidr), i, cnt_space_vectors, internal_spacedimension, ...
                                length(vector_entries))
                    end
                    for j = 1:length(vector_entries)
                        vector_entry = sscanf(vector_entries{j}, '%f');
                        if isempty(vector_entry)
                            error(['%s\n in field specification "space directions:",  ' ...
                                'vector for data axis %d (space axis %d) too short. ' ...
                                'Detected %d entries instead of expected %d corresponding to ' ...
                                'space (world) dimension.'],...
                                fopen(fidr), i, cnt_space_vectors, j-1, internal_spacedimension)
                        end
                        SD_data(j, cnt_space_vectors) = vector_entry;
                    end
                end
            end
            % Store array of cleaned up strings
            headerInfo.spacedirections = spacedir_vecs;               % cell array of strings, ideally of the form {'(1,0,0)' '(0,2,0.1)' 'none'  '(0.1,0.1,1.1)'} if internal_spacedimension==3

            % Extract numerical data more leniently:
            SD_data_chk = strrep(space_dir_tmp, 'none', '' );                       % ignore "none" entries
            SD_data_chk = extractNumbersWithout(SD_data_chk, {'(',')',',', '"', ''''} );           % detects numbers up to first non numerical entry
            if numel(SD_data_chk) ~= (internal_spacedimension)^2
                   error(['Expected ''space directions'' to specify a %d-by-%d matrix ' ...
                       '(%d elements in total) in agreement with world space dimension.' ...
                       ' Found %d element(s) instead.\n'],...
                            internal_spacedimension, internal_spacedimension,...
                            (internal_spacedimension)^2, numel(SD_data_chk));
            end

            % Sanity check
            if ~isequal(SD_data_chk(:), SD_data(:))
               error(['%s:\n ''space directions'' field specifier: couldn''t' ...
                   ' read space vectors. Please refer to the NRRD format definition.'],...
                   fopen(fidr))
            end
            headerInfo.spacedirections_matrix = SD_data;
            % Correctness of field specification is checked below after
            % whole header is parsed because we need to be sure that the
            % "dimension" basic field specification was set

        elseif foundKeyword('SPACE UNITS:', cs )
            % Always optional, must come after space or space dimension
            if define_orientation ~= 1
                error('Field specification ''space units'' cannot be specified before ''space'' or ''space dimension''.')
            end

            space_units_tmp = strrep( cs(length('SPACE UNITS:')+1:end), 'none', '');        % ignore none entries
            % FIXME:  ideally, check that the sum of elements including none and
            % " " matches headerInfo.dimension, i.e. the total dimension, as
            % specified in the standard. Standard a bit unclear: should
            % unknown units be specified with "???", "none" or "" ? (quotes
            % seem to be required).

            headerInfo.spaceunits = extract_spaceunits_list( space_units_tmp );

            if length(headerInfo.spaceunits) ~= internal_spacedimension
               error(['Expected ''space units'' to contain %d elements enclosed in double quotes' ...
                   ' to match the ''space'' or ''space dimension'' field but found the following ' ...
                   '%d element(s):\n%s'],...
                   internal_spacedimension, length(headerInfo.spaceunits), sprintf('%s\t',headerInfo.spaceunits{:}))
            end

        elseif foundKeyword('SPACE ORIGIN:', cs )
            % Always optional, must come after space or space dimension

            assert(define_orientation==1, ...
                    sprintf(['%s: field ''space origin'' cannot be specified' ...
                            ' before ''space'' or ''space dimension''.'], ...
                            fopen(fidr)));

            iSO = extractNumbersWithout( cs(length('SPACE ORIGIN:')+1:end), {'(',')',','} );

            assert(numel(iSO) == internal_spacedimension,...
                    sprintf(['%s: expected ''space origin'' to specify a ' ...
                            '%d-element vector to match the ''space'' or ' ...
                            '''space dimension'' field but found %d ' ...
                            'element(s).'],...
                            fopen(fidr),internal_spacedimension, numel(iSO)) );

            headerInfo.spaceorigin = iSO(:);

        elseif foundKeyword('MEASUREMENT FRAME:', cs )
            % Always optional, must come after space or space dimension

            assert(define_orientation==1,...
                    sprintf(['%s: field ''measurement frame'' cannot be ' ...
                            'specified before ''space'' or ''space dimension''.'],...
                            fopen(fidr)));

            measframe_str = strrep( cs(length('MEASUREMENT FRAME:')+1:end), 'none', '');

            iMF = extractNumbersWithout( measframe_str, {'(',')',','} ); % fails if non number entries are not previously removed

            assert(numel(iMF) == (internal_spacedimension)^2,...
                sprintf(['%s: expected ''measurement frame'' to specify a ' ...
                        '%d-by-%d matrix (%d total elements) but found %d ' ...
                        'element(s).'],...
                        fopen(fidr), internal_spacedimension, ...
                        internal_spacedimension, (internal_spacedimension)^2,...
                        numel(iMF)));

            headerInfo.measurementframe = reshape(iMF(:), [internal_spacedimension, internal_spacedimension]);

            % FIXME: this will gladly accept '1 0 0 0 1 0 0 0 1' instead of
            % '(1,0,0) (0,1,0)  (0,0,1)' for instance. But it should have
            % length spacedimension according to standard so this is not too bad.

        elseif foundKeyword('THICKNESSES:', cs )

            sThicknesses = extractStringList( cs(length('THICKNESSES:')+1:end) ); % fixed bug with extractStringList where 2 entries are present
            iThicknesses = [];
            lenThicknesses = length( sThicknesses );
            for iI=1:lenThicknesses
                iThicknesses = [iThicknesses, str2double(sThicknesses{iI}) ];
            end
            headerInfo.thicknesses = iThicknesses;

        elseif foundKeyword('CENTERINGS:', cs )

            headerInfo.centerings = extractStringList( cs(length('CENTERINGS:')+1:end ) ); % fixed bug with extractStringList where 2 entries are present

        elseif foundKeyword('LINE SKIP:', cs) || foundKeyword('LINESKIP', cs)

            if foundKeyword('LINE SKIP:', cs)
                headerInfo.lineskip = sscanf( cs( length('LINE SKIP:')+1:end ), '%d' );
            else
                headerInfo.lineskip = sscanf( cs( length('LINESKIP:')+1:end ), '%d' );
            end
            assert(headerInfo.lineskip >= 0,...
                    sprintf(['Field lineskip or line skip should be greater' ...
                            ' than or equal to zero, detected %d.'], ...
                            headerInfo.lineskip));

        elseif foundKeyword('BYTE SKIP:', cs) || foundKeyword('BYTESKIP:', cs)

            if foundKeyword('BYTE SKIP:', cs)
                headerInfo.byteskip = sscanf( cs( length('BYTE SKIP:')+1:end ), '%d' );
            else
                headerInfo.byteskip = sscanf( cs( length('BYTESKIP:')+1:end ), '%d' );
            end
            assert(headerInfo.byteskip >= -1, ...
                    sprintf(['Field byteskip or byte skip can only take ' ...
                            'non-negative integer values or -1, detected %d.'], ...
                            headerInfo.byteskip));

        elseif foundKeyword('MODALITY', cs )

            headerInfo.modality = strtrim( extractKeyValueString( cs(length('MODALITY')+1:end ) ) );

        elseif foundKeyword('DWMRI_B-VALUE', cs )

            headerInfo.bvalue = str2double( extractKeyValueString( cs(length('DWMRI_B-VALUE')+1:end ) ) );

        elseif foundKeyword('DWMRI_GRADIENT_', cs )

            [iGNr, dwiGradient] = extractGradient(cs(length('DWMRI_GRADIENT_')+1:end ));
            headerInfo.gradients(iGNr+1,:) = dwiGradient;
            % FIXME: make it more general than numbering limited to 4 digits?

        elseif foundKeyword('DATA FILE:', cs ) || foundKeyword('DATAFILE:', cs)
            % This tells us it is a detached header and ends it. 3 possibilities here

            field_value = strtrim( cs(length('DATA FILE:')+1:end) );
            if foundKeyword('DATAFILE:', cs)
                field_value = strtrim( cs(length('DATAFILE:')+1:end) );
            end

            [filelist, LIST_mode, subdim] = extract_datafiles(field_value);

            headerInfo.datafiles = filelist;

            % In LIST mode, filelist is empty and is filled by reading the
            % rest of the header
            if LIST_mode
                datafile_cnt = 0;
                while ~feof(fidr)
                    cs = fgetl(fidr);
                    if isempty(cs)
                        break;
                    end
                    datafile_cnt = datafile_cnt + 1;
                    headerInfo.datafiles{datafile_cnt} = cs;
                end
            end

            % We could technically break the loop here bc data file/datafile is
            % supposed to close a header; however I think it does no harm to
            % keep parsing (in LIST_mode, end of file reached anyways)

        else

            % see if we are dealing with a comment

            csTmp = strtrim( cs );
            if csTmp(1)~='#' && ~strcmp(cs(1:4),'NRRD')
                fprintf('WARNING nhdr_nrrd_read: Could not parse input line: ''%s'' \n', cs );
                % REVIEW: is it better to blindly write it to the output
                % structure in order to be able to write the exact same file
                % later on ?
            end
        end

    end

    % Check for required fields
    % REVIEW: should this be checked only when the data is read? People
    % might be interested in just parsing the header no matter what is
    % in it...)
    assert(isfield(headerInfo, 'sizes'), 'Missing required ''sizes'' field in header');
    assert(isfield(headerInfo, 'dimension'), 'Missing required ''dimension'' field in header');
    assert(isfield(headerInfo, 'encoding'), 'Missing required ''encoding'' field in header');
    assert(isfield(headerInfo, 'type'), 'Missing required ''type'' field in header');

    % Check other cross-field dependencies, etc.
    % TODO: assert lengths of all
    % fields specified, check that kinds and axis sizes match, etc.
    if isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1
        assert( strcmpi(headerInfo.encoding, 'raw'), ...
            sprintf('byte skip value of -1 is only valid with raw encoding. See definition of NRRD File Format for more information.\n'));
    end

    matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);

    %  endian field required if data defined on 2 bytes or more
    if ~(any(strcmpi(headerInfo.encoding,{'txt', 'text', 'ascii'})) || any(strcmpi(matlabdatatype, {'int8', 'uint8'})))
        assert(isfield(headerInfo, 'endian'), 'Missing required ''endian'' field in header');
    else
        if ~isfield(headerInfo,'endian')
            headerInfo.endian = 'little';    % useless but harmless, just for code compatibility because endian field may be accessed later on while reading data
        end
    end

    % Space and orientation information
    if define_orientation
        assert(isfield(headerInfo, 'spacedirections'), ...
                ['Missing field ''space directions'', required if either' ...
                ' ''space'' or ''space dimension'' is set.']);
        if length(headerInfo.spacedirections) ~= headerInfo.dimension
            fprintf(['WARNING nhdr_nrrd_read %s:\n Unexpected format found for ''space directions'' specifier.\n',...
                    ' Expected none entries and vectors delimited by parentheses such as (1.2,0.2,0.25), ', ...
                    'the number of entries being equal to the dimension field specification.\n',...
                    ' See definition of NRRD File Format for more information.\n'], fopen(fidr));
        end
    end

    % TODO: add support for positive line skips
    if isfield(headerInfo, 'lineskip') && headerInfo.lineskip > 0
        assert(isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1, ...
            sprintf(['lineskip option is currently not supported and can ' ...
                    'only be set to zero unless raw encoding is used and ' ...
                    'byte skip is set to -1 (which cancels the effect of ' ...
                    'lineskip altogether).']));
    end

    % TODO: add support for positive byte skips (see line skip)
    if isfield(headerInfo, 'byteskip') && ~strcmpi(headerInfo.encoding,'raw')
        assert( headerInfo.byteskip == 0, ...
            sprintf('byte skip option with non raw encoding is currently not supported and can only be set to zero.\n'));
    end
    if isfield(headerInfo, 'byteskip') && strcmpi(headerInfo.encoding, 'raw')
        assert( headerInfo.byteskip == -1, ...
            sprintf('non-negative byte skip values with raw encoding are currently not supported; byte skip can only be set to -1.\n'));
    end

catch me
    % Clean up before raising error
    fclose(fidr);
    rethrow(me);
end

% Read the data. Detect if data is in the file or in a detached data file
% and check what type of detached file it is if need be.
if bReadData
    N_data_tot = prod(headerInfo.sizes);
    if isfield(headerInfo, 'datafiles')
        % This is a detached header file

        % We no longer need to read the header (closing it before reading
        % all the detached data files is safer)
        fclose(fidr);

        % TODO: we currently only read slices along the slowest axis (i.e.,
        % last coordinate)
        if ~isempty(subdim) && subdim~=headerInfo.dimension
            error(['(detached header): reading data from slices along axis other than the last' ...
                ' (i.e. slowest) one, is currently not supported.\n' ...
                'Last argument [<subdim>] in ''data file'' or ''datafile'' field ' ...
                'should be removed or set equal to %d, which is the detected' ...
                ' ''dimension'' field value, for now.'],...
                        headerInfo.dimension);
        end

        % Read data chunk by chunk from detached data files
        N_data_files = length(headerInfo.datafiles);
        assert(mod(N_data_tot, N_data_files)==0, ...
            sprintf(['Number of detected data files (%d) does not divide total' ...
                    ' number of values contained in data %d obtained from prod(sizes=[%s]).\n'],...
                    N_data_files, N_data_tot, sprintf('%d ',headerInfo.sizes)));

        N_data_per_file = N_data_tot/N_data_files;

        headerInfo.data = zeros(headerInfo.sizes, matlabdatatype);      % specify right type of zeros, otherwise double by default

        for i = 1:N_data_files

            % Check type of detached data file
            [~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i});        % redundant because done above

            data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file;

            if strcmpi(ext_data,'.nhdr')

                error(['datafile %di/%d: nhdr file should not be used as ' ...
                        'detached data file.'], ...
                        i, length(headerInfo.datafiles));

            elseif strcmpi(ext_data, '.nrrd')
                % Detached NRRD file

                bRead_Detached_Data = true;
                tmp_struct = nhdr_nrrd_read(fullfile(mainFpath, ...
                                            [fname_data, ext_data]), ...
                                            bRead_Detached_Data);   % recursive call
                % TODO: check the rest of the structure for
                % inconsistencies with metadata from header
                assert(N_data_files==1 || headerInfo.dimension == tmp_struct.dimension +1, ...
                    sprintf(['Detached header %s: the number of dimensions in detached ' ...
                            'nrrd data file %d/%d (%s) should be one fewer than that of' ...
                            ' the header file.\nDetected %d instead of %d.\nDifferent ' ...
                            'datafile dimensions as specified by the nrrd standard are' ...
                            ' not supported as of now.\n'],...
                            fopen(fidr), i, N_data_files, headerInfo.datafiles{i},...
                            tmp_struct.dimension, headerInfo.dimension-1 ));

                headerInfo.data(data_ind) = tmp_struct.data(:);

                % Store detached data header for last detached file seen,
                % assuming that all are the same. Do not store data though
                % as it would be redundant.
                tmp_struct = rmfield(tmp_struct, 'data');
                headerInfo.detached_header = tmp_struct;

            else
                % e.g., detached .raw file

                fid_data = fopen( fullfile(mainFpath, [fname_data, ext_data]), 'r');
                if( fid_data < 1 )
                    error(['While reading detached header file %s:\ndetached ' ...
                           'data file number %d/%d (%s) could not be opened.'], ...
                           nrrdFileName, i, N_data_files, headerInfo.datafiles{i});
                end
                try
                    tmp_data = readData(fid_data, N_data_per_file, headerInfo.encoding, matlabdatatype);
                    fclose(fid_data);
                catch me_detached
                    fclose(fid_data);
                    rethrow(me_detached);
                end

                tmp_data = adjustEndian(tmp_data, headerInfo.endian);

                headerInfo.data(data_ind) = tmp_data(:);
            end

        end

    else
        % This is a NRRD standalone file: read data directly from it
        try
            headerInfo.data = readData(fidr, N_data_tot, headerInfo.encoding, matlabdatatype);
            fclose(fidr);
        catch me_detached
            fclose(fidr);
            rethrow(me_detached);
        end
        headerInfo.data = adjustEndian(headerInfo.data, headerInfo.endian);
        headerInfo.data = reshape(headerInfo.data, headerInfo.sizes(:)');   % data into expected form. Transpose required by reshape function bc size vectors need to be row vectors
    end
else
    fclose(fidr);
end

end

% -------------------------------------------------------------------%

% ====================================================================
% --- Auxiliary functions -------------------------------------------%
% ====================================================================

function [iGNr, dwiGradient] = extractGradient( st )

% first get the gradient number

iGNr = str2num( st(1:4) ); % FIX numbering limited to 4 digits (synchronize with writer module)

% find where the assignment is

assgnLoc = strfind( st, ':=' );

if ( isempty(assgnLoc) )
    dwiGradient = [];
    return;
else

    dwiGradient = sscanf( st(assgnLoc+2:end), '%f' );

end

end

% Return part of string after :=
function kvs = extractKeyValueString( st )

assgnLoc = strfind( st, ':=' );

if ( isempty(assgnLoc) )
    kvs = [];
    return;
else

    kvs = st(assgnLoc(1)+2:end);

end

end

% Turn space-separated list into cell array of strings.
function sl = extractStringList( strList )
sl = strsplit(strtrim(strList)); % old Matlab file exchange version had a strange bug with lists of length 2
end

% Store in an array the list of numbers separated by the tokens listed in
% the withoutTokens cell array
function iNrs = extractNumbersWithout( inputString, withoutTokens )

auxStr = inputString;

for iI=1:length( withoutTokens )

    auxStr = strrep( auxStr, withoutTokens{iI}, ' ' );

end

iNrs = sscanf( auxStr, '%f' );

end

% Return true if the string keyWord is the beginning of the content string
% cs (ignoring case)
function fk = foundKeyword( keyWord, cs )
lenKeyword = length( keyWord );
fk = (lenKeyword <= length(cs)) && strcmpi( cs(1:lenKeyword), keyWord);
end

% Return cell array of strings with space units in the form {mm, km, m}
% from input of the form "mm " "mm" " m" where undesired blank spaces may
% have crept in. The double quotes are removed in the process but will be
% added by the nhdr/nrrd writer function.
function su_ca = extract_spaceunits_list( fieldValue )
fv_trimmed = strtrim( fieldValue );
su_ca = strsplit(fv_trimmed, '"');                              % units are delimited by double quotes
su_ca = su_ca(~ ( strcmp(su_ca, '') | strcmp(su_ca, ' ') ) );   % remove empty or blank space strings
for i = 1:length(su_ca)
    su_ca{i} = strtrim( su_ca{i} );
end
end

% Extract data file list into a cell array from the datafile or data file
% field in a NHDR file, stripped off its leading and trailing spaces;
% LIST_mode is set to one if a list of files closes the header;
% subdim is empty by default and may be set through either of these field
% specifications:
% data file: <format> <min> <max> <step> [<subdim>]
% data file: LIST [<subdim>]

function [filelist, LIST_mode, subdim] = extract_datafiles( field_string)

field_string = strtrim(field_string);
filelist = {};
LIST_mode = 0;
subdim = [];
if length(field_string)>= 4 && strcmpi(field_string(1:4), 'LIST')
    % LIST mode
    subdim = sscanf(field_string(5:end), '%d'); % assert 1<=dimensions
    LIST_mode = 1;
    return;
else
    % single detached data file or multiple files written in concise form, non LIST mode
    str_lst = strsplit(field_string); % without delimiters specified, splits at any sequence in the set  {' ','\f','\n','\r','\t','\v'}
    if length(str_lst) == 1
        % Single detached data file:
        filelist{1} = str_lst{1};
        % subdim does not make sense here since all of the data is
        % contained in the deatached data file
    elseif any(length(str_lst)==[0, 2, 3]) || length(str_lst)>5
        % Invalid cases (2, 3 or striclty more than 5 entries)
        error(['error in ''data list'' (or ''data list'') field, in non' ...
                ' ''LIST'' mode:\nexpected  ''<filename>'' or ''<format> ' ...
                '<min> <max> <step> [<subdim>]'' but found instead ''%s'' ' ...
                '(%d elements instead of 1, 4 or 5).'],...
                sprintf('%s ',str_lst{:}), length(str_lst));
    else
        % Multiple detached data files with filenames including integral
        % values generated by min:step:max

        str_format = str_lst{1};
        id_min = sscanf(str_lst{2}, '%d');
        id_max = sscanf(str_lst{3}, '%d');
        step = sscanf(str_lst{4}, '%d');

        % check format and indexing values
        assert(step~=0,...
                sprintf(['detached data files with names specified by ' ...
                        'sprintf()-like format:  step should be strictly positive' ...
                        ' or negative, not zero.']));
        if step < 0
            assert(id_min >= id_max, ...
                    sprintf(['detached data files with names specified by ' ...
                            'sprintf()-like format:  when step is <0, min ' ...
                            'should be larger than or equal to max, here we ' ...
                            'found min=%d < max=%d.'],id_min, id_max));
        else
            assert(id_min <= id_max,...
                    sprintf(['detached data files with names specified by ' ...
                            'sprintf()-like format: : when step is >0, min ' ...
                            'should be smaller than or equal to max, here we ' ...
                            'found min=%d > max=%d.'],id_min, id_max));
        end

        % populate data file list
        fileIDs = id_min:step:id_max;
        filelist = cell(length(fileIDs), 1);
        for i = 1:length(fileIDs)  % does not necessarily incude <max>, it cannot be larger than <max>
            expanded_fname = sprintf(str_format, fileIDs(i));
            filelist{i} = expanded_fname;
        end

        if length(str_lst) == 5
            subdim = sscanf(str_lst{5}, '%d');
        end
    end
end

end

% Read data in a column vector from open file with pointer fidIn at the
% beginning of the data to read. The other arguments specify the number of
% elements to read (determined in main body of reader function above), data
% file encoding (specified in NHDR/NRRD header) and the Matlab type the
% data should be converted to (determined in main body of reader function
% above).
function data = readData(fidIn, Nelems, meta_encoding, matlabdatatype)

switch (meta_encoding)
    case {'raw'}
        % This implicitely assumes that byteskip was set to -1 in raw
        % encoding
        data_bytes = (Nelems)*sizeOf(matlabdatatype);   % bytes of useful info to return in data array
        fseek(fidIn,0,'eof');                           % set position pointer to end of file
        tot_bytes_file = ftell(fidIn);                  % current location of the position pointer (I believe as a byte offset from beginning of file)
        fseek(fidIn,tot_bytes_file-data_bytes,'bof');   % make sure position pointer is right before useful data starts

        % Just read binary file (original version of the code)
        data = fread(fidIn, inf, [matlabdatatype '=>' matlabdatatype]);

    case {'gzip', 'gz'}

        tmp = fread(fidIn, Inf, 'uint8=>uint8'); % byte per byte

        tmpBase = tempname(pwd);
        tmpFile = [tmpBase '.gz'];
        fidTmp = fopen(tmpFile, 'wb');  % this creates tmpFile
        assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression')

        try
            fwrite(fidTmp, tmp, 'uint8');
        catch me
            fclose(fidTmp);
            delete(tmpFile);
            rethrow(me);
        end
        fclose(fidTmp);

        try
            gunzip(tmpFile); % this creates tmpBase (tmpFile without the .gz)
        catch me
            delete(tmpFile);
            rethrow(me);
        end
        delete(tmpFile);

        fidTmp = fopen(tmpBase, 'rb');
        % cleaner = onCleanup(@() fclose(fidTmp));
        try
            data = readData(fidTmp, Nelems, 'raw', matlabdatatype); % recursive call
        catch me
            fclose(fidTmp);
            delete(tmpBase);
            rethrow(me);
        end
        fclose(fidTmp);
        delete(tmpBase);

    case {'txt', 'text', 'ascii'}

        data = fscanf(fidIn, '%f');
        data = cast(data, matlabdatatype);

    otherwise
        assert(false, 'Unsupported encoding')
end
% Check number of read elements (sometimes there is an offest of about 1)
assert(Nelems == numel(data),...
        sprintf(['Error reading binary content of %s: detected %d elements' ...
                ' instead of %d announced in header.'],...
                fopen(fidIn), numel(data), Nelems));
end

% Switch byte ordering according to endianness specified in nhdr/nrrd file.
% The metadata should just contain a field 'endian' set to 'little' or
% 'big'.
function data = adjustEndian(data, meta_endian)

[~,~,endian] = computer();

needToSwap = (isequal(endian, 'B') && isequal(lower(meta_endian), 'little')) || ...
    (isequal(endian, 'L') && isequal(lower(meta_endian), 'big'));

if (needToSwap)
    data = swapbytes(data);
end
end

% Size of Matlab's numeric classes in bytes
% See: https://stackoverflow.com/questions/16104423/size-of-a-numeric-class

function numBytes = sizeOf(dataClass)

    % create temporary variable of data type specified
    eval(['var = ' dataClass '(0);']);

    % Use the functional form of whos, and get the number of bytes
    W = whos('var');
    numBytes = W.bytes;

end


text:
function clbvalpl(var1)
    %clbvalpl.m                      A.Allmann
    %
    %   Calculates Freq-Mag functions (b-value) for a catalog
    %   works on cluscat or newclcat

    % Last modification 8/95

    global newclcat cluscat mess bfig backcat
    global ttcat ttm text3 text4 newcat txt1 txt2 txt3

    if var1==1
        if isempty(ttcat)
            if ~isempty(newclcat)  &&  ~isempty(backcat)
                if length(newclcat(:,1))>length(backcat(:,1))
                    newcat=cluscat;
                else
                    newcat=newclcat;
                end
            elseif isempty(newclcat)              %set catalog for bvalue-plot
                newcat=cluscat;
            else
                newcat=newclcat;
            end
        else
            newcat=ttcat;
        end

        [existFlag,figNumber]=figure_exists('b-value curve',1);
        if existFlag
            figure_w_normalized_uicontrolunits(bfig);
            clf reset;
            set(bfig,'visible','off')
        else
            bfig=figure;                     %build figure for plot
        end
        set(bfig,'Units','normalized','NumberTitle','off','Name','b-value curve');
        set(gcf,'pos',[ 0.435  0.8 0.5 0.5])
        matdraw
        uicontrol('Style','Pushbutton',...
            'Callback','myprint',...
            'Units','normalized',...
            'String','Print','Position',[0.02 .68 .08 .05]);

        uicontrol('Style','Pushbutton',...
            'Callback','set(bfig,''visible'',''off'');welcome;done',...
            'Units','normalized',...
            'String','Close','Position',[0.02 .88 .08 .05]);
        uicontrol('Style','Pushbutton',...
            'Callback','clinfo(8)',...
            'Units','normalized',...
            'String','Info','Position',[0.02 .78 .08 .05]);

        set(gcf,'visible','on');
    end
    maxmag = max(newcat.Magnitude);

    % number of mag units
    nmagu = (maxmag*10)+1;

    bval = zeros(1,nmagu);
    % bvalsum = zeros(1,nmagu);
    bvalsum3 = zeros(1,nmagu);

    [bval,xt2] = hist(newcat.Magnitude,(0:0.1:maxmag));
    % bvalsum = cumsum(bval);                        % N for M <=
    bvalsum3 = cumsum(bval(length(bval):-1:1));    % N for M >= (counted backwards)
    xt3 = (maxmag:-0.1:0);

    % backg_be = log10(bvalsum);
    backg_ab = log10(bvalsum3);
    orient tall
    rect = [0.2,  0.3, 0.70, 0.6];           % plot Freq-Mag curves
    axes('position',rect);
    % semilogy(xt2,bvalsum,'om')
    % hold on
    % semilogy(xt2,bvalsum,'-.m')
    semilogy(xt3,bvalsum3,'-.m')
    hold on
    semilogy(xt3,bvalsum3,'om')
    if var1==1
        xlabel('Magnitude ')
        ylabel('Cumulative Number')
    end
    figure_w_normalized_uicontrolunits(mess);
    clf;
    str=['Please select two magnitudes \newlineto be used in the calculation \newlineof the straightline fit.\newlineWait to push Info or Close\newlineafter the selection'];
    te = text(0.01,0.9,str) ;

    set(te,'FontSize',14);
    set(gca,'visible','off');

    disp('Please select two magnitudes to be used in the caclulation of the straight   line fit')

    figure_w_normalized_uicontrolunits(bfig)
    if var1==2
        delete(ttm);delete(text3);delete(text4);delete(txt1);delete(txt2);delete(txt3);
    end
    seti = uicontrol('Units','normal',...
        'Position',[.4 .01 .2 .05],'String','Select Mag1 ');

    pause(1)

    par2 = 0.1 * max(bvalsum3);
    par3 = 0.12 * max(bvalsum3);
    M1b = [];
    M1b = ginput(1);
    tt3=num2str(fix(100*M1b(1))/100);
    text3=text( M1b(1),M1b(2),['|: M1=',tt3] );
    set(seti,'String','Select Mag2');

    pause(0.1)

    M2b = [];
    M2b = ginput(1);
    tt4=num2str(fix(100*M2b(1))/100);
    text4=text( M2b(1),M2b(2),['|: M2=',tt4] );

    pause(0.1)
    delete(seti)

    eqnumber=length(find(newcat.Magnitude>M1b(1) & newcat.Magnitude<M2b(1)));
    tt6=num2str(eqnumber);

    ll = xt3 > M1b(1) & xt3 < M2b(1);
    x = xt3(ll);
    y = backg_ab(ll);
    [p,s] = polyfit(x,y,1);                   % fit a line to background
    f = polyval(p,x);
    f = 10.^f;
    hold on
    ttm= semilogy(x,f,'b');                         % plot linear fit to backg
    set(ttm,'LineWidth',2)
    r = corrcoef(x,y);
    r = r(1,2);
    std_backg = std(y - polyval(p,x));      % standard deviation of fit

    p=-p(1,1);
    p=fix(100*p)/100;
    std_backg=fix(100*std_backg)/100;
    tt2=num2str(std_backg);
    tt1=num2str(p);

    rect=[0 0 1 1];
    h2=axes('position',rect);
    set(h2,'visible','off');
    txt1=text(.16, .18,['B-Value: ',tt1]);
    txt2=text(.16, .12,['Standard Deviation: ',tt2]);
    txt3=text(.16, .06,['Eqs in limits: ',tt6]);
    uicontrol('Style','Pushbutton',...
        'Callback','clbvalpl(2)',...
        'Units','normalized',...
        'String','Repeat','Position',[0.7 .1 .12 .08]);
    welcome;

author: CelsoReyes


text:
function c = digit_inc ( c )

%*****************************************************************************80
%
%% DIGIT_INC increments a digit.
%
%  Discussion:
%
%    0 goes to 1, 1 goes to 2 and so on, but 9 goes back to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character.
%
%    Output, character C, the incremented character.
%
  if ( '0' <= c & c < '9' )
    c = char ( 1 + double ( c ) );
  elseif ( c == '9' )
    c = '0';
  end

  return
end


text:
% generating channels with fliter method and spectrum method
% evaluate autocorrelation functions
%%
clear all
clc
close all

% doppler
fDTs=0.01;
fD=100;
% sampling time
Ts=fDTs/fD;

Ns=1E5;

%% spectrum method
h1=flat_spec(Ns,fD,fDTs);

% autocorrelation
Acn=xcorr(h1,'biased');
l=length(Acn);
% truncate
w=ceil(4/fDTs);
TAc=Acn(ceil(l/2)-w:ceil(l/2)+w);

lt=length(TAc);
t=-Ts*floor(lt/2):Ts:Ts*floor(lt/2);
% theoritical value of autocorrelation
Ac_th=besselj(0,2*pi*fD*t);

figure;
plot(t,real(Ac_th),t,real(TAc),'r');
title('Autocorrelation');
xlabel('time');
legend('theory','spec. method');
grid on

%% filter method
h1=flat_filter(Ns,fD,fDTs);

Acn=xcorr(h1,'biased');
l=length(Acn);
w=ceil(4/fDTs);
TAc=Acn(ceil(l/2)-w:ceil(l/2)+w);

l=length(TAc);
t=-Ts*floor(l/2):Ts:Ts*floor(l/2);
% theory
Ac_th=besselj(0,2*pi*fD*t);

figure;
plot(t,real(Ac_th),t,real(TAc),'r');
title('Autocorrelation');
xlabel('time');
legend('theory','filter method');
grid on


text:
function [score,outputs] = evaluate(CPD, fam, data, ns, cnodes)
% Evaluate evaluate the performance of the classification/regression tree on given complete data
% score = evaluate(CPD, fam, data, ns, cnodes)
%
% fam(i) is the node id of the i-th node in the family of nodes, self node is the last one
% data(i,m) is the value of node i in case m (can be cell array).
% ns(i) is the node size for the i-th node in the whold bnet
% cnodes(i) is the node id for the i-th continuous node in the whole bnet
%
% Output
% score is the classification accuracy (for classification)
%          or mean square deviation (for regression)
%            here for every case we use the mean value at the tree leaf node as its predicted value
% outputs(i) is the predicted output value for case i
%
% Author: yimin.zhang@intel.com
% Last updated: Jan. 19, 2002

if iscell(data)
  local_data = cell2num(data(fam,:));
else
  local_data = data(fam, :);
end

%get local node sizes and node types
node_sizes = ns(fam);
node_types = zeros(1,size(ns,2)); %all nodes are disrete
node_types(cnodes)=1;
node_types=node_types(fam);

fam_size=size(fam,2);
output_type = node_types(fam_size);

num_cases=size(local_data,2);
total_error=0;

outputs=zeros(1,num_cases);
for i=1:num_cases
  %class one case using the tree
  cur_node=CPD.tree.root;  % at the root node of the tree
  while (1)
    if (CPD.tree.nodes(cur_node).is_leaf==1)
      if (output_type==0) %output is discrete
        %use the class with max probability as the output
        [maxvalue,class_id]=max(CPD.tree.nodes(cur_node).probs);
        outputs(i)=class_id;
        if (class_id~=local_data(fam_size,i))
          total_error=total_error+1;
        end
      else   %output is continuous
        %use the mean as the value
        outputs(i)=CPD.tree.nodes(cur_node).mean;
        cur_deviation = CPD.tree.nodes(cur_node).mean-local_data(fam_size,i);
        total_error=total_error+cur_deviation*cur_deviation;
      end
      break;
    end
    cur_attr = CPD.tree.nodes(cur_node).split_id;
    attr_val = local_data(cur_attr,i);
    if (node_types(cur_attr)==0)  %discrete attribute
        % goto the attr_val -th child
        cur_node = CPD.tree.nodes(cur_node).children(attr_val);
    else
        if (attr_val <= CPD.tree.nodes(cur_node).split_threshhold)
          cur_node = CPD.tree.nodes(cur_node).children(1);
        else
          cur_node = CPD.tree.nodes(cur_node).children(2);
        end
    end
    if (cur_node > CPD.tree.num_node)
      fprintf('Fatal error: Tree structure corrupted.\n');
      return;
    end
  end
  %update the classification error number
end
if (output_type==0)
  score=1-total_error/num_cases;
else
  score=total_error/num_cases;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function runpenl(w0)

% runpenl(w0)
% See Article 2.5
% This example integrates the pendulum equation
%   theta"(t) + 0.1*theta'(t)+sin(theta) = 0
% and animates the motion

% The equation of motion of a pendulum is
% described by two first order equations:
% theta'(t)=w; w'(t)=-.1*w-sin(theta). Let
% z=[theta; w], then
% z'(t)=[z(2); 0.2*z(2)-sin(z(1))]

% Create a function defining the differential
% equation in matrix form

% pendulum=inline('[z(2);-sin(z(1))]','t','z');
pendulum=inline(...
   '[z(2);-0.2*z(2)-sin(z(1))]','t','z');

% Choose time limits for the solution
tmax=30; n=501; t=linspace(0,tmax,n);

% Set integration tolerances
opts=odeset('reltol',1e-5,'abstol',1e-5);

% Input angular velocity repeatedly
if nargin==0
  disp(' ')
  disp('   DAMPED PENDULUM MOTION DESCRIBED BY')
  disp(' theta"(t)+0.1*theta''(t)+sin(theta) = 0')
  while 1
    disp(' ')
    disp('Select the angular velocity at the lowest')
    disp('point. Values over 2.42 push the pendulum')
    disp('over the top. Input w0 (or press return')
    disp('to stop)')
    w0=input('w0 = ? > ');
    if isempty(w0)
      disp(' '), disp('All Done'), disp(' '), return
    end

    % Specify the initial conditions
    theta0=0; z0=[theta0;w0];

    % Solve the differential equation using ode45
    [t,th]=ode45(pendulum,t,z0,opts);

    % Plot the angular deflection of the pendulum
    plot(t,180/pi*th(:,1)), xlabel('time')
    ylabel('angular deflection (degrees)')
    title('ANGULAR DEFLECTION OF THE PENDULUM')
    disp(' ')
    disp('Press return to see the animation')
    grid on, shg, pause

    % Animate the motion
    for j=1:n
      T=th(j,1); X=[0;sin(T)]; Y=[0;-cos(T)];
      plot(X,Y,X(2),Y(2),'o','markersize',12)
      axis([-1,1,-1,1]), axis square, axis off
      drawnow, shg, pause(.05)
    end
    clf
  end
  return
else
  theta0=0; z0=[theta0;w0];
  [t,th]=ode45(pendulum,t,z0,opts);
  if max(th(:,1))>pi
     titl='OVER THE TOP';
  else
     titl='LARGE AMPLITUDE MOTION';
  end

  for j=1:n
     T=th(j,1); X=[0;sin(T)]; Y=[0;-cos(T)];
     plot(X,Y,X(2),Y(2),'o','markersize',12)
     axis([-1,1,-1,1]), axis square, axis off
     title(titl), drawnow, shg
  end
  pause(1), close, return
end


text:
function avg = AvgFilter(x)
%
%
persistent prevAvg k
persistent firstRun

if isempty(firstRun)
  k = 1;
  prevAvg = 0;

  firstRun = 1;
end

alpha = (k - 1) / k;
avg   = alpha*prevAvg + (1 - alpha)*x;

prevAvg = avg;
k       = k + 1;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function c=cross7(a,b)
%%CROSS7 Find the cross product of two 7-dimensional vectors. The cross
%        product has the properties:
%        1) cross7(a,b) is a bilinear function of a and b.
%        2) cross7(a,b) is perpendicular to both a and b.
%           Thus, dot(cross7(a,b),a)=0 and dot(cross7(a,b),b)=0
%        3) norm(cross7(a,b))^2=norm(a)^2*norm(b)^2-dot(a,b)^2
%        Just as two forms of the 3D cross product exist (left-handed and
%        right-handed) 480 forms of the 7-dimensional cross product exist.
%        This function implements a cross product defined in terms of an
%        orthonormal basis using asymmetry as in [1].
%
%INPUTS: a A 7 by numVecs matrix of numVecs vectors.
%        b A 7 by numVecs matrix of numVecs vectors.
%
%OUTPUTS: c The 7 by numVecs cross products of a and b, aXb, for each of
%           the vectors in a and b.
%
%As noted in [2], it is only possible for a cross product between two
%vectors to have all of the properties listed above in three and seven
%dimensions.
%
%As noted in 1, the Jacobi identity does not hold. This means that
%cross7(cross7(a,b),c)+cross7(cross7(b,c),a)+cross7(cross7(c,a),b) does
%not necessarily equal zero. However, the following identities can be
%proven:
%cross7(a,b)=-cross7(b,a)
%dot(a,cross7(b,c))=dot(b,cross7(c,a))=dot(c,cross7(a,b))
%cross7(a,cross7(a,b))=dot(a,b)*a-norm(a)^2*b
%cross7(cross7(a,b),cross7(a,c))=cross7(cross7(cross7(a,b),c),a)+cross7(cross7(cross7(b,c),a),a)+cross7(cross7(cross7(c,a),a),b)
%
%REFERENCES:
%[1] P. Lounesto, "Octonians and triality," Advances in Clifford Algebras,
%    vol. 11, no. 2, pp. 191-213, Dec. 2001.
%[2] W. S. Massey, "Cross products of vectors in higher dimensional
%    Euclidean spaces," The American Mathematical Monthly, vol. 90, no. 10,
%    pp. 697-701, Dec. 1983.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numVecs=size(a,2);
c=zeros(7,numVecs);

c(1,:)=+0              +a(2,:).*b(4,:)  +a(3,:).*b(7,:)  -a(4,:).*b(2,:)  +a(5,:).*b(6,:)  -a(6,:).*b(5,:)  -a(7,:).*b(3,:);
c(2,:)=-a(1,:).*b(4,:) +0               +a(3,:).*b(5,:)  +a(4,:).*b(1,:)  -a(5,:).*b(3,:)  +a(6,:).*b(7,:)  -a(7,:).*b(6,:);
c(3,:)=-a(1,:).*b(7,:) -a(2,:).*b(5,:)  +0               +a(4,:).*b(6,:)  +a(5,:).*b(2,:)  -a(6,:).*b(4,:)  +a(7,:).*b(1,:);
c(4,:)=+a(1,:).*b(2,:) -a(2,:).*b(1,:)  -a(3,:).*b(6,:)  +0               +a(5,:).*b(7,:)  +a(6,:).*b(3,:)  -a(7,:).*b(5,:);
c(5,:)=-a(1,:).*b(6,:) +a(2,:).*b(3,:)  -a(3,:).*b(2,:)  -a(4,:).*b(7,:)  +0               +a(6,:).*b(1,:)  +a(7,:).*b(4,:);
c(6,:)=+a(1,:).*b(5,:) -a(2,:).*b(7,:)  +a(3,:).*b(4,:)  -a(4,:).*b(3,:)  -a(5,:).*b(1,:)  +0               +a(7,:).*b(2,:);
c(7,:)=+a(1,:).*b(3,:) +a(2,:).*b(6,:)  -a(3,:).*b(1,:)  +a(4,:).*b(5,:)  -a(5,:).*b(4,:)  -a(6,:).*b(2,:)  +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:
function tsim=get_simtime(host, port)
% TSIM = GET_SIMTIME(HOST, PORT)
%
% Read simulation time from Flightgear over socket
% connection (HOST, PORT).
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);
in = BufferedReader(InputStreamReader(fg_socket.getInputStream()));

% simulation time in seconds
out.println(sprintf('get /sim/time/elapsed-sec\r\n'));
t_sec_str = char(in.readLine());
j = findstr(t_sec_str,'''');
tsim = str2double(t_sec_str(j(1)+1:j(2)-1));

out.close();
in.close();
fg_socket.close();

end


text:
function H=Ray_model(L)
% Rayleigh Channel Model
%  Input : L  : # of channel realization
%  Output: H  : Channel vector
H = (randn(1,L)+j*randn(1,L))/sqrt(2);

author: LyricYang


text:
function sc = sortcomp(compin)
%nejak to moc nefunguje....
compd = dip_image(compin);
comp = double(gaussf(compd,1));

amax = squeeze(max(max(abs(comp))));
amean = squeeze(mean(mean((abs(comp)))));
v = (amax-amean)./(amax+amean);
[a, b] = sort(v, 'descend');
sc = comp(:,:,b);

imax = squeeze(max(max(sc)));
imin = squeeze(min(min(sc)));
msub = find(abs(imin)>imax); %giggest nonzero value negative
sc(:,:,msub) = -sc(:,:,msub); %flip negative to positive

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


function [  ] = salientNumCnt( filesDir, totalNum )
%SALIENTNUMCNT Summary of this function goes here
%   Detailed explanation goes here

    allPoints = zeros(1,totalNum);

    fileExt = '*.txt';

    files = dir(fullfile(filesDir,fileExt));

    for i=0:length(files)-1
        i
        fileCnt = num2str(i);
        fileName = [filesDir, fileCnt, '.txt'];

        file_t = fopen(fileName);
        pointID = fscanf(file_t, '%d');
        fclose(file_t);

        allPoints(sub2ind(size(allPoints), ones(1, length(pointID)), (pointID)')) = 1;

    end

    disp('salient predicted:');
    disp(sum(allPoints));

end


text:
function [Q, R] = gs_m(A)
%GS_M    Modified Gram-Schmidt QR factorization.
%        [Q, R] = GS_M(A) uses the modified Gram-Schmidt method to compute the
%        factorization A = Q*R for m-by-n A of full rank,
%        where Q is m-by-n with orthonormal columns and R is n-by-n.

%        Reference:
%        N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%        Second edition, Society for Industrial and Applied Mathematics,
%        Philadelphia, PA, 2002; sec 19.8.

[m, n] = size(A);
Q = zeros(m,n);
R = zeros(n);

for k=1:n
    R(k,k) = norm(A(:,k));
    Q(:,k) = A(:,k)/R(k,k);
    R(k,k+1:n) = Q(:,k)'*A(:,k+1:n);
    A(:,k+1:n) = A(:,k+1:n) - Q(:,k)*R(k,k+1:n);
end

author: opencobra

repo: cobratoolbox


text:
function [k, gk, kern] = pathKernDiagCompute(kern, x)

% PATHKERNDIAGCOMPUTE Compute diagonal of PATH kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the path kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% RETURN gk : evaluation of ground kernel
% RETURN kern: updated kernel structure
% computed at the given points.
%
% SEEALSO : pathKernParamInit, kernDiagCompute, kernCreate, pathKernCompute
%
% COPYRIGHT : Andrea Baisero, Carl Henrik Ek, 2013

% SHEFFIELDML

maxl=max(cellfun(@(x)size(x,1),x));
kern=pathKernUpdateWMat(kern,maxl);

num=length(x);

gk=cell(1,num);
k=zeros(1,num);
for i=1:num
    li=size(x{i},1);
    w=kern.wmat(1:li,1:li);
    w=.5*(w+w(end:-1:1,end:-1:1));

    gk{i}=kernCompute(kern.gkern,x{i});
    k(i)=sum(sum(w.*gk{i}));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% recompile mex code to have actual version of the NN code

clear functions

eval('mex -I.. -I../.. ../fnearneigh.cpp -O -output fnearneigh');
eval('mex -I.. -I../.. ../fnearneigh.cpp -O -output fnearneigh_partial -DPARTIAL_SEARCH');
eval('mex -I.. -I../.. ../fnearneigh.cpp -O -DBRUTE -output nn_brute');
eval('mex -I.. -I../.. ../fnearneigh.cpp -O -DBRUTE -output nn_brute_partial -DPARTIAL_SEARCH');

eval('mex -I.. -I../.. ../nn_prepare.cpp -O -output nn_prepare');

eval('mex -I.. -I../.. ../nn_search.cpp -O -output nn_search');
eval('mex -I.. -I../.. ../nn_search.cpp -O -output nn_search_partial -DPARTIAL_SEARCH');

eval('mex -I.. -I../.. ../range_search.cpp -O -output range_search');
eval('mex -I.. -I../.. ../range_search.cpp -O -output range_search_partial -DPARTIAL_SEARCH');

eval('mex -I.. -I../.. ../corrsum.cpp -O -output corrsum');
eval('mex -I.. -I../.. ../corrsum.cpp -O -output corrsum_partial -DPARTIAL_SEARCH');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function Xsol=Bouguet2Tsai()
%
% JMM Montiel. 6 sep 2005

% load matlab_calibration.mat
% fc = fc_frontal;
% kc = kc_frontal;
% cc = cc_frontal;
% alpha_c = alpha_c_frontal;
% cam = initialize_cam;
% Huai{
load H:\relaylatest\toolbox_calib\calib_casios\Calib_Results.mat
cam.k1 =    -kc(1);
cam.k2 =    kc(2);
cam.nRows = 720;
cam.nCols = 1280;
cam.Cx =    cc(1);
cam.Cy =    cc(2);
d=0.0112;
cam.f =     d*fc(1);
cam.dx =    d;
cam.dy =    d;
cam.model = 'two_distortion_parameters';

cam.K =     sparse( [ fc(1)   0     cc(1);
                0  fc(2)    cc(2);
                0    0     1] );

% Huai}

% % generate distorted table
% nCols = 320;
% nRows = 240;
% inc=round(min(nCols,nRows)/20)
% [u_grid,v_grid]=meshgrid(-150:inc:nCols+150,-150:inc:nRows+150);
% [nGrid,mGrid]=size(u_grid);
% undistorted=[reshape(u_grid,1, nGrid*mGrid); reshape(v_grid,1,nGrid*mGrid)];
% distorted = distort_bouguet(undistorted,kc,fc,cc)

addpath TOOLBOX_Calib;

Xini = [cam.k1; cam.k2; cam.f; cam.Cx; cam.Cy]

% Image size
nCols = cam.nCols;
nRows = cam.nRows;

% Step for the image grid
inc = round(min(nCols,nRows)/30);
% Image grid
[u_grid,v_grid] = meshgrid(3*inc:inc:nCols-3*inc,3*inc:inc:nRows-3*inc);
[nGrid,mGrid] = size(u_grid);
uv_d = [reshape(u_grid,1, nGrid*mGrid); reshape(v_grid,1,nGrid*mGrid)];

xy_u = normalize(uv_d,fc,cc,kc,alpha_c);
x_distort = [(xy_u(1,:)*fc(1) + cc(1));(xy_u(2,:)*fc(2) + cc(2))];

% Non-linear minimization, minimizes the error in the image grid between the
% Bouguet's calibration model and the one we use in our code. Xsol should
% return the calibration parameters for our code in the following order:
% k1, k2, f, cx, cy. dx and dy are set here to 0.0112
Xsol = lsqnonlin('matlab2tsai_error',...
    Xini,[0;0;0;0;0],[3*cam.k1;3*cam.k2;3*cam.f;nCols;nRows],...
    optimset('LargeScale','on','Display','Iter','TolFun',1e-14),cam,uv_d,x_distort)

repo: ekfmonoslam


text:
function [x,omega] = plotFibre(odf,f,varargin)
% plot odf along a fibre
%
% Syntax
%   plotFibre(odf,f);
%
% Input
%  odf - @SO3Fun
%  f   - @fibre
%
% Options
%  resolution - resolution of each plot
%
% Example
%   odf = SantaFe;
%   f = fibre.gamma(odf.CS,odf.SS)
%   plotFibre(SantaFe,f)
%
% See also
% S2Grid/plot savefigure Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

% get axis
[mtexFig,isNew] = newMtexFigure(varargin{:});

% extract fibre
if isa(f,'fibre')
  [ori,omega] = orientation(f,odf.CS,odf.SS);
elseif isa(varargin{1},'quaternion')
  omega = angle(f(1),f);
  ori = orientation(f,odf.CS,odf.SS);
end

% find loop
delta = angle(ori(2:end),ori(1));
fz = find(delta(:)<1e-2);

% remove values to close together
fz = fz([true;diff(fz)>1]);

if any(fz) && round(numel(delta) / fz(1)) == numel(fz) && ...
    ~check_option(varargin,'comlete')
  ind = 1:find(delta<1e-2,1,'first');
  ori = ori(ind);
  omega = omega(ind);
end

% evaluate the ODF
x = eval(odf,ori,varargin{:});

% plot the fibre
optiondraw(plot(omega./degree,x,'parent',mtexFig.gca),varargin{:});

if isNew
  xlim(mtexFig.gca,[min(omega),max(omega)]./degree);
  ylabel(mtexFig.gca,'Frequency (mrd)')
  xlabel(mtexFig.gca,['misorientation angle (degree) to ' char(ori(1))]);
  drawNow(mtexFig,varargin{:})
end

if nargout == 0, clear x omega; end

author: mtex-toolbox

repo: mtex

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


text:
%GENDATD Generation of 'difficult' normally distributed classes
%
%   A = GENDATD(N,K,D1,D2,LABTYPE)
%
% INPUT
%   N        Number of objects in each of the classes (default: [50 50])
%   K        Dimensionality of the dataset (default: 2)
%   D1       Difference in mean in feature 1 (default: 3)
%   D2       Difference in mean in feature 2 (default: 3)
%   LABTYPE  'crisp' or 'soft' labels (default: 'crisp').
%
% OUTPUT
%   A        Generated dataset
%
% DESCRIPTION
% Generation of a K-dimensional 2-class dataset A of N objects.
% Class variances are very different for the first two dimensions.
% Separation is thereby, for small sample sizes, 'difficult'.
%
% D1 is the difference between the means for the first feature, D2
% is the difference between the means for the second feature. In all
% other directions the means are equal. The two covariance matrices
% are equal with a variance of 1 in all directions except for the
% second feature, which has a variance of 40. The first two feature
% are rotated over 45 degrees to construct a strong correlation.
% Class priors are P(1) = P(2) = 0.5.
%
% If N is a vector of sizes, exactly N(I) objects are generated
% for class I, I = 1,2.
%
% LABTYPE defines the desired label type: 'crisp' or 'soft'. In the
% latter case true posterior probabilities are set for the labels.
%
% SEE ALSO (PRTools Guide)
% DATASETS, PRDATASETS

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

% $Id: gendatd.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function A = gendatd(N,k,d1,d2,labtype)

	if nargin < 5, labtype = 'crisp'; end
	if nargin < 4, d2 = 3; end
	if nargin < 3, d1 = 3; end
	if nargin < 2,  k = 2; end
	if nargin < 1, N = [50 50]; end

	if k < 2,
		error('Number of features should be larger than 1'),
	end
	V = ones(1,k); V(2) = 40; V = sqrt(V);
	p = [0.5 0.5];
	N = genclass(N,p);
	ma = zeros(1,k);
	mb = zeros(1,k); mb(1:2) = [d1, d2];
	A = [randn(N(1),k).*V(ones(1,N(1)),:) + ma(ones(1,N(1)),:); ...
	     randn(N(2),k).*V(ones(1,N(2)),:) + mb(ones(1,N(2)),:)];
	A(:,1:2) = A(:,1:2)*[1 -1; 1 1]./sqrt(2);
	lab = genlab(N);
	A = prdataset(A,lab,'name','Difficult Dataset','prior',p);

	switch labtype
	 case 'crisp'
	  ;
	 case 'soft'
	  U = prdataset([ma(1:2);mb(1:2)],getlablist(A));
	  G = diag(V(1:2));
	  W = nbayesc(U,G);
	  targets = A(:,1:2)*W*classc;
	  A = setlabtype(A,'soft',targets);
	 otherwise
	  error(['Label type ' labtype ' not supported'])
	end

	return

author: marianux

repo: ecg-kit

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


text: %AND Mapping overload


text:
function y = vl_nnspnorm(x, param, dzdy)
%VL_NNSPNORM CNN spaital normalization.
%   Y = VL_NNSPNORM(X, PARAM) computes the spatial normalization of
%   the data X with parameters PARAM = [PH PW ALPHA BETA]. Here PH and
%   PW define the size of the spatial neighbourhood used for
%   nomalization.
%
%   For each feature channel, the function computes the sum of squares
%   of X inside each rectangle, N2(i,j). It then divides each element
%   of X as follows:
%
%      Y(i,j) = X(i,j) / (1 + ALPHA * N2(i,j))^BETA.
%
%   DZDX = VL_NNSPNORM(X, PARAM, DZDY) computes the derivative of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

pad = floor((param(1:2)-1)/2) ;
pad = [pad ; param(1:2)-1-pad] ;

n2 = vl_nnpool(x.*x, param(1:2), 'method', 'avg', 'pad', pad) ;
f = 1 + param(3) * n2 ;

if nargin <= 2 || isempty(dzdy)
  y = f.^(-param(4)) .* x ;
else
  t = vl_nnpool(x.*x, param(1:2), f.^(-param(4)-1) .* dzdy .* x, 'method', 'avg', 'pad', pad) ;
  y = f.^(-param(4)) .* dzdy - 2 * param(3)*param(4) * x .* t ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [Raw_data Im_data] = ssfpReadData(fname, Frames, frac_ks)
% function [Raw_data Im_data] = ssfpReadData(fname, Frames, [frac_ks=5/8]);
%
%	Function reads selected data from a P-file
%
%
%	INPUT:
%		fname   = P file name.
%	    Frames  = Number of temporal frames collected
%       frac_ks = fraction of K-space traversed. [default 5/8]
%
%	EXAMPLE:
%	  [R I] = getrawdata('P00000.7', 101);	% Read full p-file
%
%	Thomas S. John -- Feb 2009.
if notDefined('frac_ks'),   frac_ks = 0.625;        end

verbose = prefsVerboseCheck;

%% load raw data from Pfile
if verbose >= 1
    fprintf('[%s]: Loading Raw Data from %s. \t(%s) \n', mfilename, fname, datestr(now));
end

A = rawloadX(fname);

%% organize/recon data into image space
if verbose >= 1
    fprintf('[%s]: Reconning %s. \t(%s) \n', mfilename, fname, datestr(now));
end

[Raw_data Im_data] = organize_data(A,Frames);

%% construct final images from partial K-space acquisitions
if verbose >= 1
    fprintf('[%s]: Homodyne recon for %s. \t(%s) \n', ...
            mfilename, fname, datestr(now));
end

[RO PE Slices_per_frame Frames] = size(Raw_data);

for f=1:Frames
    for s=1:Slices_per_frame
        Im_data(:,:,s,f) = partial_k(Raw_data(:,:,s,f),frac_ks,RO,PE);
    end
end

% up/down rows flip
Im_data = Im_data(RO:-1:1,:,:,:);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%Program for random permuting the share generation

%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/
%s_athi1983@yahoo.co.in

%Program Description
%This program generates the share 1 & share 2 randomly for each pixel.

function out = generateShare(a,b)

a1 = a(1);
a2 = a(2);
b1 = b(1);
b2 = b(2);

in = [a
      b];
out = zeros(size(in));
randNumber = floor(1.9*rand(1));

if (randNumber == 0)
    out = in;
elseif (randNumber == 1)
    a(1) = a2;
    a(2) = a1;
    b(1) = b2;
    b(2) = b1;
    out = [a
           b];
end


text:
function varargout = uminus(varargin)
%UMINUS   Unary minus for a DISKFUN.

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

[varargout{1:nargout}] = uminus@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function M = nfchoa_order(nls,conf)
%NFCHOA_ORDER maximum order of spatial band-limited NFC-HOA
%
%   Usage: M = nfchoa_order(nls,conf)
%
%   Input parameters:
%       nls     - number of secondary sources
%
%   Output parameters:
%       M       - spherical harmonics order
%       conf    - configuration struct (see SFS_config)
%
%   NFCHOA_ORDER(nls,conf) returns the maximum order of spherical harmonics for
%   the given number of secondary sources in order to avoid spectral repetitions
%   (spatial aliasing) of the dirving signals. The order is
%
%        / nls/2 - 1,   even nls
%   M = <
%        \ (nls-1)/2    odd nls
%
%   for a circular array and
%         _____
%   M = \|nls/2
%
%   for a spherical array.
%
%   See also: driving_function_imp_nfchoa, driving_function_mono_nfchoa
%
%   References:
%       Ahrens (2012) - "Analytic Methods of Sound Field Synthesis", Springer,
%       ISBN 978-3-642-25743-8

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking input parameters =======================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargpositivescalar(nls);
isargstruct(conf);

%% ===== Configuration ===================================================
if conf.nfchoa.order
    M = conf.nfchoa.order;
    return;
end
dimension = conf.dimension;

%% ===== Computation =====================================================
% Get maximum order of spherical harmonics to avoid spatial aliasing
if strcmp('2D',dimension) || strcmp('2.5D',dimension)
    % Ahrens (2012), p. 132
    if isodd(nls)
        M = (nls-1)/2;
    else
        M = nls/2 - 1;
    end
elseif strcmp('3D',dimension)
    % Ahrens (2012), p. 125
    M = floor(sqrt(nls/2));
end

author: sfstoolbox

repo: sfs-matlab


text:
% update model
function [model, id_max] = L1APG_update(img, model, x1, y1, x2, y2)

% initialization
para = model.para;
para.Lambda = model.Lambda;
rel_std_afnv = para.rel_std_afnv;
sz_T = para.sz_T;
nT = para.nT;
T = model.T;
T_mean = model.T_mean;
norms = model.norms;
A = model.A;
alpha = 50;

% construct template
num = numel(x1);
map_aff = zeros(num, 6);
for i = 1:num
    pos = [y1(i), y2(i), y1(i); x1(i) x1(i) x2(i)];
    aff_obj = corners2affine(pos, sz_T);
    map_aff(i,:) = aff_obj.afnv;
end

% draw samples
n_sample = 20;
aff_samples = zeros(n_sample*num, 6);
index_samples = zeros(n_sample*num, 1);
for i = 1:num
    samples = ones(n_sample,1)*map_aff(i,:);
    sc = sqrt(sum(map_aff(i,1:4).^2)/2);
    std_aff	= 0.1 * rel_std_afnv .* [1, sc, sc, 1, sc, sc];
    samples = draw_sample(samples, std_aff); %draw transformation samples from a Gaussian distribution
    aff_samples((i-1)*n_sample+1:i*n_sample, :) = samples;
    index_samples((i-1)*n_sample+1:i*n_sample) = i;
end

aff_samples = [aff_samples; model.map_aff];
n_sample = size(aff_samples,1);
[Y, Y_inrange] = crop_candidates(double(img), aff_samples(:,1:6), sz_T);
if(sum(Y_inrange == 0) == n_sample)
    fprintf('Target is out of the frame!\n');
end

[Y, Y_crop_mean, Y_crop_std] = whitening(Y);	 % zero-mean-unit-variance
[Y, Y_crop_norm] = normalizeTemplates(Y);        % norm one

p = zeros(n_sample,1); % observation likelihood initialization
% reconstruction coefficients
for i = 1:n_sample
    c = APGLASSOup(model.Temp'*Y(:,i), model.Dict, para);
    D_s = (Y(:,i) - A(:,1:nT)*c(1:nT)).^2;  % reconstruction error
    p(i) = exp(-alpha*(sum(D_s)));          % probability w.r.t samples
end

% find the tempalte to be replaced
[~, indW] = min(c(1:nT));

% insert new template
[~, id_max] = max(p(1:end-1));
T(:,indW)	= Y(:,id_max);
T_mean(indW)= Y_crop_mean(id_max);
norms(indW) = Y_crop_std(id_max) * Y_crop_norm(id_max);

[T, ~] = normalizeTemplates(T);
A(:,1:nT)	= T;

%Temaplate Matrix
Temp = A;
Dict = Temp'*Temp;
Temp1 = T*pinv(T);

model.T = T;
model.T_mean = T_mean;
model.norms = norms;
model.A = A;
model.Temp = Temp;
model.Dict = Dict;
model.Temp1 = Temp1;
model.map_aff = 0.2*aff_samples(id_max,:) + 0.8*model.map_aff;
model.aff_samples = ones(para.n_sample,1)*model.map_aff;
model.occlusionNf = 0;

id_max = index_samples(id_max);

author: yuxng


text:
% DEMCMU35VARGPLVMRECONSTRUCTTAYLOR Reconstruct right leg and body of CMU 35.
%
% DESC
%  Load a model learned with vargplvm and reconstruct the missing indices
%  using the Taylor framework developed for fgplvm (angle errors)
% Based on demCmu35fgplvmReconstructTaylor.m

% See also: vargplvmTaylorAngleErrors
% VARGPLVM

randn('seed', 1e5);
rand('seed', 1e5);

% dynUsed=1;
%
% % Get the sequence numbers.
% [Y, lbls] = lvmLoadData('cmu35WalkJog'); %%%
% seq = cumsum(sum(lbls)) - [1:31]; %%%

dataSetName = 'cmu35gplvm';

if ~exist('experimentNo') experimentNo = 1; end

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);

dataSetName = 'cmu35gplvm';

capName = dataSetName;
capName(1) = upper(capName(1));

fileName=['dem' capName 'Vargplvm' num2str(experimentNo)];

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);

%---- Remove some sequences %%%%%%
%{
seqFrom=7;
seqEnd=7;

if seqFrom ~= 1
    Yfrom = seq(seqFrom-1)+1;
else
    Yfrom = 1;
end
Yend=seq(seqEnd);
Y=Y(Yfrom:Yend,:);
seq=seq(seqFrom:seqEnd);
seq=seq - ones(1,length(seq)).*(Yfrom-1);
%}
%-----

origBias = mean(Y);
origScale = 1./sqrt(var(Y));
% %scale = ones(size(scale));
% Y = Y - repmat(origBias, size(Y, 1), 1);
% Ytest = Ytest - repmat(origBias, size(Ytest, 1), 1);
% Y = Y.*repmat(origScale, size(Y, 1), 1);
% Ytest = Ytest.*repmat(origScale, size(Ytest, 1), 1);

% REMOVE LEG TEST
% Indices associated with right leg.
legInd = [8:14];
% Indices associated with upper body.
bodyInd = [21:50];

dt=0.05;
timeStampsTest = ([0:size(Ytest,1)-1].*dt)';

%for experimentNo = 1:3;

% Load saved model.
%load(['dem' capName num2str(experimentNo) '.mat']);
load(fileName); % From demCmu35gplvmVargplvm1.m
fprintf(1,'# Loaded %s\n',fileName);
if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    fprintf(1,'# Dynamics kernel:');
    kernDisplay(model.dynamics.kern);
    model.dynamics.t_star = timeStampsTest;
    kernName = model.dynamics.kern.comp{1}.type;
else
    kernName = 'static';
end
fprintf(1,'# experimentNo=%d\n',experimentNo);
fprintf(1,'# Reconstruction Iters=%d\n',model.reconstrIters);

 model.reconstrIters = 2000; %%%TEMP

startInd = 63;

legErrs = vargplvmTaylorAngleErrors(model, Y, Ytest, startInd, origBias, origScale, legInd, [dataSetName ...
    'Leg'], experimentNo);

save(['dem' capName 'ReconstructLegs' num2str(experimentNo) '.mat'], 'legErrs');

%----
figPath = ['Results/CMU/' num2str(experimentNo) kernName '/Leg/'];
%try
%    saveAllOpenFigures(figPath, [],1)
%catch
%    %
%end
%----

legErrs

bodyErrs = vargplvmTaylorAngleErrors(model, Y, Ytest, startInd, origBias, origScale, bodyInd, [dataSetName ...
    'Body'], experimentNo);

save(['dem' capName 'ReconstructBody' num2str(experimentNo) '.mat'], 'bodyErrs');

figPath = ['Results/CMU/' num2str(experimentNo) kernName '/Body/'];
%try
%    saveAllOpenFigures(figPath, [],1)
%catch
%    %
%end

%----
bodyErrs

% bar(model.kern.comp{1}.inputScales)

%end

repo: vargplvm


text:
% readbvconf() - read Brain Vision Data Exchange format configuration
%                file
%
% Usage:
%   >> CONF = readbvconf(pathname, filename);
%
% Inputs:
%   pathname  - path to file
%   filename  - filename
%
% Outputs:
%   CONF      - structure configuration
%
% Author: Andreas Widmann, University of Leipzig, 2007

% Copyright (C) 2007 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% 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 2 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

% $Id: readbvconf.m 44 2009-11-12 02:00:56Z arnodelorme $

function CONF = readbvconf(pathname, filename)

if nargin < 2
    error('Not enough input arguments');
end

% Open and read file
[IN, message] = fopen(fullfile(pathname, filename));
if IN == -1
    [IN, message] = fopen(fullfile(pathname, lower(filename)));
    if IN == -1
        error(message)
    end;
end
raw={};
while ~feof(IN)
    raw = [raw; {fgetl(IN)}];
end
fclose(IN);

% Remove comments and empty lines
raw(strmatch(';', raw)) = [];
raw(cellfun('isempty', raw) == true) = [];

% Find sections
sectionArray = [strmatch('[', raw)' length(raw) + 1];
for iSection = 1:length(sectionArray) - 1

    % Convert section name
    fieldName = lower(char(strread(raw{sectionArray(iSection)}, '[%s', 'delimiter', ']')));
    fieldName(isspace(fieldName) == true) = [];

    % Fill structure with parameter value pairs
    switch fieldName
        case {'commoninfos' 'binaryinfos'}
            for line = sectionArray(iSection) + 1:sectionArray(iSection + 1) - 1
                splitArray = strfind(raw{line}, '=');
                CONF.(fieldName).(lower(raw{line}(1:splitArray(1) - 1))) = raw{line}(splitArray(1) + 1:end);
            end
        case {'channelinfos' 'coordinates'}
            for line = sectionArray(iSection) + 1:sectionArray(iSection + 1) - 1
                splitArray = strfind(raw{line}, '=');
                tmp = str2num(strtrim(raw{line}(3:splitArray(1) - 1)));
                if isempty(tmp)
                    tmp = NaN; end
                CONF.(fieldName)(tmp) = {raw{line}(splitArray(1) + 1:end)};
            end
        case {'markerinfos'} % Allow discontinuity for markers (but not channelinfos and coordinates!)
            for line = sectionArray(iSection) + 1:sectionArray(iSection + 1) - 1
                splitArray = strfind(raw{line}, '=');
                tmp = str2num(strtrim(raw{line}(3:splitArray(1) - 1)));
                if isempty(tmp)
                    tmp = NaN; end
                CONF.(fieldName)(line - sectionArray(iSection), :) = {raw{line}(splitArray(1) + 1:end) tmp};
            end
            if ~all(1:size(CONF.(fieldName), 1) == [CONF.(fieldName){:, 2}])
                warning('Marker number discontinuity.')
            end
        case 'comment'
            CONF.(fieldName) = raw(sectionArray(iSection) + 1:sectionArray(iSection + 1) - 1);
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function [timestamps] = ufmf_read_timestamps(header,t0,t1)

FRAME_CHUNK = 1;

fp = header.fid;

if nargin < 2,
  t0 = 1;
end
if nargin < 3,
  t1 = header.nframes;
end

timestamps = nan(1,t1-t0+1);
for t = t0:t1,
  fseek(fp,header.frame2file(t),'bof');

  % read in the chunk type: 1
  chunktype = fread(fp,1,'uchar');
  if chunktype ~= FRAME_CHUNK,
    error('Expected chunktype = %d at start of frame, got %d',FRAME_CHUNK,chunktype);
  end
  % read in timestamp: 8
  timestamps(t-t0+1) = fread(fp,1,'double');
end

author: kristinbranson

repo: JAABA