text:
function r = randp(lambda)
% randp(lambda) returns Poisson distributed Vector with mean lambda

try
  r = poissrnd(lambda);
  return
end

lambda   = reshape(lambda,1,[]);
aktiv    = find(lambda > 1e-10);
ll       = log(lambda(aktiv));
r(aktiv) = rand(1,length(aktiv));
r(find(lambda < 1e-10)) = 0;

i = 0;
while length(aktiv) > 0
		r(aktiv) = r(aktiv) - exp(i*ll - lambda(aktiv) - gammaln(i+1));
		ind = find(r(aktiv)<=0);
		r(aktiv(ind)) = i;
		aktiv(ind) = [];
		ll(ind) = [];
		i = i + 1;
end


text:
function rts_retrieveSolarPanelsAndStore(~, ~,tcpipClient, vesselGUID, dblsTcpipClient)
%rts_retrieveSolarPanelsAndStore Summary of this function goes here
%   Detailed explanation goes here

    try
        data = get(dblsTcpipClient,'UserData');
        numDbls = data{1}(2);
        if(~isempty(numDbls))
            rts_retrieveDataAndStore(tcpipClient, numDbls, 'd');
        end
    catch ME
%         disp(ME.message);
%         for k=1:length(ME.stack)
%             ME.stack(k)
%         end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function detX = det_2x2(X)

[d1,d2,~] = size(X);
assert(d1==2 && d2==2);

detX = X(1,1,:).*X(2,2,:)-X(1,2,:).*X(2,1,:);

author: kristinbranson

repo: JAABA


text:
%int2zstr - Convert integer to zero-padded string.
%
%  USAGE
%
%    s = int2zstr(n,digits)
%
%    n              integer to convert
%    digits         optional total number of digits
%

% Copyright (C) 2012 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function s = int2zstr(n,digits)

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

% Convert to string
s = int2str(n);

% Pad if necessary
if nargin < 2, return; end
l = digits-length(s);
if l <= 0, return; end
s = [ones(1,l)*'0' s];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function o = setField(o, value)
% o = setField(o, source)
% For each field of o such that it is a field in source,
% we set o.(field) = source.(field)

key = fieldnames(o);
for i = 1:length(key)
   if isfield(value, key{i}) || isprop(value, key{i})
      o.(key{i}) = value.(key{i});
   end
end

end

author: opencobra

repo: cobratoolbox


text:
function [positionSLR,marksSLR,ind,messages]=SLR2(marks,messages)
% post processing rule for delineation
%
% [positionSLR,marksSLR,ind,messages]=SLR2(marks,messages)
%
% INPUT:
% marks: structure with a matrix of annotations in samples; in each field
%        each matrix of annotations has one line per beat and 10 colunms
%        corresponding respectively to P onset, P peak, P end, QRS onset, QRS
%        main peak, QRS end, T onset, T peak, T prima, T end
% messages.setup.SLR.freq: sampling frequency of the recording.
% messages.setup.SLR.QRSlimit: admited diference between marks of the same beat in ms (optional)
%           by default QRSlimit= 250 ms.
% messages.setup.SLR.k: number of neighbours to consider in SLR rule (optional)
%    by default k=3
% messages.setup.SLR.delta: neighbourhood to consider in SLR rule in ms (optional)
%    by default delta = [4 nan 4 12 nan 10 12 nan nan 12] ms for 10 marks
%    or  delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12] ms for 15 marks
% messages.setup.wavedet.refrper: minimum time between two admisible beats, by default 257ms
% messages.setup.SLR.tolerance = time in ms to define each individual QRS wave, by default 10ms
%
% OUTPUT:
%positionSLR: struture with multilead based on marks obtained by SLR rule
% marksSLR: matrix multilead based on marks obtained by SLR rule
% ind: indexes of the initial column vectors for each set, corresponding to the sincronized beats
%      each line of ind corresponds to a beat and each column to a set
%      if a beat is absent of a set
% messages.errors: errors during the processing
% messages.errors_desc:errors description
% messages.status: sucess of the procesing (1) ou fail (0)
% messages.status: parameters considered
% Last update: 27JAN2012
positionSLR=[];
marksSLR=[];
ind=[];
if nargin<2
    messages.errors= 'Fatal error in SLR.';
    messages.errors_desc= 'Inputs marks and messages.setup.SLR.freq required';
    messages.status=0;
    return
else
    if size(marks{1},2)~=10 && size(marks{1},2)~=15
        messages.errors= 'Fatal error in SLR.';
        messages.errors_desc= 'Input marks should include 10 or 15 marks.';
        messages.status=0;
        return
    end
    if ~isfield(messages,'setup') || ((~isfield(messages.setup,'SLR') ||...
            (~isfield(messages.setup.SLR,'freq') && ~isfield(messages.setup.SLR,'fa')))  &&...
            (~isfield(messages.setup.wavedet,'wavedet') || ( isfield(messages.setup.wavedet,'wavedet') &&...
            ~isfield(messages.setup.wavedet,'freq')))) %
        messages.errors= 'Fatal error in SLR.';
        messages.errors_desc= 'Inputs messages.setup.SLR.freq required';
        messages.status=0;
        return
    else
        try
            fa=messages.setup.SLR.freq;
        catch me
            me.message = 'Not exist messages.setup.SLR.freq';
            try
                fa=messages.setup.SLR.fa;
            catch me
                me.message = 'Not exist messages.setup.SLR.fa';
                fa= messages.setup.wavedet.freq;
            end
        end
    end
end
messages.status=1;
if ~isfield(messages,'errors')
    messages.errors=[];
end
if ~isfield(messages,'errors_desc')
    messages.errors_desc=[];
end
if ~isfield(messages,'warnings')
    messages.warnings=[];
end
if ~isfield(messages,'status')
    messages.status=1;
end
if ~isfield(messages.setup.SLR,'k')
    messages.setup.SLR.k=min(3,size(marks{1},1)-1);
end
if ~isfield(messages.setup.SLR,'QRSlimit')
    messages.setup.SLR.QRSlimit=250;
end
if ~isfield(messages.setup.SLR,'delta')
    if size(marks{1},2)==10
        messages.setup.SLR.delta = [4 nan 4 12 nan 10 12 nan nan 12];
    elseif size(marks{1},2)==15
        messages.setup.SLR.delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12];
    end
end

QRSlimit= messages.setup.SLR.QRSlimit*fa/1000;
delta= messages.setup.SLR.delta*fa/1000;
k= messages.setup.SLR.k;
if isempty(k) || ~isnumeric(k) || isnan(k)
    k=min(3,size(marks{1},1)-1);
    messages.setup.SLR.k=k;
end
if isempty(QRSlimit) || ~isnumeric(QRSlimit) || isnan(QRSlimit)
    messages.setup.SLR.QRSlimit=250;
    QRSlimit=messages.setup.SLR.QRSlimit*fa/1000;
end
if isempty(delta) | size(delta)~=size(marks{1},2) %#ok<OR2>
    if size(marks{1},2)==10
        messages.setup.SLR.delta = [4 nan 4 12 nan 10 12 nan nan 12];
    elseif size(marks{1},2)==15
        messages.setup.SLR.delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12];
    end
    delta= messages.setup.SLR.delta*fa/1000;
end

if ~isfield(messages.setup.SLR,'refrper')
    if isfield(messages.setup,'wavedet') && isfield(messages.setup.wavedet,'refrper')
        messages.setup.SLR.refrper=messages.setup.wavedet.refrper*fa/1000;
    else
        messages.setup.SLR.refrper=275*fa/1000;
    end
end
if ~isfield(messages.setup.SLR,'tolerance')
    messages.setup.SLR.tolerance=10*fa/1000;
end

if isstruct(marks)
    leads = fieldnames(marks);
    messages.setup.SLR.nleads=length(leads);
else
   messages.setup.SLR.nleads = size(marks,2);
end

if length(delta)==10
    messages.setup.SLR.marks_desc={'(','P',')','(','QRS',')','(','T','T''',')'};
    peaks=[2 5 8 9];
    onsets=[1 4 6];
    ends=[3 6 10];
    refmark=5;
elseif length(delta)==15
    messages.setup.SLR.marks_desc={'(','P','P''',')','QRS_first_peak','(','Q','R','S','R''',')','(','T','T''',')'};
    peaks=[2 3 8 13 14];
    onsets=[1 6 12];
    ends=[4 11 15];
    refmark=[8 9 7 10];%RUTE % porque la onda R puede ni existir en una derivacion especifica!
end
nmarks=length(delta);
beats=0;
refmark_n=refmark;

for g=1:messages.setup.SLR.nleads
    if isstruct(marks)
        f = getfield(marks, char(leads(g)));  %#ok<GFLD>
    else
        f = marks{g};
    end
    if length(f)> size(beats,1)
        beats((size(beats,1)+1):length(f),1:g)=zeros(length((size(beats,1)+1):length(f)),g);
    end
    for n=1:size(f,1)
        refmark(n)=refmark_n(find(~isnan(f(n,refmark_n)),1,'first')); %#ok<AGROW>
        beats(n,g)=f(n,refmark(n));  %#ok<AGROW>
    end
end

beats(beats==0)=NaN;
[ind,messages1]=sincronizabeats(beats,QRSlimit);
if messages1.status == 0
    return;
end
ind(sum(~isnan(ind),2)<k+1,:) = [];

beats=NaN*ones(size(ind));
for g=1:size(ind,2)
    if isstruct(marks)
        f = getfield(marks, char(leads(g)));  %#ok<GFLD>
    else
        f = marks{g};
    end
    beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),5);
end
ind(find(diff(nanmedian(beats,2)) < messages.setup.SLR.refrper)+1,:) = [];

nbeats=size(ind,1);
marksSLR=NaN.*ones(nbeats,nmarks);
marksSLR_lead=NaN.*ones(nbeats,nmarks);
for mark=peaks
    beats=NaN*ones(size(ind));
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g)));  %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end
    marksSLR(:,mark)=nanmedian(beats,2);
end

if  nmarks==15;
    mark=7:10;
    QRS{1}=NaN*ones(size(ind));
    QRS{2}=QRS{1};
    QRS{3}=QRS{1};
    QRS{4}=QRS{1};
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g)));  %#ok<GFLD>
        else
            f = marks{g};
        end
        QRS{1}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(1));  %#ok<AGROW>
        QRS{2}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(2)); %#ok<AGROW>
        QRS{3}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(3)); %#ok<AGROW>
        QRS{4}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(4)); %#ok<AGROW>
    end
    for beat=1:nbeats
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        %look for small R waves that should go with Q waves
        %RS that is an inverted QR or RSR' that is an inverted QRS
        ii=find(AA(2,:)-nanmedian(AA(2,:))<-messages.setup.SLR.tolerance);
        QRS{1}(beat,ii)=QRS{2}(beat,ii); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{3}(beat,ii); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{4}(beat,ii); %#ok<AGROW>
        QRS{4}(beat,ii)=NaN; %#ok<AGROW>
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        %look for Q waves that should go with R waves
        %QS complexes or QR that is an inverted RS
        ii=find(AA(1,:)-nanmedian(AA(1,:))>messages.setup.SLR.tolerance);
        QRS{4}(beat,ii)=QRS{3}(beat,ii); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{2}(beat,ii); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{1}(beat,ii); %#ok<AGROW>
        QRS{1}(beat,ii)=NaN; %#ok<AGROW>
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        %look for R' waves that should go with S waves
        ii=find(AA(4,:)-nanmedian(AA(4,:))<-messages.setup.SLR.tolerance);
        QRS{1}(beat,ii)=QRS{2}(beat,ii); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{3}(beat,ii); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{4}(beat,ii); %#ok<AGROW>
        QRS{4}(beat,ii)=NaN;  %#ok<AGROW>
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        marksSLR(beat,mark)=nanmedian(AA,2)';
    end
end
if size(marksSLR,2) == 15
    positionSLR.P = marksSLR(:,2);
    positionSLR.Pprima = marksSLR(:,3);
    positionSLR.qrs = marksSLR(:,8);
    positionSLR.QRSon = marksSLR(:,6);
    positionSLR.Q = marksSLR(:,7);
    positionSLR.R = marksSLR(:,8);
    positionSLR.S = marksSLR(:,9);
    positionSLR.Rprima = marksSLR(:,10);
    positionSLR.T = marksSLR(:,13);
    positionSLR.Tprima = marksSLR(:,14);
else
    positionSLR.P = marksSLR(:,2);
    positionSLR.qrs = marksSLR(:,5);
    positionSLR.T = marksSLR(:,8);
    positionSLR.Tprima = marksSLR(:,9);
end
for mark=onsets, %onsets
    beats=NaN.*ones(nbeats,nmarks);
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end
    aux= sum(~isnan(beats),2);
    for beat=1:nbeats
        if aux(beat)>k, % if the mark was found in more than k of the leads
            kk = beats(beat,:);
            kk(isnan(kk))=[];
            [kk kki] =sort(kk); % sorted the mark found in the leads in beat
            while (length(kk)>k && kk(k+1)>(kk(1)+delta(mark))), %eliminate the first mark until kk satisfy the rule or length(kk)<k
                kk(1)=[];
                kki(1)=[];
            end
            %if length(kk)<k there is not a multilead mark
            if length(kk)==k
                marksSLR(beat,mark)=nan;
            else
                marksSLR(beat,mark)=kk(1);
                marksSLR_lead(beat,mark)=kki(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end
end
positionSLR.Pon = marksSLR(:,onsets(1));
positionSLR.QRSon = marksSLR(:,onsets(2));
positionSLR.Ton = marksSLR(:,onsets(3));
positionSLR.Pon_lead = marksSLR_lead(:,onsets(1));
positionSLR.QRSon_lead = marksSLR_lead(:,onsets(2));
positionSLR.Ton_lead = marksSLR_lead(:,onsets(3));

for mark=ends, % ends
    beats=NaN.*ones(nbeats,nmarks);
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end
    aux= sum(~isnan(beats),2);
    for beat=1:nbeats
        if aux(beat)>k, %
            kk = beats(beat,:);
            kk(isnan(kk))=[];
            [kk kki] =sort(kk,'descend');
           % kk =flipud(sort(kk)')';
            while (length(kk)>k && kk(k+1)<(kk(1)-delta(mark))),
                kk(1)=[];
                kki(1)=[];
            end
            if length(kk)==k
                marksSLR(beat,mark)=nan;
            else
                marksSLR(beat,mark)=kk(1);
                marksSLR_lead(beat,mark)=kki(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end
end
positionSLR.Poff = marksSLR(:,ends(1));
positionSLR.QRSoff = marksSLR(:,ends(2));
positionSLR.Toff = marksSLR(:,ends(3));
positionSLR.Poff_lead = marksSLR_lead(:,ends(1));
positionSLR.QRSoff_lead = marksSLR_lead(:,ends(2));
positionSLR.Toff_lead = marksSLR_lead(:,ends(3));

author: marianux

repo: ecg-kit

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


text:
function [varargout] = ndgrid(varargin)

% NDGRID Generation of arrays for N-D functions and interpolation.
% [X1,X2,X3,...] = NDGRID(x1,x2,x3,...) transforms the domain
% specified by vectors x1,x2,x3, etc. into arrays X1,X2,X3, etc. that
% can be used for the evaluation of functions of N variables and N-D
% interpolation.  The i-th dimension of the output array Xi are copies
% of elements of the vector xi.
%
% [X1,X2,...] = NDGRID(x) is the same as [X1,X2,...] = NDGRID(x,x,...).
%
% For example, to evaluate the function  x2*exp(-x1^2-x2^2-x^3) over the
% range  -2 < x1 < 2,  -2 < x2 < 2, -2 < x3 < 2,
%
%     [x1,x2,x3] = ndgrid(-2:.2:2, -2:.25:2, -2:.16:2);
%     z = x2 .* exp(-x1.^2 - x2.^2 - x3.^2);
%     slice(x2,x1,x3,z,[-1.2 .8 2],2,[-2 -.2])
%
% NDGRID is like MESHGRID except that the order of the first two input
% arguments are switched (i.e., [X1,X2,X3] = NDGRID(x1,x2,x3) produces
% the same result as [X2,X1,X3] = MESHGRID(x2,x1,x3)).  Because of
% this, NDGRID is better suited to N-D problems that aren't spatially
% based, while MESHGRID is better suited to problems in cartesian
% space (2-D or 3-D).
%
% This is a drop-in replacement for the MATLAB version in elmat, which is
% relatively slow for big grids. Note that this function only works up
% to 5 dimensions
%
% See also MESHGRID, INTERPN.

% Copyright(C) 2010, Jan-Mathijs Schoffelen, DCCN
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin==0
  ft_error('MATLAB:ndgrid:NotEnoughInputs', 'Not enough input arguments.');
end
if nargin==1, varargin = repmat(varargin,[1 max(nargout,2)]); end

ndims = numel(varargin);
switch ndims
case 2
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';

  varargout{1} = x(:, ones2);
  varargout{2} = y(ones1, :);
case 3
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};

  varargout{1} = x(:, ones2, ones3);
  varargout{2} = y(ones1, :, ones3);
  varargout{3} = z(ones1, ones2, :);
case 4
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));
  ones4 = ones(1,numel(varargin{4}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};
  xx   = zeros(1,1,1,numel(varargin{4}));
  xx(:) = varargin{4};

  varargout{1} = x(:, ones2, ones3, ones4);
  varargout{2} = y(ones1, :, ones3, ones4);
  varargout{3} = z(ones1, ones2, :, ones4);
  varargout{4} = xx(ones1, ones2, ones3, :);
case 5
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));
  ones4 = ones(1,numel(varargin{4}));
  ones5 = ones(1,numel(varargin{5}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};
  xx   = zeros(1,1,1,numel(varargin{4}));
  xx(:) = varargin{4};
  yy   = zeros(1,1,1,1,numel(varargin{5}));
  yy(:) = varargin{5};

  varargout{1} = x(:, ones2, ones3, ones4, ones5);
  varargout{2} = y(ones1, :, ones3, ones4, ones5);
  varargout{3} = z(ones1, ones2, :, ones4, ones5);
  varargout{4} = xx(ones1, ones2, ones3, :,ones5);
  varargout{5} = yy(ones1, ones2, ones3, :,ones5);
otherwise
  ft_error('this version of ndgrid supports inputs up to 5 dimensions');
  %call the ndgrid from elmat
  %FIXME this has to be done
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [t,s]=zerocros(x,m)
%ZEROCROS finds the zeros crossings in a signal [T,S]=(X,M)% find zero crossings in a signal
% Inputs:  x = input waveform
%          m = mode string containing:
%              'p' - positive crossings only
%              'n' - negative crossings only
%              'b' - both (default)
%              'r' - round to integer values
%
% Outputs: t = sample positions of zero crossings (not necessarily integers)
%          s = estimated slope of x at the zero crossing
%
% This routine uses linear interpolation to estimate the position of a zero crossing
% A zero crossing occurs between x(n) and x(n+1) iff (x(n)>=0) ~= (x(n+1)>=0)

% Example: x=sin(2*pi*(0:1000)/200); x(1:100:1001)=0; zerocros(x);
% Note that we get a zero crossing at the end but not at the start.

%	   Copyright (C) Mike Brookes 2003-2015
%      Version: $Id: zerocros.m 6077 2015-04-20 07:08:39Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
    m='b';
end
s=x>=0;
k=s(2:end)-s(1:end-1);
if any(m=='p')
    f=find(k>0);
elseif any(m=='n')
    f=find(k<0);
else
    f=find(k~=0);
end
s=x(f+1)-x(f);
t=f-x(f)./s;
if any(m=='r')
    t=round(t);
end
if ~nargout
    n=length(x);
    plot(1:n,x,'-',t,zeros(length(t),1),'o');
end

repo: ekfmonoslam


text:
% All Pass Filter Sensitivities, RSS & WCA
% File: c:\M_files\short_updates\allpassrss.m
% Circuit function: G4.m
% updated 11/16/06
%
clear;clc
K=1e3;u=1e-6;
R1=626.25;R2=22.55*K;R3=25.05*K;R4=225.45*K;
C1=0.1*u;C2=C1;
Nom=[R1 R2 R3 R4 C1 C2]; % vector of nominal component values
BF=450;LF=550;NP=101; % linear freq sweep
F=linspace(BF,LF,NP);
% Form symmetric tolerance array T
Tr=0.01;Tc=0.05;
%
T=[-Tr -Tr -Tr -Tr -Tc -Tc;Tr Tr Tr Tr Tc Tc];
Nc=size(T,2); % Nc = number of components
% For assymetric tolerances if present
Mr=1+(T(2,:)+T(1,:))/2;
Tv=(T(2,:)-T(1,:))./(2*Mr); % Used only in RSS
Nav=Nom.*Mr; % Shift components to average value if
% tolerances are assymetric.  Mr = all 1's if symmetric.
%
[An,Bn,Dn,En,I]=G4(Nom); % Nominal SS arrays
%
dpf=0.0001; % derivative perturbation factor
rd=180/pi; % convert radians to degrees
%
% Q & R perturbation vectors are sequenced below:
% Reset; Q = [1 1 1 1 1]; R = [1 1 1 1 1]
% p = 1; Q = [1.0001 1 1 1 1];R = [0.9999 1 1 1 1];
% p = 2; Q = [1 1.0001 1 1 1];R = [1 0.9999 1 1 1];
% and so forth up to p = Nc.
%
Q=1+dpf;R=1-dpf;
%
for i=1:NP % Begin frequency sweep
   Qx=ones(1,Nc);Rx=ones(1,Nc); % Reset perturbation vectors
   s=2*pi*F(i)*j;
   Vo(i)=rd*angle(Dn*((s*I-An)\Bn)+En); % Nominal output
   for p=1:Nc % Begin component loop
      Qx(p)=Q;Rx(p)=R;
      if p > 1;Qx(p-1)=1;Rx(p-1)=1;end; % Reset previous
% Perturbate components forward with Q = 1+dpf
      [A,B,D,E]=G4(Nom.*Qx);
      Vr=rd*angle(D*((s*I-A)\B)+E);
% Perturbate components backward with R = 1-dpf
      [A,B,D,E]=G4(Nom.*Rx);
      Vb=rd*angle(D*((s*I-A)\B)+E);
      %
      Sen(i,p)=(Vr-Vb)/(2*Vo(i)*dpf); % Centered difference approximation
      %
      % For EVA
      %
      if Sen(i,p) > 0
         Lo(p)=1+T(1,p);Hi(p)=1+T(2,p);
      else
         Lo(p)=1+T(2,p);Hi(p)=1+T(1,p);
      end
   end % end component loop
   %
   % Get EVA VL and EVA VH
   %
   [A,B,D,E]=G4(Nom.*Lo);
   VL(i)=rd*angle(D*((s*I-A)\B)+E);
   [A,B,D,E]=G4(Nom.*Hi);
   VH(i)=rd*angle(D*((s*I-A)\B)+E);
   %
   % Get RSS using norm function
   %
   STn=norm(Sen(i,:).*Tv);
   Vrss1(i)=Vo(i)*(1-STn);Vrss2(i)=Vo(i)*(1+STn);
end % close frequency sweep loop i
%
subplot(2,2,1)
h=plot(F,Sen(:,5),'r',F,Sen(:,6),'b');
set(h,'LineWidth',2);
grid on
axis auto
set(gca,'FontSize',8)
%xlabel('Freq (Hz)');
ylabel('%/%')
title('C Sensitivities')
legend('C1','C2',0);
%
subplot(2,2,2)
h=plot(F,Sen(:,1),'k',F,Sen(:,2),'r',F,Sen(:,3),'b',F,Sen(:,4),'g');
set(h,'LineWidth',2);
grid on
axis auto
set(gca,'FontSize',8)
%xlabel('Freq (Hz)');
ylabel('%/%')
title('R Sensitivities')
legend('R1','R2','R3','R4',0);
%
subplot(2,2,4)
m=plot(F,VL,'b',F,VH,'r',F,Vo,'k--');
set(m,'LineWidth',2)
grid on
axis auto
set(gca,'FontSize',8)
%axis([BF LF 30 180])
%YT=linspace(30,180,6);
%set(gca,'ytick',YT);
xlabel('Freq (Hz)');
ylabel('Degrees')
title('EVA')
legend('EVLo','EVHi','Nom',0)
%
subplot(2,2,3)
m=plot(F,Vrss1,'b',F,Vrss2,'r',F,Vo,'k--');
set(m,'LineWidth',2)
grid on
axis auto
set(gca,'FontSize',8)
%axis([BF LF 30 180])
%YT=linspace(30,180,6);
%set(gca,'ytick',YT);
xlabel('Freq (Hz)');
ylabel('Degrees')
title('RSS')
legend('RSSLo','RSSHi','Nom',0)
%
figure(1);

author: Sable

repo: mcbench-benchmarks


text:
function varargout=readcnx(varargin)
%
% For calling details please see v_readcnx.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_readcnx.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: readcnx.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_readcnx(varargin{:});
else
    v_readcnx(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
classdef ImageFolder < datasets.Dataset
%IMAGEFOLDER Dataset created from a directory of JPEG images
%   Encapsulates a directory of JPEG images for training. The MatConvNet
%   function VL_IMREADJPEG is used, which supports fast multi-threaded
%   reading from disk.
%
%   D = datasets.ImageFolder('dataDir', '/folder', 'imageSize', [H, W])
%   loads a dataset from the directory '/folder', resizing images to have
%   height H and width W. Any subfolders are recursively scanned for JPEG
%   files. The resulting list of image files will be stored in D.filenames.
%
%   D.train() returns a cell array with mini-batches, from the shuffled
%   training set. Each mini-batch consists of a set of indexes drawn from
%   the vector D.trainSet.
%
%   D.val() returns a cell array with mini-batches, from the validation set
%   (without shuffling). Each mini-batch consists of a set of indexes drawn
%   from the vector D.valSet.
%
%   IMAGES = D.get(BATCH) returns a tensor of images for the given
%   mini-batch BATCH.
%
%   datasets.ImageFolder(..., 'option', value, ...) sets the following
%   properties, in addition to 'dataDir' and 'imageSize':
%
%   `batchSize`:: 128
%     The batch size.
%
%   `numThreads`:: 1
%     Number of image reading threads. If > 1, multi-threading is enabled.
%
%   `useGpu`:: false
%     Whether the IMAGES tensor is returned as a gpuArray.
%
%   `keepAspect`:: true
%     Whether aspect ratio is maintained when resizing images. If true,
%     images that do not conform to the aspect ratio of D.imageSize will be
%     cropped.
%
%   `removeMean`:: true
%     Whether the dataset's mean is subtracted from every pixel.
%
%   `whiten`:: true
%     Whether the dataset's standard deviation is removed from every pixel.
%
%   `augmentation`:: none
%     Struct that specifies data augmentation options. The valid fields
%     are: flip, location, aspect, scale, brightness, contrast, saturation,
%     crop. They correspond to similarly-named arguments of VL_IMREADJPEG,
%     with the exception of 'aspect' which corresponds to VL_IMREADJPEG's
%     'CropAnisotropy'. See 'help vl_imreadjpeg' for more details.
%
%   `partialBatches`:: false
%     Whether partial batches are returned (which can happen for the last
%     batch in a set, if the batch size does not divide the set size).
%
%   `trainSet`:: []
%     The training set. Note that subclasses (e.g. datasets.ImageNet)
%     usually override this property with the default training set.
%
%   `valSet`:: []
%     The validation set. Note that subclasses (e.g. datasets.ImageNet)
%     usually override this property with the default validation set.
%
%   `augmentImage`:: training only
%     Whether data augmentation is applied to each image in the dataset
%     (logical vector).
%
%   See 'autonn/examples/cnn/image_folder_example.m' for a full example.

% Copyright (C) 2018 Joao F. Henriques, 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).

  properties
    dataDir
    filenames

    numThreads = 1  % higher to enable multithreading/prefetching
    useGpu = false
    removeMean = true
    whiten = true

    imageSize
    keepAspect = true

    augmentation = struct('flip', false, 'location', false, 'aspect', 1, ...
      'scale', 1, 'brightness', 0, 'contrast', 0, 'saturation', 0, 'crop', 1)

    augmentImage  % vector, for each image, false disables augmentation (e.g. val images)

    % pixel statistics over dataset (computed automatically)
    rgbMean
    rgbCovariance
    rgbDeviation  % independent standard deviation of each channel (3x1)
    rgbDeviationFull  % standard deviation as a 3x3 matrix
    initialized = false
  end

  methods
    function o = ImageFolder(varargin)
      % parse generic ImageFolder arguments on non-subclassed
      % construction. allows use as a stand-alone class.
      varargin = o.parseGenericArgs(varargin) ;

      opts.skipInitialization = false ;
      opts = vl_argparse(opts, varargin) ;

      if ~opts.skipInitialization
        o.initialize() ;
      end
    end

    function initialize(o)
      % finishes initializing the dataset, must be called by subclasses
      assert(~o.initialized) ;
      assert(~isempty(o.dataDir), 'Must specify dataset.dataDir.') ;
      assert(~isempty(o.imageSize), 'Must specify dataset.imageSize.') ;
      o.initialized = true ;

      % list image files if they weren't set
      if isempty(o.filenames)
        o.filenames = o.listImages(o.dataDir) ;
      end

      % apply data augmentation to training images only by default
      if isempty(o.augmentImage)
        o.augmentImage = true(size(o.filenames)) ;
        o.augmentImage(o.valSet) = false ;
      end

      % load or recompute RGB statistics over all images
      if o.removeMean || o.whiten
        cache = [o.dataDir '/datasetstats.mat'] ;
        if exist(cache, 'file')
          s = load(cache) ;
          o.rgbMean = s.rgbMean ;
          o.rgbCovariance = s.rgbCovariance ;
        else
          o.computeImageStats() ;
          rgbMean = o.rgbMean ;  %#ok<*PROP,NASGU>
          rgbCovariance = o.rgbCovariance ;  %#ok<NASGU>
          save(cache, 'rgbMean', 'rgbCovariance') ;
        end

        % standard deviation
        o.rgbDeviation = sqrt(diag(o.rgbCovariance)) ;
        [v, d] = eig(o.rgbCovariance) ;
        o.rgbDeviationFull = v * sqrt(d) ;
      end

      % reset MEX file state
      clear('vl_imreadjpeg') ;
    end

    function files = listImages(o, folder, prefix)
      % returns all JPEG image files recursively from the given directory

      if nargin < 3  % no leading slash on top-level prefix/partial folder
        prefix = '' ;
      else
        prefix = [prefix '/'] ;
      end
      folder = [folder '/'] ;

      % list folder contents
      s = dir([folder '*']) ;

      % iterate files/folders
      files = cell(1, numel(s)) ;
      for i = 1:numel(s)
        name = s(i).name ;
        if strcmpi(name(max(1,end-3) : end), '.jpg') || ...
           strcmpi(name(max(1,end-4) : end), '.jpeg')
          % a file, store it (wrapped around a scalar cell)
          files{i} = {[prefix name]} ;

        elseif s(i).isdir && ~any(strcmp(name, {'.', '..'}))
          % a folder, recurse (stored as a cell array)
          files{i} = o.listImages([folder name], [prefix name]) ;
        end
      end

      % flatten nested cell arrays. this also removes empty entries, which
      % were not assigned in the For loop (e.g. invalid files/folders).
      files = [files{:}] ;
    end

    function args = parseGenericArgs(o, args)
      % called by subclasses to parse generic Dataset arguments.
      % start by parsing parent class's arguments
      args = o.parseGenericArgs@datasets.Dataset(args) ;

      % now parse ImageFolder arguments
      args = vl_parseprop(o, args, {'dataDir', 'filenames', 'removeMean', ...
        'whiten', 'numThreads', 'useGpu', 'imageSize', 'keepAspect', ...
        'augmentation', 'augmentImage'}) ;
    end

    function batches = partition(o, idx, batchSz)
      % partition indexes into batches (stored in a cell array).
      % if IDX is a matrix, each column is a distinct sample.
      if nargin < 3  % allow overriding batch size
        batchSz = o.batchSize ;
      end
      batches = o.partition@datasets.Dataset(idx, batchSz) ;  % call parent class

      % for prefetching, append an extra row with the next batch for each
      % batch. this way, get() can know which batch comes next to prefetch.
      if o.numThreads > 1
        batches = [batches; batches(1,2:end), {[]}] ;
      end
    end

    function [images, idx] = get(o, batch)
      % return a single batch of images, with optional prefetching.
      % also returns the image indexes.
      if ~iscell(batch)
        batch = {batch};
      end
      assert(~isempty(batch) && numel(batch) <= 2, ...
        'Malformed batch; use dataset.partition/train/val.');

      % get current batch
      idx = batch{1} ;
      images = o.getImageBatch(o.filenames(idx), o.augmentImage(idx(1))) ;

      % start prefetching next batch of images
      if o.numThreads > 1 && numel(batch) == 2 && ~isempty(batch{2})
        o.getImageBatch(o.filenames(batch{2}), o.augmentImage(batch{2}(1))) ;
      end
    end

    function data = getImageBatch(o, images, augmentImages)
      % load image batch, with data augmentation and whitening
      assert(o.initialized, 'Subclass did not call ImageFolder.initialize.') ;

      imagePaths = strcat([o.dataDir filesep], images) ;

      % for validation, disable all augmentation except any fixed cropping
      if augmentImages
        augment = o.augmentation ;  % full augmentation
      else
        augment = struct('crop', o.augmentation.crop) ;  % crop only
      end

      % fill in missing augmentation options with defaults
      defaults = struct('flip', false, 'location', false, 'aspect', 1, ...
        'scale', 1, 'brightness', 0, 'contrast', 0, 'saturation', 0, 'crop', 1) ;
      augment = vl_argparse(defaults, augment) ;

      % make brightness jitter proportional to dataset's color deviation,
      % unless the data is whitened
      if ~isequal(augment.brightness, 0) && ~o.whiten
        if ~isempty(o.rgbDeviationFull)
          augment.brightness = double(augment.brightness * o.rgbDeviationFull)  ;
        else  % no color deviation to use
          augment.brightness = 0 ;
        end
      end

      args{1} = {imagePaths, ...
                 'NumThreads', o.numThreads, ...
                 'Pack', ...
                 'Interpolation', 'bicubic', ...
                 'Resize', o.imageSize(1:2), ...
                 'CropSize', augment.crop * augment.scale, ...
                 'CropAnisotropy', augment.aspect, ...
                 'Brightness', augment.brightness, ...
                 'Contrast', augment.contrast, ...
                 'Saturation', augment.saturation} ;

      if ~o.keepAspect  % Squashing effect
        args{end+1} = {'CropAnisotropy', 0} ;
      end

      if augment.flip
        args{end+1} = {'Flip'} ;
      end

      if augment.location
        args{end+1} = {'CropLocation', 'random'} ;
      else
        args{end+1} = {'CropLocation', 'center'} ;
      end

      if o.useGpu
        args{end+1} = {'Gpu'} ;
      end

      if o.removeMean && ~isempty(o.rgbMean)
        args{end+1} = {'SubtractAverage', o.rgbMean} ;
      end

      args = horzcat(args{:}) ;

      if o.numThreads > 1 && nargout == 0
        vl_imreadjpeg(args{:}, 'prefetch') ;
      else
        data = vl_imreadjpeg(args{:}) ;
        data = data{1} ;

        if o.whiten && ~isempty(o.rgbDeviation)
          data = data ./ reshape(o.rgbDeviation, 1, 1, 3) ;
        end
      end
    end

    function computeImageStats(o)
      % computes RGB mean/covariance over whole dataset
      bs = o.batchSize ;
      images = o.filenames ;
      [avg, rgbm1, rgbm2] = deal({}) ;

      maxImages = 10000 ;  % max number of images, should be enough for pixel-wise mean
      skip = floor(max(1, numel(images) / maxImages)) ;  % number of batches to skip
      for t = 1 : skip * bs : numel(images)
        time = tic ;
        batch = t : min(t+bs-1, numel(images)) ;
        fprintf('collecting image stats: batch starting with image %d ...', batch(1)) ;

        data = o.getImageBatch(images(batch), false) ;

        z = reshape(shiftdim(data,2),3,[]) ;
        rgbm1{end+1} = mean(z,2) ;
        rgbm2{end+1} = z*z'/size(z,2) ;
        avg{end+1} = mean(data, 4) ;
        time = toc(time) ;
        fprintf(' %.1f Hz\n', numel(batch) / time) ;
      end

      % averageImage = gather(mean(cat(4,avg{:}),4)) ;
      rgbm1 = gather(mean(cat(2,rgbm1{:}),2)) ;
      rgbm2 = gather(mean(cat(3,rgbm2{:}),3)) ;
      o.rgbMean = rgbm1 ;
      o.rgbCovariance = rgbm2 - rgbm1*rgbm1' ;

      fprintf('computeImageStats: all done\n') ;
    end
  end
end

repo: MFT


text:
function planC = generate_DICOM_UID_Relationships(planC,structRefForC)
%"generate_DICOM_UID_Relationships"
%   Creates new fields in the passed planC containing the UIDs required for
%   DICOM export.  These fields are not a permenant part of the planC, so
%   this function should NOT be called on a global planC so that the new
%   fields are destroyed upon return.
%
%   As additional export capabilities are added, new UID fields may need to
%   be defined in this function.
%
%JRA 07/05/06
%
%Usage:
%   planC = generate_DICOM_UID_Relationships(planC)
%
% 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/>.

% Organization root for DOCOM UIDs
orgRoot = '1.3.6.1.4.1.9590.100.1.2';

% Initialize Series Number
seriesNumber = 7700;

% Read CERROptions.json
optS = getCERROptions();

indexS = planC{end};

%Generate a master study UID for all parts of this planC.
%Study_Instance_UID = dicomuid;
%if isfield(planC{indexS.scan}(1).scanInfo(1),'DICOMHeaders') && ...
%            ~isempty(planC{indexS.scan}(1).scanInfo(1).DICOMHeaders)
if isfield(planC{indexS.scan}(1).scanInfo(1),'studyInstanceUID') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).studyInstanceUID)
    Study_Instance_UID = planC{indexS.scan}(1).scanInfo(1).studyInstanceUID;
else
    %Study_Instance_UID = dicomuid;
    Study_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
end

if isfield(planC{indexS.scan}(1).scanInfo(1),'patientID') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).patientID)
    Patient_ID = planC{indexS.scan}(1).scanInfo(1).patientID;
else
    %Patient_ID = dicomuid;
    Patient_ID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
end

if isfield(planC{indexS.scan}(1).scanInfo(1),'studyDate') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).studyDate)
    Study_Date = planC{indexS.scan}(1).scanInfo(1).studyDate;
else
    Study_Date = datestr(now,'yyyymmdd');
end

if isfield(planC{indexS.scan}(1).scanInfo(1),'studyTime') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).studyTime)
    Study_Time = planC{indexS.scan}(1).scanInfo(1).studyTime;
else
    Study_Time = datestr(now,'hhmmss');
end

if isfield(planC{indexS.scan}(1).scanInfo(1),'patientName') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).patientName)
    patientName = planC{indexS.scan}(1).scanInfo(1).patientName;
else
    patientName = '';
end

if isfield(planC{indexS.scan}(1).scanInfo(1),'patientBirthDate') && ...
        ~isempty(planC{indexS.scan}(1).scanInfo(1).patientBirthDate)
    Patient_Birth_Date = planC{indexS.scan}(1).scanInfo(1).patientBirthDate;
else
    Patient_Birth_Date = '';
end

%% SCAN UIDs
%Iterate over CT scans
for i = 1:length(planC{indexS.scan})

    %Set the study instance UID.
    planC{indexS.scan}(i).Patient_ID = Patient_ID;

    %Set patient birth date.
    planC{indexS.scan}(i).Patient_Birth_Date = Patient_Birth_Date;

    %Set the study instance UID.
    planC{indexS.scan}(i).Study_Instance_UID = Study_Instance_UID;

    %Set the study date
    planC{indexS.scan}(i).Study_Date = Study_Date;

    %Set the study time
    planC{indexS.scan}(i).Study_Time = Study_Time;

    %Set the series number
    seriesNumber = seriesNumber + 1;
    planC{indexS.scan}(i).SeriesNumber = seriesNumber;

    %Generate a series instance UID for each scan;
    %planC{indexS.scan}(i).Series_Instance_UID = dicomuid;
    %if isfield(planC{indexS.scan}(i).scanInfo(1), 'DICOMHeaders') && ...
    %        ~isempty(planC{indexS.scan}(i).scanInfo(1).DICOMHeaders)
    if strcmpi(optS.retainOriginalUIDonExport,'yes') && ...
            isfield(planC{indexS.scan}(i).scanInfo(1),'seriesInstanceUID') && ...
            ~isempty(planC{indexS.scan}(i).scanInfo(1).seriesInstanceUID)
        planC{indexS.scan}(i).Series_Instance_UID = ...
            planC{indexS.scan}(i).scanInfo(1).seriesInstanceUID;
    else
        Series_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
        %planC{indexS.scan}(i).Series_Instance_UID = dicomuid;
        planC{indexS.scan}(i).Series_Instance_UID = Series_Instance_UID;
    end

    %Generate a frame of reference UID for each scan.
    %planC{indexS.scan}(i).Frame_Of_Reference_UID = dicomuid;
    %if isfield(planC{indexS.scan}(i).scanInfo(1), 'DICOMHeaders') && ...
    %        ~isempty(planC{indexS.scan}(i).scanInfo(1).DICOMHeaders)
    if strcmpi(optS.retainOriginalUIDonExport,'yes') && ...
            isfield(planC{indexS.scan}(i).scanInfo(1),'frameOfReferenceUID') && ...
            ~isempty(planC{indexS.scan}(i).scanInfo(1).frameOfReferenceUID)
        planC{indexS.scan}(i).Frame_Of_Reference_UID = ...
            planC{indexS.scan}(i).scanInfo(1).frameOfReferenceUID;
    else
        Frame_Of_Reference_UID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot));
        planC{indexS.scan}(i).Frame_Of_Reference_UID = Frame_Of_Reference_UID;
        %planC{indexS.scan}(i).Frame_Of_Reference_UID = dicomuid;
    end

    % %     %Generate a SOP ClassUID for each scan.
    % %     modality = planC{indexS.scan}(i).scanInfo(1).imageType;
    % %
    % %     switch lower(modality)
    % %         case {'ct scan'}
    % %             SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.2';
    % %         case {'mr scan'}
    % %             SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.4';
    % %         case {'pet scan'}
    % %             SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.129';
    % %         otherwise
    % %             %wy error('Currently unsupported scanType.  Only CT Scan is supported.');
    % %             SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.2';
    % %     end

    if ~isempty(planC{indexS.scan}(i).scanInfo(1).imageType)
        modality =planC{indexS.scan}(i).scanInfo(1).imageType;
    else
        modality = planC{indexS.scan}.scanInfo(1).imageType;
    end

    if any(strfind(upper(modality), 'CT'))
        SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.2';

    elseif any(strfind(upper(modality), 'MR'))
        SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.4';

    else
        SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.2';
    end

    %Iterate over slices
    if strcmpi(optS.retainOriginalUIDonExport,'yes') && ...
            isfield(planC{indexS.scan}(i).scanInfo(1),'sopInstanceUID') && ...
            ~isempty(planC{indexS.scan}(i).scanInfo(1).sopInstanceUID)
        for j=1:length(planC{indexS.scan}(i).scanInfo)
            sopInstanceUID = planC{indexS.scan}(i).scanInfo(j).sopInstanceUID;
            sopClassUID = planC{indexS.scan}(i).scanInfo(j).sopClassUID;
            planC{indexS.scan}(i).scanInfo(j).SOP_Class_UID      = sopClassUID;
            planC{indexS.scan}(i).scanInfo(j).SOP_Instance_UID   = sopInstanceUID;
            if isempty(planC{indexS.scan}(i).scanInfo(j).imageOrientationPatient)
                planC{indexS.scan}(i).scanInfo(j).imageOrientationPatient = [1,0,0,0,1,0];
            end
        end
    else
        for j=1:length(planC{indexS.scan}(i).scanInfo)
            planC{indexS.scan}(i).scanInfo(j).SOP_Class_UID      = SOP_Class_UID;
            Frame_Of_Reference_UID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot));
            planC{indexS.scan}(i).scanInfo(j).SOP_Instance_UID   = Frame_Of_Reference_UID;
            %planC{indexS.scan}(i).scanInfo(j).SOP_Instance_UID   = dicomuid;
            if isempty(planC{indexS.scan}(i).scanInfo(j).imageOrientationPatient)
                planC{indexS.scan}(i).scanInfo(j).imageOrientationPatient = [1,0,0,0,1,0];
            end
        end
    end

end

%% DOSE UIDs
%Iterate over doses
for i = 1:length(planC{indexS.dose})

    %Set the patient id.
    planC{indexS.dose}(i).Patient_ID = Patient_ID;

    %Set patient birth date.
    planC{indexS.dose}(i).Patient_Birth_Date = Patient_Birth_Date;

    % Set the patient name
    planC{indexS.dose}(i).patientName = patientName;

    %Set the study instance UID.
    planC{indexS.dose}(i).Study_Instance_UID = Study_Instance_UID;

    %Set the series number
    seriesNumber = seriesNumber + 1;
    planC{indexS.dose}(i).SeriesNumber = seriesNumber;

    planC{indexS.dose}(i).Study_Date = Study_Date;

    planC{indexS.dose}(i).Study_Time = Study_Time;

    %Generate a series instance UID for each dose;
    %planC{indexS.dose}(i).Series_Instance_UID = dicomuid;
    Series_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
    planC{indexS.dose}(i).Series_Instance_UID = Series_Instance_UID;

    %Set the frame of reference UID to that of the associated scan, if it
    %exists.
    aS = getDoseAssociatedScan(i, planC);
    if ~isempty(aS)
        planC{indexS.dose}(i).Frame_Of_Reference_UID = planC{indexS.scan}(aS).Frame_Of_Reference_UID;
        planC{indexS.dose}(i).imageOrientationPatient = planC{indexS.scan}(aS).scanInfo(1).imageOrientationPatient;
    else
        %planC{indexS.dose}(i).Frame_Of_Reference_UID = dicomuid;
        Frame_Of_Reference_UID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot));
        planC{indexS.dose}(i).Frame_Of_Reference_UID = Frame_Of_Reference_UID;
        warning('No associated scan found. Assuming HFS orientation.');
        planC{indexS.dose}(i).imageOrientationPatient = [1 0 0 0 1 0];
    end

    %Generate a SOP ClassUID for each dose.
    SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.481.2';

    planC{indexS.dose}(i).SOP_Class_UID = SOP_Class_UID;
    %planC{indexS.dose}(i).SOP_Instance_UID = dicomuid;
    SOP_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
    planC{indexS.dose}(i).SOP_Instance_UID = SOP_Instance_UID;

    %Generate a Referenced RT Plan Sequence SOP Class UID for each dose.
    %When plan support is added this must reference actual plans, but for
    %now they are dummy plans.
    planC{indexS.dose}(i).Referenced_RT_Plan_Sequence_SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.481.5';
    %planC{indexS.dose}(i).Referenced_RT_Plan_Sequence_SOP_Instance_UID = dicomuid;
    Referenced_RT_Plan_Sequence_SOP_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
    planC{indexS.dose}(i).Referenced_RT_Plan_Sequence_SOP_Instance_UID = Referenced_RT_Plan_Sequence_SOP_Instance_UID;

end

%% STRUCTURE UIDs

%Generate a single series instance UID for all structures.
%Structure_Set_Series_UID = dicomuid;
Structure_Set_Series_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);

%Generate a single SOP_Class and SOP_Instance ID for all structures;
Structure_Set_SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.481.3';
%Structure_Set_SOP_Instance_UID = dicomuid;
Structure_Set_SOP_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);

%Iterate over structures
seriesNumber = seriesNumber + 1;
for i = 1:length(planC{indexS.structures})

    %Set the patient id.
    planC{indexS.structures}(i).Patient_ID = Patient_ID;

    %Set patient birth date.
    planC{indexS.structures}(i).Patient_Birth_Date = Patient_Birth_Date;

    % Set the patient name
    planC{indexS.structures}(i).patientName = patientName;

    %Set the study instance UID.
    planC{indexS.structures}(i).Study_Instance_UID = Study_Instance_UID;

    %Set the series number
    planC{indexS.structures}(i).SeriesNumber = seriesNumber;

    planC{indexS.structures}(i).Study_Date = Study_Date;

    planC{indexS.structures}(i).Study_Time = Study_Time;

    %Set the series instance UID
    planC{indexS.structures}(i).Series_Instance_UID = Structure_Set_Series_UID;

    %Set the default frame of reference UID to that of the associated scan.
    assocScanNum = getStructureAssociatedScan(i, planC);

    % Set referenced frame of reference UID
    structRefFrameOfReferenceUID = planC{indexS.scan}(assocScanNum).Frame_Of_Reference_UID;

    % Set referenced series UID
    refSeriesInstanceUID = planC{indexS.scan}(assocScanNum).Series_Instance_UID;
    %planC{indexS.scan}(assocScanNum).scanInfo(1).seriesInstanceUID;

    % Set sereferced sopInstanceUIDs
    sopClassUidc = {planC{indexS.scan}(assocScanNum).scanInfo(:).sopClassUID};
    sopInstanceUidc = {planC{indexS.scan}(assocScanNum).scanInfo(:).sopInstanceUID};

%     % Handle special case of assignig reference UID fromanother structure
%     % (e.g. exporting registered images from MIM assistant changes their frameOfreferenceUID)
    if exist('structRefForC','var') && ~isempty(structRefForC)
        ind = ismember(planC{indexS.structures}(i).structureName,structRefForC(:,1));
        if sum(ind) == 1
            structRefFrameOfReferenceUID = structRefForC{ind,2};
            refSeriesInstanceUID = structRefForC{ind,3};
            sopClassUidc = structRefForC{ind,4};
            sopInstanceUidc = structRefForC{ind,5};
        elseif any(ismember('all',structRefForC(:,1)))
                structRefFrameOfReferenceUID = structRefForC{1,2};
                refSeriesInstanceUID = structRefForC{1,3};
                sopClassUidc = structRefForC{1,4};
                sopInstanceUidc = structRefForC{1,5};
        end
    end

    planC{indexS.structures}(i).Frame_Of_Reference_UID = structRefFrameOfReferenceUID;
    planC{indexS.structures}(i).Referenced_Series_Instance_UID = refSeriesInstanceUID;

    %Set the SOP UIDs.
    planC{indexS.structures}(i).SOP_Class_UID    = Structure_Set_SOP_Class_UID;
    planC{indexS.structures}(i).SOP_Instance_UID = Structure_Set_SOP_Instance_UID;

    % set the referenced SOP UID for each slice
    for slc = 1:length(planC{indexS.structures}(i).contour)
        planC{indexS.structures}(i).contour(slc).SOP_Class_UID = sopClassUidc{slc};
        planC{indexS.structures}(i).contour(slc).SOP_Instance_UID = sopInstanceUidc{slc};
    end

    aS = getStructureAssociatedScan(i, planC);
    if ~isempty(aS)
        planC{indexS.structures}(i).imageOrientationPatient = planC{indexS.scan}(aS).scanInfo(1).imageOrientationPatient;
    else
        warning('No associated scan found. Assuming HFS orientation.');
        planC{indexS.structures}(i).imageOrientationPatient = [1 0 0 0 1 0];
    end

end

%% DVH UIDs

%Set the associated structure set UIDs for all DVHs to the same values.
Referenced_Structure_Set_SOP_Class_UID    = Structure_Set_SOP_Class_UID;
Referenced_Structure_Set_SOP_Instance_UID = Structure_Set_SOP_Instance_UID;

%Iterate over DVHs
seriesNumber = seriesNumber + 1;
for i = 1:length(planC{indexS.DVH})
    planC{indexS.DVH}(i).Referenced_Structure_Set_SOP_Class_UID     = Referenced_Structure_Set_SOP_Class_UID;
    planC{indexS.DVH}(i).Referenced_Structure_Set_SOP_Instance_UID  = Referenced_Structure_Set_SOP_Instance_UID;
    planC{indexS.DVH}(i).SeriesNumber = seriesNumber;
end

%% GSPS UIDs
%Iterate over GSPS
for scanNum = 1:length(planC{indexS.scan})
    originalSliceSOPInstanceUIDc = {planC{indexS.scan}(scanNum).scanInfo.sopInstanceUID};
    seriesNumber = seriesNumber + 1;
    for i = 1:length(planC{indexS.GSPS})

        % scanNum associated with this GSPS (expand to handle multiple scans)
        %scanNum = 1;

        %Set the patient id
        planC{indexS.GSPS}(i).Patient_ID = Patient_ID;

        %Set the patient birth date
        planC{indexS.GSPS}(i).Patient_Birth_Date = Patient_Birth_Date;

        % Set the patient name
        planC{indexS.GSPS}(i).patientName = patientName;

        %Set the study instance UID.
        planC{indexS.GSPS}(i).Study_Instance_UID = Study_Instance_UID;

        planC{indexS.GSPS}(i).SeriesNumber = seriesNumber;

        planC{indexS.GSPS}(i).Study_Date = Study_Date;

        planC{indexS.GSPS}(i).Study_Time = Study_Time;

        %Generate a series instance UID for each dose;
        %planC{indexS.GSPS}(i).Series_Instance_UID = dicomuid;
        %Series_Instance_UID = javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot);
        %planC{indexS.GSPS}(i).Series_Instance_UID = Series_Instance_UID;

        %Set the frame of reference UID to that of the associated scan, if it
        %exists.
        planC{indexS.GSPS}(i).Frame_Of_Reference_UID = planC{indexS.scan}(1).Frame_Of_Reference_UID;

        %Generate a SOP ClassUID for each GSPS.
        SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.11.1';

        planC{indexS.GSPS}(i).SOP_Class_UID = SOP_Class_UID;
        %planC{indexS.GSPS}(i).SOP_Instance_UID = dicomuid;
        SOP_Instance_UID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot));
        planC{indexS.GSPS}(i).SOP_Instance_UID = SOP_Instance_UID;

        if strcmpi(optS.retainOriginalUIDonExport,'no')
            sopInstanceUID = planC{indexS.GSPS}(i).SOPInstanceUID;
            sliceNum = strncmp(sopInstanceUID,originalSliceSOPInstanceUIDc,length(sopInstanceUID));
            referenced_SOP_instance_uid = ...
                planC{indexS.scan}(scanNum).scanInfo(sliceNum).SOP_Instance_UID;
            referenced_SOP_class_uid = ...
                planC{indexS.scan}(scanNum).scanInfo(sliceNum).SOP_Class_UID;
            Series_Instance_UID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',orgRoot));
            planC{indexS.GSPS}(i).Series_Instance_UID = Series_Instance_UID;
        else
            referenced_SOP_instance_uid = planC{indexS.GSPS}(i).SOPInstanceUID;
            referenced_SOP_class_uid = ...
                planC{indexS.scan}(scanNum).scanInfo(1).SOP_Class_UID;
            series_Instance_UID = planC{indexS.scan}(scanNum).scanInfo(1).seriesInstanceUID;
            planC{indexS.GSPS}(i).Series_Instance_UID = series_Instance_UID;
        end
        planC{indexS.GSPS}(i).referenced_SOP_instance_uid = referenced_SOP_instance_uid;
        planC{indexS.GSPS}(i).referenced_SOP_class_uid = referenced_SOP_class_uid;

        %Generate a Referenced Study Sequence SOP Class UID for each dose.
        %planC{indexS.GSPS}(i).Referenced_Study_Sequence_SOP_Class_UID = '1.2.840.10008.5.1.4.1.1.481.5';
        %planC{indexS.GSPS}(i).Referenced_Study_Sequence_SOP_Instance_UID = dicomuid;

        %Generate a Referenced Series Sequence SOP Class UID for each dose.

    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ calmars ] = calmar( cumReturns, mdds, frequency )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Doyen Sahoo
% Contributors: Steven Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [r c] = size(cumReturns);
    calmars = zeros(1,c);
    den = 252/frequency;
    Y = r/den;
    annualisedReturns = ((cumReturns(r,:)./cumReturns(1,:)) .^ (1/Y)) - 1;
    calmars = annualisedReturns ./ mdds;
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function stroud_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests BALL_F1_ND, BALL_F3_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_ND_INDEX;

  n_max = 3;
  num = function_nd_num ( );

  xc = [ 1.0, -1.0, 2.0 ];
  r = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  For integrals in a ball in ND:\n' );
  fprintf ( 1, '  BALL_F1_ND approximates the integral;\n' );
  fprintf ( 1, '  BALL_F3_ND approximates the integral.\n' );
  fprintf ( 1, '\n' );

  for n = 2 : n_max

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Spatial dimension N = %d\n', n );
    fprintf ( 1, '  Ball center:\n' );
    for i = 1 : n
      fprintf ( 1, '%12f  ', xc(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Ball radius = %f\n', r );
    fprintf ( 1, '  Ball volume = %f\n', ball_volume_nd ( n, r ) );
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '    Rule:      F1          F3\n' );
    fprintf ( 1, '    F(X)\n' );
    fprintf ( 1, '\n' );

    for i = 1 : num

      FUNC_ND_INDEX = i;

      result1 = ball_f1_nd ( 'function_nd', n, xc, r );
      result2 = ball_f3_nd ( 'function_nd', n, xc, r );

      fname = function_nd_name ( i );

      fprintf ( 1, '  %7s  %12f  %12f\n', fname, result1, result2 );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Example
% Controlling KUKA iiwa from V-rep simulation using MATLAB and KST as a
% middle-ware.
% This script works together with V-rep scene: iiwaFromVrep.ttt

%% To run this example do the following:
% 1- Run the Vrep simulation: iiwaFromVrep.ttt
% 2- Run the (MarlabToolboxServer) application on the robot
% 3- Run this script from MATLAB.

% Copyright: Mohammad SAFEEA, 19th-April-2018

%% Initiation part of the code
disp('Program started');
kst_Path=getTheKSTDirectory(pwd);
addpath(kst_Path);
vrep=remApi('remoteApi');
vrep.simxFinish(-1);
clientID=vrep.simxStart('127.0.0.1',19999,true,true,5000,5);
ip='172.31.1.147';
global t_Kuka;
if (clientID>-1)
    t_Kuka=net_establishConnection( ip );
    if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
      warning('Connection could not be establised, script aborted');
      return;
    end
    jHandles=zeros(7,1);
    for i=1:7
            s=['LBR_iiwa_7_R800_joint',num2str(i)];
            [res, daHandle]=vrep.simxGetObjectHandle(clientID,s,vrep.simx_opmode_oneshot_wait);
            jHandles(i)=daHandle;
    end
else
    return;
end
jPos={0., pi / 180 * 20., 0., -pi / 180 * 70., 0.,...
                        pi / 180 * 90., 0.};
relVel=0.25;
movePTPJointSpace( t_Kuka , jPos, relVel);
%% Control loop
realTime_startDirectServoJoints( t_Kuka );
pause(10);
counter=0;
tic;
totalTimeSecs=60; % 60 seconds of execution time
vrep.simxStartSimulation(clientID,vrep.simx_opmode_oneshot_wait);
for i=1:7
        [res,tempPos]=vrep.simxGetJointPosition(clientID,jHandles(i),...
            vrep.simx_opmode_streaming);
end
vrep.simxWriteStringStream(clientID,'path_velocity','0.07',vrep.simx_opmode_oneshot_wait);
jpos=zeros(7,25000);
t_0=toc;
while(t_0<totalTimeSecs)
    for i=1:7
            [res,tempPos]=vrep.simxGetJointPosition(clientID,jHandles(i),...
                vrep.simx_opmode_buffer);
            jPos{i}=tempPos;
    end

    if(toc-t_0>0.003)
        counter=counter+1;
        sendJointsPositionsf( t_Kuka ,jPos);
        for i=1:7
            jpos(i,counter)=jPos{i};
        end
        t_0=toc;
	end
end
%% Turn off KST server
net_turnOffServer( t_Kuka );
%% Stop simulation
vrep.simxStopSimulation(clientID,vrep.simx_opmode_oneshot_wait);
%% End simulation
vrep.simxFinish(clientID);
vrep.delete();

repo: KST-Kuka-Sunrise-Toolbox


text:
% Question No: 4

% Detect the line segments in a binary image using Hough transform

function houghtr(x)
im=imread(x);
f=edge(im,'canny');
figure,imshow(f),title('Image after applying Canny Filter');
[h,theta,rho]=hough1(f,0.5);
imshow(theta,rho,h,[],'notruesize'),axis on,axis normal;
xlabel('\theta'),ylabel('\rho');
[r,c]=houghpeaks1(h,5);
hold on;
plot(theta(c),rho(r),'linestyle','none','marker','s','color','w');
lines=houghlines1(f,theta,rho,r,c);
figure,imshow(f),hold on;
for k=1:length(lines)
    xy=[lines(k).point1;lines(k).point2];
    plot(xy(:,2),xy(:,1),'Linewidth',4,'Color',[.6 .6 .6]);
end
end

function [h,theta,rho]=hough1(f,dtheta,drho)
if nargin<3
    drho=1;
end
if nargin<2
    dtheta=1;
end
f=double(f);
[m,n]=size(f);
theta=linspace(-90,0,ceil(90/dtheta)+1);
theta=[theta -fliplr(theta(2:end-1))];
ntheta=length(theta);
d=sqrt((m-1)^2+(n-1)^2);
q=ceil(d/drho);
nrho=2*q-1;
rho=linspace(-q*drho,q*drho,nrho);
[x,y,val]=find(f);
x=x-1;y=y-1;
h=zeros(nrho,length(theta));
for k=1:ceil(length(val)/1000)
    first=(k-1)*1000+1;
    last=min(first+999,length(x));
    x_matrix=repmat(x(first:last),1,ntheta);
    y_matrix=repmat(y(first:last),1,ntheta);
    val_matrix=repmat(val(first:last),1,ntheta);
    theta_matrix=repmat(theta,size(x_matrix,1),1)*pi/180;
    rho_matrix=x_matrix.*cos(theta_matrix)+y_matrix.*sin(theta_matrix);
    slope=(nrho-1)/(rho(end)-rho(1));
    rho_bin_index=round(slope*(rho_matrix-rho(1))+1);
    theta_bin_index=repmat(1:ntheta,size(x_matrix,1),1);
    h=h+full(sparse(rho_bin_index(:),theta_bin_index(:),val_matrix(:),nrho,ntheta));
end
end

function [r,c,hnew]=houghpeaks1(h,numpeaks,threshold,nhood)
if nargin<4
    nhood=size(h)/50;
    nhood=max(2*ceil(nhood/2)+1,1);
end
if nargin<3
    threshold=0.5*max(h(:));
end
if nargin<2
    numpeaks=1;
end
done=false;
hnew=h;r=[];c=[];
while ~done
    [p,q]=find(hnew==max(hnew(:)));
    p=p(1);q=q(1);
    if hnew(p,q)>=threshold
        r(end+1)=p;
        c(end+1)=q;
        p1=p-(nhood(1)-1)/2;
        p2=p+(nhood(1)-1)/2;
        q1=q-(nhood(2)-1)/2;
        q2=q+(nhood(2)-1)/2;
        [pp,qq]=ndgrid(p1:p2,q1:q2);
        pp=pp(:);qq=qq(:);
        badrho=find((pp<1)|(pp>size(h,1)));
        pp(badrho)=[];
        qq(badrho)=[];
        theta_too_low=find(qq<1);
        qq(theta_too_low)=size(h,2)+qq(theta_too_low);
        pp(theta_too_low)=size(h,1)-pp(theta_too_low)+1;
        theta_too_high=find(qq>size(h,2));
        qq(theta_too_high)=qq(theta_too_high)-size(h,2);
        pp(theta_too_high)=size(h,1)-pp(theta_too_high)+1;
        hnew(sub2ind(size(hnew),pp,qq))=0;
        done=length(r)==numpeaks;
    else
        done=true;
    end
end
end

function [r,c]=houghpixels1(f,theta,rho,rbin,cbin)
[x,y,val]=find(f);
x=x-1;
y=y-1;
theta_c=theta(cbin)*pi/180;
rho_xy=x*cos(theta_c)+y*sin(theta_c);
nrho=length(rho);
slope=(nrho-1)/(rho(end)-rho(1));
rho_bin_index=round(slope*(rho_xy-rho(1))+1);
idx=find(rho_bin_index==rbin);
r=x(idx)+1;
c=y(idx)+1;
end

function lines=houghlines1(f,theta,rho,rr,cc,fillgap,minlength)
if nargin<6
    fillgap=20;
end
if nargin<7
    minlength=40;
end
numlines=0;
lines=struct;
for k=1:length(rr)
    rbin=rr(k);
    cbin=cc(k);
    [r,c]=houghpixels1(f,theta,rho,rbin,cbin);
    if isempty(r)
        continue
    end
    omega=(90-theta(cbin))*pi/180;
    t=[cos(omega) sin(omega);-sin(omega) cos(omega)];
    xy=[r-1 c-1]*t;
    x=sort(xy(:,1));
    diff_x=[diff(x);Inf];
    idx=[0;find(diff_x>fillgap)];
    for p=1:length(idx)-1
        x1=x(idx(p)+1);
        x2=x(idx(p+1));
        linelength=x2-x1;
        if linelength>=minlength
            point1=[x1 rho(rbin)];
            point2=[x2 rho(rbin)];
            tinv=inv(t);
            point1=point1*tinv;
            point2=point2*tinv;
            numlines=numlines+1;
            lines(numlines).point1=point1+1;
            lines(numlines).point2=point2+1;
            lines(numlines).length=linelength;
            lines(numlines).theta=theta(cbin);
            lines(numlines).rho=rho(rbin);
        end
    end
end
end


text:
function [T, Y]=FOChuaM(parameters, orders, TSim, Y0)
%
% Numerical Solution of the Fractional-Order Chua's System
% with piecewise linear nonlinearity defined by a memristor.
%
%   D^q1 x(t) = alpha(y(t) - x(t) + epsilon x(t) - W(w)x(t))
%   D^q2 y(t) = x(t) - y(t) + z(t)
%   D^q3 z(t) = -beta y(t) - gamma z(t)
%   D^q4 w(t) = x(t)
%   where  W(x) is defined as W_w(a, b, w);
%
% function [T, Y] = FOChuaM(parameters, orders, TSim, Y0)
%
% Input:    parameters - model parameters [alpha, beta, gamma, epsilon, a, b]
%           orders - derivatives orders [q1, q2, q3, q4]
%           TSim - simulation time (0 - TSim) in sec
%           Y0 - initial conditions [Y0(1), Y0(2), Y0(3), Y0(4)]
%
% Output:   T - simulation time (0 : Tstep : TSim)
%           Y - solution of the system (x=Y(1), y=Y(2), z=Y(3), w=Y(4))
%
% Author:  (c) Ivo Petras (ivo.petras@tuke.sk), 2010.
%

% time step:
h=0.005;
% number of calculated mesh points:
n=round(TSim/h);
%orders of derivatives, respectively:
q1=orders(1); q2=orders(2); q3=orders(3); q4=orders(4);
% constants of Chua's system:
alpha=parameters(1); beta=parameters(2); gamma=parameters(3);
epsilon=parameters(4); a=parameters(5); b=parameters(6);
% binomial coefficients calculation:
cp1=1; cp2=1; cp3=1; cp4=1;
for j=1:n
    c1(j)=(1-(1+q1)/j)*cp1;
    c2(j)=(1-(1+q2)/j)*cp2;
    c3(j)=(1-(1+q3)/j)*cp3;
    c4(j)=(1-(1+q4)/j)*cp4;
    cp1=c1(j); cp2=c2(j);
    cp3=c3(j); cp4=c4(j);
end
% initial conditions setting:
x(1)=Y0(1); y(1)=Y0(2); z(1)=Y0(3); w(1)=Y0(4);
% calculation of phase portraits /numerical solution/:
for i=2:n
    x(i)=(alpha*(y(i-1)-x(i-1)+epsilon*x(i-1)-W_w(a,b,w(i-1))*x(i-1)))*h^q1 - memo(x, c1, i);
    y(i)=(x(i)-y(i-1)+z(i-1))*h^q2 - memo(y, c2, i);
    z(i)=(-beta*y(i)-gamma*z(i-1))*h^q3 - memo(z, c3, i);
    w(i)=(x(i))*h^q4 - memo(w, c4, i);
end
for j=1:n
    Y(j,1)=x(j);
    Y(j,2)=y(j);
    Y(j,3)=z(j);
    Y(j,4)=w(j);
end
T=h:h:TSim;
%


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

for i=1:size(X, 1),
    for j=1:K,
        x = X(i, :)';
        projection_k = x' * U(:, j);
        Z(i, j) = projection_k;
    end
end

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

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 8 Assignments/K-Means Clustering and PCA/mlclass-ex7/projectData.m


text:
function K = levelsetCurvature(phi,mode)
%levelsetCurvature Computes the curvature of a level set function.
%   K = levelsetCurvature(PHI,MODE) computes the curvature, K, of level
%   set function PHI. If MODE = 'Cu' the curvature is not normalized
%   (this is the default). If MODE = 'Cn', the curvature is normalized
%   by dividing it by its maximum value. This function is based on Eqs.
%   (12-61) and (12-62) of DIPUM3E.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% PRELIMINARIES
% Set default.
if nargin == 1
    mode = 'Cu';
end
[M,N] = size(phi);
phi = double(phi);

% PAD WITH 1'S TO BE ABLE TO COMPUTE THE DERIVATIVE AT IMAGE BORDERS.
phip = padarray(phi,[1 1],1,'both');

% COMPUTE THE CENTRAL DIFFERENCES.
phix = 0.5*(phip(3:end,2:N+1) - phip(1:M,2:N+1));
phiy = 0.5*(phip(2:M+1,3:end) - phip(2:M+1,1:N));
phixx = phip(3:end,2:N+1) + phip(1:M,2:N+1) - 2*phi;
phiyy = phip(2:M+1,3:end) + phip(2:M+1,1:N) - 2*phi;
phixy = 0.25.*(phip(3:end,3:end) - phip(1:M,3:end)...
                            - phip(3:end,1:N) + phip(1:M,1:N));
% COMPUTE THE CURVATURE.
DEN = (phix.^2 + phiy.^2 + eps).^1.5;
K = (phixx.*(phiy.^2) - 2*(phix.*phiy.*phixy) + phiyy.*(phix.^2))./DEN;

% CHECK FOR NORMALIZATION.
if isequal(mode,'Cn')
    K = K/max(abs(K(:)));
end

% COMPENSATE FOR BORDER EFFECTS BY SETTING THE FIRST AND LAST ROWS AND
% COLUMNS EQUAL TO THEIR IMMEDIATE PREDECESSORS.
% CURVATURE
% Rows.
K(:,1) = K(:,2);
K(:,end) = K(:,end-1);
% Columns.
K(1,:) = K(2,:);
K(end,:) = K(end-1,:);

author: dipum

repo: dipum-toolbox

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


text:
classdef MimProjectsSidePanel < GemListBoxWithTitle
    % MimProjectsSidePanel. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimProjectsSidePanel is part of the side panel and contains the sliding list
    %     box showing project names
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        PatientDatabase
        GuiCallback
    end

    properties (SetAccess = private)
        CurrentProjectId
    end

    methods
        function obj = MimProjectsSidePanel(parent, patient_database, gui_callback)
            obj = obj@GemListBoxWithTitle(parent, 'PROJECT', 'Import images', 'Delete project');
            obj.PatientDatabase = patient_database;
            obj.GuiCallback = gui_callback;
        end

        function RepopulateSidePanel(obj, project_id)
            obj.AddProjectsToListBox(project_id);
        end

        function project_has_changed = UpdateSidePanel(obj, project_id)
            project_has_changed = ~strcmp(project_id, obj.CurrentProjectId);
            if project_has_changed
                obj.ListBox.SelectItem(project_id, true);
            end
            obj.CurrentProjectId = project_id;
        end

        function SelectProject(obj, project_id, selected)
            obj.ListBox.SelectItem(project_id, selected);
        end

    end

    methods (Access = protected)

        function AddButtonClicked(obj, ~, event_data)
        end

        function DeleteButtonClicked(obj, ~, event_data)
        end
    end

    methods (Access = private)
        function AddProjectsToListBox(obj, selected_project_id)
            [project_names, project_ids] = obj.PatientDatabase.GetListOfProjects;
            obj.ListBox.ClearItems;

            for index = 1 : numel(project_ids)
                project_name = project_names{index};
                project_id = project_ids{index};
                patient_item = MimProjectNameListItem(obj.ListBox.GetListBox, project_name, project_id, obj.GuiCallback);
                obj.ListBox.AddItem(patient_item);
            end

            obj.ListBox.SelectItem(selected_project_id, true);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test_tutorial_multivariateanalysis(datadir, dmltdir)

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockstatistics ft_topoplotER ft_freqstatistics ft_topoplotTFR

% this is a test script that I made following the report of Matt on
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1585 and realizing that the
% wiki page http://www.fieldtriptoolbox.org/tutorial/multivariateanalysis was
% not included in a test script.

% This test script corresponds to the documentation on the wiki from 3 July 2012

if nargin==0
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/classification');
  dmltdir = dccnpath('/home/common/matlab/fieldtrip/external/dmlt');
end

addpath(genpath(dmltdir));

filename = dccnpath(fullfile(datadir, 'covatt'));
load(filename);

cfg             = [];
cfg.parameter   = 'trial';
cfg.keeptrials  = 'yes'; % classifiers operate on individual trials
cfg.channel     = {'MLO' 'MRO'}; % occipital channels only

tleft   = ft_timelockanalysis(cfg,left);
tright  = ft_timelockanalysis(cfg,right);

cfg         = [];
cfg.layout  = 'CTF275.lay';
cfg.method  = 'crossvalidate';

cfg.design  = [ones(size(tleft.trial,1),1); 2*ones(size(tright.trial,1),1)]';
cfg.latency = [2.0 2.5]; % final bit of the attention period
stat = ft_timelockstatistics(cfg,tleft,tright);

disp(stat.statistic)

% this is the part that Matt mentioned in bug 1585
cfg.statistic = {'accuracy' 'binomial' 'contingency'};

stat = ft_timelockstatistics(cfg,tleft,tright);

disp(stat.statistic.contingency)

stat.mymodel = stat.model{1}.primal;

cfg              = [];
cfg.parameter    = 'mymodel';
cfg.layout       = 'CTF275.lay';
cfg.xlim         = [2.0 2.5];
cfg.comments     = '';
cfg.colorbar     = 'yes';
cfg.interplimits = 'electrodes';
ft_topoplotER(cfg,stat);

cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 8:2:14;
cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;
cfg.channel      = {'MLO' 'MRO'};
cfg.toi          = 2.0:0.1:2.5;
cfg.keeptrials   = 'yes'; % classifiers operate on individual trials

tfrleft          = ft_freqanalysis(cfg, left);
tfrright         = ft_freqanalysis(cfg, right);

cfg         = [];
cfg.layout  = 'CTF275.lay';
cfg.method  = 'crossvalidate';
cfg.design  = [ones(size(tfrleft.powspctrm,1),1); 2*ones(size(tfrright.powspctrm,1),1)]';
stat        = ft_freqstatistics(cfg,tfrleft,tfrright);

disp(stat.statistic)

stat.mymodel = stat.model{1}.primal;

cfg              = [];
cfg.layout       = 'CTF275.lay';
cfg.parameter    = 'mymodel';
cfg.comment      = '';
cfg.colorbar     = 'yes';
cfg.interplimits = 'electrodes';
ft_topoplotTFR(cfg,stat);

cfg         = [];
cfg.layout  = 'CTF275.lay';
cfg.method  = 'crossvalidate';
cfg.design  = [ones(size(tfrleft.powspctrm,1),1); 2*ones(size(tfrright.powspctrm,1),1)]';
cfg.mva     = {dml.standardizer dml.enet('family','binomial','alpha',0.2)};

stat = ft_freqstatistics(cfg,tfrleft,tfrright);

disp(stat.statistic)

stat.mymodel     = stat.model{1}.weights;

cfg              = [];
cfg.layout       = 'CTF275.lay';
cfg.parameter    = 'mymodel';
cfg.comment      = '';
cfg.colorbar     = 'yes';
cfg.interplimits = 'electrodes';
ft_topoplotTFR(cfg,stat);


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             MI
%   - pre-registration     rigid2D
%   - regularizer          mbCurvature
%   - optimization         lBFGS
% ===============================================================================

close all, help(mfilename);

setup2DPETCTData
imgModel('reset','imgModel','splineInter','regularizer','none','theta',1e-3);
distance('reset','distance','MI','nT',32,'nR',32);
trafo('reset','trafo','rigid2D');
regularizer('reset','regularizer','mfCurvature','alpha',1e-1);

PIRpara = optPara('lBFGS','solver','backslash');
NPIRpara = optPara('lBFGS','solver',regularizer('get','solver'),'maxIter',40);

[yc,wc,his] = MLIR(ML,'PIRobj',@PIRBFGSobjFctn,'PIRpara',PIRpara,...
  'NPIRobj',@NPIRBFGSobjFctn,'NPIRpara',NPIRpara,...
  'minLevel',4,'maxLevel',7,'parametric',1,'plotMLiter',0);

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


text:
function p=v_sone2phon(s)
%V_PHON2SONE convert SONE loudness values to PHONs p=(s)
%Inputs:    s is a matrix of sone values
%
%Outputs:   p is a matrix, the same size as s, of phon values
%
% The phon scale measures perceived loudness in dB; at 1 kHz it is identical to dB SPL
% relative to 20e-6 Pa sound pressure. The sone scale is proportional to apparent loudness
% and, by definition, equals 1 at 40 phon. The form of the loudness curve is taken from [1].
% The hearing threshold at 1 kHz for 18 to 25 year olds with normal hearing is taken from [2].
%
% Refs: [1]	J. Lochner and J. Burger. Form of the loudness function in the presence of masking noise.
%           The Journal of the Acoustical Society of America, 33: 1705, 1961.
%       [2]	ISO/TC43. Acoustics, Normal equal-loudness-level contours.
%           Standard ISO 226:2003, Aug. 2003.

%      Copyright (C) Mike Brookes 2012-2013
%      Version: $Id: v_sone2phon.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b d
if isempty(a)
    b=1/(log(10)*0.1*0.27); % 0.27 is the exponent from [1] and [2]
    d=exp(2.4/b); % 2.4 dB is teh hearing threshold from [2]
    a=exp(40/b)-d; % scale factor to make p=40 give s=1
end
if nargout>0

    p=b*log(a*s+d);
else
    if nargin<1 || isempty(s)
        pp=linspace(5,90,100)'; % phon values
        ss=v_phon2sone(pp);
    else
        ss=s;
    end
    semilogx(ss,v_sone2phon(ss));
    v_axisenlarge(-1);
    v_xticksi;
    ylabel('phon = dB SPL @ 1 kHz');
    xlabel('sone');
end


text:
% SCRIPT TEST THE DIRECT DYNAMICS OF A 3 DOF PLANAR ROBOT ROBOT

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

fprintf('\nTHE SIMULATION PRESENTS THE ROBOT AT AN INITIAL POSITION WHEN NO TORQUES ARE APPLIED\n')

%load robot parameters
robot=load_robot('example', '3dofplanar');

%simulate for 50 seconds, change this depending on your computer speed and
%the total time that you want to simulate
total_simulation_time = 10;

%initial position and joint speed
q0 = [0 0 0]';
qd0 = [0 0 0]';

drawrobot3d(robot, q0);
adjust_view(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The student should try different combinations of the following parameters:
%   g: the direction of the gravity vector. In this case, if we select g=[0  0 9.81]';
%       the movement of the arm is not affected by the gravity, since it
%       moves in a plane perpendicular to the gravity vector.
%   tau: the torques applied to each joint. The student should observe the effects of selecting
%       tau = [0 0 0]' or different values in combination with the
%       direction of the vector g.
%   robot.dynamics.friction = 0 selects no friction at the joints, whereas
%       robot.dynamics.friction = 1 considers that there exists friction. This
%       friction is modelled by robot.motors.Viscous (viscous friction) and
%       robot.motors.Coulomb (Coulomb friction). The student should observe
%       that selecting g=[0  9.81 0]' and tau = [0 0 0]' and
%       robot.dynamics.friction = 0 turns into an infinite triple pendulum
%       movement. In addition, selecting selecting g=[0  9.81 0]' and tau = [0 0 0]' and
%       robot.dynamics.friction = 1 simulates the case in which the triple
%       pendulum converges to a steady solution with the three links
%       hanging along the Y direction.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%you may redefine the gravity vector
%in this case you may one of the next two lines, that define
%the gravity acting along the Y axis or the Z axis, respectively.
g=[0  -9.81 0]'; %y axis
%g=[0  0 9.81]'; % Z axis

%tau = [3 2 1]';
tau = [0 0 0]';%no torques applied

%select friction or not
robot.dynamics.friction = 0;

fprintf('\nCOMPUTING FORWARD DYNAMICS (this may take a while)')

%this may take a while, since it requires integration
%of the acceleration at each time step
[t, q, qd] = forwarddynamic(robot, total_simulation_time, q0, qd0, tau, g, []);

%animate it!!
animate(robot, q)

figure, plot(t, q), grid, title('Position vs. time')
xlabel('time (s)'), ylabel('Position (rad)')
legend('q_1', 'q_2', 'q_3', 'q_4', 'q_5', 'q_6');

figure, plot(t, qd), grid, title('Speed vs. time')
xlabel('time (s)'), ylabel('Speed (rad/s)')
legend('qd_1', 'qd_2', 'qd_3', 'qd_4', 'qd_5', 'qd_6');

repo: ARTE


text:
function p=v_normcdflog(x,m,s)
%V_NORMCDFLOG calculates log of Normal Cumulative Distribution function p=(x,m,s)
%
% Inputs:
%
%     X        Input data (vector or matrix)
%     M        Mean of Normal distribution [default 0]
%     S        Std deviation of Normal distribution [default 1]
%
% Outputs:
%
%     P        P = log(normcdf(X)); same size as X
%
% The routine gives accurate values even if X is large and negative

%      Copyright (C) Mike Brookes 2016
%      Version: $Id: v_normcdflog.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b
if isempty(a)
    a=0.996;
%   f=@(x) -0.5*(x.^2+log(2*pi))-real(log(-(a+x.^2)./x))-real(log(normcdf(x)));
%   b=fzero(f,-22);         % cutoff value for conventional formula
    b=-22.2491306156561;    % precalculated value
end
if nargin>2
    x=(x-m)/s;
elseif nargin>1
    x=x-m;
end
t=x<b;                                                  % mask for large negative values
p=zeros(size(x));
p(~t)=real(log(0.5*erfc(-x(~t).*sqrt(0.5))));           % use this formula normally
p(t)=-0.5*(x(t).^2+log(2*pi))-real(log(-x(t)-a./x(t))); % use this approximation for large negative x


text:
function mf = meanfilt1(data,n)

for i=1:length(data)
    mf(i)=mean(data(max(1,i-round((n-1)/2)):min(i+round((n-1)/2),length(data))));
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
Network VGG16 {

	Layer CONV1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C:3, R: 3, S: 3, Y: 224, X: 224 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV2 {
		Type: CONV
		Dimensions { K: 64, C: 64, R: 3, S: 3, Y: 224, X: 224 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}

	}

	Layer CONV3 {
		Type: CONV
		Dimensions { K: 128, C: 64, R: 3, S: 3, Y: 112, X: 112 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV4 {
		Type: CONV
		Dimensions { K: 128, C: 128, R: 3, S: 3, Y: 112, X: 112 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV5 {
		Type: CONV
		Dimensions { K: 256, C:128, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV6 {
		Type: CONV
		Dimensions { K: 256, C: 256, Y: 56, X: 56, R: 3, S: 3 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV7 {
		Type: CONV
		Dimensions { K: 256, C: 256, R: 3, S: 3, Y: 56, X: 56 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV8 {
		Type: CONV
		Dimensions { K: 512, C: 256, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV9 {
		Type: CONV
		Dimensions { K: 512, C:512, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV10 {
		Type: CONV
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV11 {
		Type: CONV
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV12 {
		Type: CONV
		Dimensions { K: 512, C: 512, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}

	Layer CONV13 {
		Type: CONV
		Dimensions { K: 512, C: 512, Y: 14, X: 14, R: 3, S: 3 }
		Dataflow {
			// Only one spatial map can be applied except on R and S
			TemporalMap(1,1) C;      // Flexible (Mapping size/offset, spatial/temporal)
			SpatialMap(1,1) K;       // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) Y';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(1,1) X';     // Flexible (Mapping size/offset, spatial/temporal)
			TemporalMap(Sz(R),Sz(R)) R; // This cannot be changed
			TemporalMap(Sz(S),Sz(S)) S; // This cannot be changed
			// Virtual Neuron of size Sz(R) x Sz(S):
			Cluster(Sz(R),P);
			SpatialMap(1,1) Y;
			SpatialMap(1,1) R;
			Cluster(Sz(S),P);
			SpatialMap(1,1) X;
			SpatialMap(1,1) S;
		}
	}
}

author: maestro-project

repo: maestro


text:
function varargout = fillSphericalTriangle(sphere, p1, p2, p3, varargin)
%FILLSPHERICALTRIANGLE Fill a triangle on a sphere.
%
%   fillSphericalTriangle(SPHERE, PT1, PT2, PT3);
%
%
%   See also
%   fillSphericalPolygon, drawSphericalTriangle, drawSphere
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 22/02/2005
%

%   HISTORY
%   27/06/2007 manage spheres other than origin
%   30/10/2008 replace intersectPlaneLine by intersectLinePlane
%   2012-02-09 rename as fillSphericalTriangle

% extract data of the sphere
ori = sphere(:, 1:3);
r   = sphere(4);

% extract direction vectors for each point
v1  = normalizeVector3d(p1 - ori);
v2  = normalizeVector3d(p2 - ori);
v3  = normalizeVector3d(p3 - ori);

% create a plane tangent to the sphere containing first point
plane = createPlane(v1, v1);

% position on the plane of the direction vectors
pp2 = planePosition(intersectLinePlane([ori v2], plane), plane);
pp3 = planePosition(intersectLinePlane([ori v3], plane), plane);

% create rough parametrization with 2 variables
nTri = 5;
s  = linspace(0, 1, nTri);
t  = linspace(0, 1, nTri);
ns = length(s);
nt = length(t);
s  = repmat(s, [nt, 1]);
t  = repmat(t', [1, ns]);

% convert to plane coordinates
xp = s * pp2(1) + t .* (1-s) * pp3(1);
yp = s * pp2(2) + t .* (1-s) * pp3(2);

% convert to 3D coordinates (still on the 3D plane)
x  = plane(1) * ones(size(xp)) + plane(4) * xp + plane(7) * yp - ori(1);
y  = plane(2) * ones(size(xp)) + plane(5) * xp + plane(8) * yp - ori(2);
z  = plane(3) * ones(size(xp)) + plane(6) * xp + plane(9) * yp - ori(3);

% project on the sphere
norm = hypot(hypot(x, y), z);
xn = x ./ norm * r + ori(1);
yn = y ./ norm * r + ori(2);
zn = z ./ norm * r + ori(3);

if nargout == 0
    % simply display the patch
    surf(xn, yn, zn, 'FaceColor', 'g', 'EdgeColor', 'none', varargin{:});

elseif nargout == 1
    % display the patch and return a handle
    h = surf(xn, yn, zn, 'FaceColor', 'g', 'EdgeColor', 'none', varargin{:});
    varargout = {h};

elseif nargout == 3
    % If 3 outputs are required, return patch vertex coordinates
    varargout = {x, y, z};
end

repo: ksptot


text:
% Test file for @chebfun/all.m.

function pass = test_all(pref)

if (nargin < 1)
    pref = chebfunpref();
end

% Test scalar valued chebfuns.
f = chebfun(@(x) sin(x), [-1 -0.5 0 0.5 1], pref);
pass(1) = ~all(f);

f = chebfun(@(x) sin(x - 0.1), [-1 -0.5 0 0.5 1], pref);
pass(2) = ~all(f);

f = chebfun(@(x) exp(2*pi*1i*x), [-1 -0.5 0 0.5 1], pref);
pass(3) = all(f);

% Test array-valued chebfun.
f = chebfun(@(x) [sin(x) sin(x - 0.1) exp(2*pi*1i*x)], [-1 -0.5 0 0.5 1], pref);
pass(4) = isequal(all(f), logical([0 0 1]));

% Test on SINGFUN:
f = chebfun(@(x) sin(x)./(x+1), 'exps', [-1 0]);
pass(5) = ~all(f);

% Test for function defined on unbounded domain:

% Blowing-up functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];

op = @(x) x.^2.*(1-exp(-x.^2))+3;
f = chebfun(op, dom, 'exps', [2 2]);
pass(6) = all(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun