text:
function [hfig,finalim,final_pforefly] = OverlayFrames_Ellipses(trx,readframe,bkgdim0,predictions,mainfly,otherflies,ts,varargin)

colorpos = [.7,0,0];
colorneg = [0,0,.7];
border = 20;

[colorpos,colorneg,border,hfig,figpos,...
  cm,fg_thresh,bg_thresh,sigma_bkgd,wmah,frac_a_back,dist_epsilon,...
  ncolors_reference,solidbkgdcolor] = ...
  myparse(varargin,'colorpos',colorpos,...
  'colorneg',colorneg,...
  'border',border,...
  'hfig',[],...
  'figpos',[],...
  'colormap',{@jet,@jet},...
  'fg_thresh',115/255,...
  'bg_thresh',10/255,...
  'sigma_bkgd',115/255-10/255,...
  'wmah',.5,...
  'frac_a_back',1,...
  'dist_epsilon',.1,...
  'ncolors_reference',[],...
  'solidbkgdcolor',[.85,.85,.85]);

if isempty(otherflies),
  flygroups = {mainfly};
else
  flygroups = {mainfly,otherflies};
end
nflygroups = numel(flygroups);
nframes = numel(ts);
allts = ts(1):ts(end);
nframes_total = numel(allts);
allflies = [flygroups{:}];

[nr0,nc0,~] = size(bkgdim0);

%% grab out relevant time points
% not rotating currently, maybe rotate in the future?

x = nan(trx.nflies,nframes_total);
y = nan(trx.nflies,nframes_total);
a = nan(trx.nflies,nframes_total);
b = nan(trx.nflies,nframes_total);
theta = nan(trx.nflies,nframes_total);
imcenter = [(nc0+1)/2,(nr0+1)/2];
for fly = 1:trx.nflies,
  idx = allts+trx(fly).off;
  i0 = find(idx>1,1);
  i1 = find(idx<trx(fly).nframes,1,'last');
  x(fly,i0:i1) = trx(fly).x(idx(i0:i1));
  y(fly,i0:i1) = trx(fly).y(idx(i0:i1));
  a(fly,i0:i1) = trx(fly).a(idx(i0:i1));
  b(fly,i0:i1) = trx(fly).b(idx(i0:i1));
  theta(fly,i0:i1) = trx(fly).theta(idx(i0:i1));
end

bkgdim = bkgdim0;

%% figure out the axes limits
xlim = [inf,-inf];
ylim = [inf,-inf];
for fly = allflies,
  for j = 1:nframes_total,
    [x1,x2,y1,y2] = ellipse_to_bounding_box(x(fly,j),y(fly,j),a(fly,j),b(fly,j),theta(fly,j));
    xlim(1) = min([xlim(1),x1,x2]);
    xlim(2) = max([xlim(2),x1,x2]);
    ylim(1) = min([ylim(1),y1,y2]);
    ylim(2) = max([ylim(2),y1,y2]);
  end
end
xlim = [floor(xlim(1))-border,ceil(xlim(2))+border];
ylim = [floor(ylim(1))-border,ceil(ylim(2))+border];
bkgdim = bkgdim(ylim(1):ylim(2),xlim(1):xlim(2));
x = x - xlim(1)+1;
y = y - ylim(1)+1;

%% colors for each frame
colors_all = nan(ts(end)-ts(1)+1,3,nflygroups);
for i = 1:nflygroups,
  cmcurr = cm{min(numel(cm),i)};
  if isempty(cmcurr),
    colors_all(:,:,i) = 1;
  elseif isnumeric(cmcurr),
    colors_all(:,:,i) = cmcurr;
  elseif ~isempty(ncolors_reference),
    colors_reference = cmcurr(ncolors_reference);
    idxcurr = round(linspace(1,ncolors_reference,ts(end)-ts(1)+1));
    colors_all(:,:,i) = colors_reference(idxcurr,:);
  else
    colors_all(:,:,i) = cmcurr(ts(end)-ts(1)+1);
  end
end
colors = colors_all(ts-ts(1)+1,:,:);

%% read in the frames, convert to grayscale
imcurr = readframe(ts(1));
if size(imcurr,3) > 1,
  imcurr = rgb2gray(imcurr);
end
%imcurr = imtransform(imcurr,tform,'XData',[1,nc0],'YData',[1,nr0],'XYScale',1);
imcurr = imcurr(ylim(1):ylim(2),xlim(1):xlim(2));
gray_ims = im2double(imcurr);

%% probability that each pixel is foreground
dbkgd = abs(bsxfun(@minus,gray_ims,bkgdim));
pback = exp(-dbkgd/sigma_bkgd^2/2);
pback = pback / max(pback(:));
minv = exp(-fg_thresh/sigma_bkgd^2/2);
maxv = exp(-bg_thresh/sigma_bkgd^2/2);
pback = min(1,max(0,(pback-minv) / (maxv-minv)));
pfore = 1-pback;

%% make a connectivity graph for computing distances within foreground
% pixels
nr = ylim(2)-ylim(1)+1;
nc = xlim(2)-xlim(1)+1;
np = nr*nc;

[cgrid,rgrid] = meshgrid(1:nc,1:nr);
rconn1 = repmat(rgrid,[1,1,4]);
cconn1 = repmat(cgrid,[1,1,4]);
rconn2 = nan(size(rconn1));
cconn2 = nan(size(cconn1));

rconn2(:,:,1) = rgrid+1; % above
cconn2(:,:,1) = cgrid;
rconn2(:,:,2) = rgrid-1; % below
cconn2(:,:,2) = cgrid;
rconn2(:,:,3) = rgrid; % right
cconn2(:,:,3) = cgrid+1;
rconn2(:,:,4) = rgrid; % left
cconn2(:,:,4) = cgrid-1;
badidx = rconn2 < 1 | rconn2 > nr | cconn2 < 1 | cconn2 > nc;
rconn1(badidx) = [];
cconn1(badidx) = [];
rconn2(badidx) = [];
cconn2(badidx) = [];
conn1 = sub2ind([nr,nc],rconn1,cconn1);
conn2 = sub2ind([nr,nc],rconn2,cconn2);

%isdy = conn2 == conn1+1 | conn2 == conn1-1;
%w = min(maxdist,-log(pfore)+dist_epsilon)/maxdist;
w = reshape((pback+dist_epsilon)/(1+dist_epsilon),[np,1]);

%% assign pixels to flies

% probability that each pixel belongs to each fly
pfly = zeros(nr*nc,trx.nflies,1);
%allflies = [flygroups{:}];

ii = 1;
t = ts(ii);
i = t-ts(1)+1;
for fly = 1:trx.nflies,
  if t > trx(fly).endframe || t < trx(fly).firstframe,
    continue;
  end
  mu = [x(fly,i)-frac_a_back*a(fly,i)*cos(theta(fly,i)),...
    y(fly,i)-frac_a_back*a(fly,i)*sin(theta(fly,i))];
  %     mu = [trx(fly).x(j)-frac_a_back*trx(fly).a(j)*cos(trx(fly).theta(j)),...
  %       trx(fly).y(j)-frac_a_back*trx(fly).a(j)*sin(trx(fly).theta(j))];

  s = sub2ind([nr,nc],min(nr,max(1,round(mu(2)))),...
    min(nc,max(1,round(mu(1)))));
  wcurr = w(conn2,ii);
  S = axes2cov(a(fly,i),b(fly,i),theta(fly,i));
  diffs = bsxfun(@minus,[cgrid(:),rgrid(:)],mu);
  c = chol(S);
  dmah = sum((diffs/c).^2',1); %#ok<UDIM>
  %     ratio = S(1,1)/S(2,2);
  %     rationorm = 2 / (ratio+1);
  %     wcurr(isdy) = wcurr(isdy)*ratio*rationorm;
  %     wcurr(~isdy) = wcurr(~isdy)*rationorm;
  G = sparse(conn1,conn2,wcurr,np,np);
  d = graphshortestpath(G,s,'Directed',true);
  %pfly(:,fly,i) = exp(-d.^2/(trx(fly).a(j)*dist_epsilon)^2);
  pfly(:,fly,ii) = exp(-(1-wmah)*d.^2/(a(fly,i)*dist_epsilon)^2 - wmah*dmah);
end

Z = sum(pfly,2);
Z(Z==0) = 1;
pfly = bsxfun(@rdivide,pfly,Z);

gray_ims = reshape(gray_ims,[nr*nc,1]);
pfore = reshape(pfore,[nr*nc,1]);
psomefly = reshape(1-prod(1-pfly(:,[flygroups{:}],:),2),[nr*nc,1]);
pforefly = pfore.*psomefly;
%pfore_any = 1 - prod(1-pfore.*psomefly,2);

%% color the images

im = zeros(nr*nc,3,1);
i = 1;
imcurr = zeros(nr*nc,3);
for flygroupi = 1:nflygroups,
  flygroup = flygroups{flygroupi};
  pcurr = 1 - prod(1-pfly(:,flygroup,i),2);
  colorcurr = colors(i,:,flygroupi);
  imcurr = imcurr + bsxfun(@times,colorcurr,pcurr);
end
if ~isempty(solidbkgdcolor),
  imcurr = bsxfun(@plus,gray_ims(:,i).*pforefly(:,i),...
    bsxfun(@times,1-pforefly(:,i),repmat(reshape(solidbkgdcolor,[1,3]),[nr*nc,1])));
else
  imcurr = bsxfun(@plus,gray_ims(:,i).*pforefly(:,i),...
    (1-pforefly(:,i)).*reshape(bkgdim,[nr*nc,1]));
end
im(:,:,i) = imcurr;
finalim = reshape(im,[nr,nc,3]);

final_pforefly = reshape(pforefly,[nr,nc]);

%% plot

if isempty(hfig),
  hfig = figure;
end

figure(hfig);
clf;
if ~isempty(figpos),
  set(hfig,'Units','pixels','Position',figpos);
end
hax = gca;

idxpos = predictions{mainfly}(ts(1):ts(end));
timestamps = trx(mainfly).timestamps(ts+trx(mainfly).off)-...
  trx(mainfly).timestamps(ts(1)+trx(mainfly).off);

i = 1;
t = ts(i);
image(finalim(:,:,:,i),'Parent',hax(i));
axis(hax(i),'image','off');
hold(hax(i),'on');
plot(hax(i),x(mainfly,:),y(mainfly,:),'k.-');
plot(hax(i),x(mainfly,idxpos),y(mainfly,idxpos),'.','Color',colorpos);
plot(hax(i),x(mainfly,~idxpos),y(mainfly,~idxpos),'.','Color',colorneg);
if predictions{mainfly}(t),
  colorcurr = colorpos;
else
  colorcurr = colorneg;
end
plot(hax(i),x(mainfly,t-ts(1)+1),y(mainfly,t-ts(1)+1),'o','color',colorcurr,'markerfacecolor',colorcurr);
text(1,1,sprintf('t = %.2fs',timestamps(i)),'HorizontalAlignment','left','VerticalAlignment','top','Parent',hax(i));

hell = nan(1,nframes);
for i = 1:nframes,
  j = ts(i)-ts(1)+1;
  hell(i) = drawellipse(x(mainfly,j),y(mainfly,j),theta(mainfly,j),2*a(mainfly,j),2*b(mainfly,j));
  set(hell(i),'LineWidth',2,'Color',colors(i,:,1));
end

colormap(colors_all(:,:,1));
tlim = timestamps(end);
set(hax,'CLim',[0,tlim]);
hcb = colorbar('East');
tticks = [0,tlim];
realylim = get(hcb,'YLim');
yticks = realylim(1)+(tticks/tlim)*(realylim(2)-realylim(1));
set(hcb,'YTick',yticks,'YTickLabel',tticks');

if isempty(figpos),
  truesize;
else
  set(hfig,'Units','pixels','Position',figpos);
end

author: kristinbranson

repo: JAABA


text:
function Power(nlp, sys)

    domains = sys.Gamma.Nodes.Domain;

    for i=1:numel(domains)
        domain = domains{i};

        u = domain.Inputs.Control.u;
        dq = domain.States.dx;
        B = domain.Gmap.Control.u;

        P = tovector(norm(dq.*(B*u)).^2);
        power_fun = SymFunction(['power_' sys.Gamma.Nodes.Name{i}],P,{u,dq});
        rs_phase = getPhaseIndex(nlp,sys.Gamma.Nodes.Name{i});
        addRunningCost(nlp.Phase(rs_phase),power_fun,{'u','dx'});

    end
    nlp.update;
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function nlev = numlevels(im_sz)
% number of pyramid levels for an image of size [r c]

% as many pyramid levels as possible, up to 1x1
min_d = min(im_sz(1:2));
nlev = 1;
while min_d>1
    nlev = nlev+1;
    min_d = floor((min_d+1)/2);
end

%nlev = floor(log(min(im_sz(1:2))) / log(2));
%nlev = ceil(log(min(im_sz(1:2))) / log(2));

end

author: drakeguan


text:
% Get number of space dimensions ("space dimension" field) according to the
% space descriptor, following the format definition at
% http://teem.sourceforge.net/nrrd/format.html#spacedirections.
%
% Returns 3 if "space" field is one of
% "right-anterior-superior" or "RAS"
% "left-anterior-superior" or "LAS"
% "left-posterior-superior" or "LPS"
% "scanner-xyz"
% "3D-right-handed"
% "3D-left-handed"
%
% Returns 4 if "space" field is one of
% "right-anterior-superior-time" or "RAST"
% "left-anterior-superior-time" or "LAST"
% "left-posterior-superior-time" or "LPST"
% "scanner-xyz-time"
% "3D-right-handed-time"
% "3D-left-handed-time"
%
% Date: October 25, 2017
% Author: Gaetan Rensonnet
function sd = nrrd_getSpaceDimensions(spacedescriptor)

if any(strcmpi(spacedescriptor,...
        {'right-anterior-superior', 'RAS',...
        'left-anterior-superior', 'LAS',...
        'left-posterior-superior', 'LPS',...
        'scanner-xyz',...
        '3D-right-handed',...
        '3D-left-handed'}...
        ))
    sd = 3;

elseif any(strcmpi(spacedescriptor,...
        {'right-anterior-superior-time', 'RAST',...
        'left-anterior-superior-time', 'LAST',...
        'left-posterior-superior-time', 'LPST', ...
        'scanner-xyz-time', ...
        '3D-right-handed-time', ...
        '3D-left-handed-time'}...
        ))
    sd = 4;

else
    sd = -1;
    % Unrecognized nrrd space descriptor (grace under fire)
end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%
%  函数imfilter的应用
f = imread('Fig0216(a).tif');
imshow(f);

w = ones(31);  %  生成滤波模板
gd = imfilter(f,w);
figure,imshow(gd,[ ])

gr = imfilter(f,w,'replicate');
figure,imshow(gr,[ ])

gs = imfilter(f,w,'symmetric');
figure,imshow(gs,[ ])

gc = imfilter(f,w,'circular');
figure,imshow(gc,[ ])

f8 = im2uint8(f);
g8r = imfilter(f,w,'replicate');
figure,imshow(g8r,[ ])

author: chengchengzi

repo: Digital-Image-Processing


text:
function sF = times(sF1,sF2)
% overloads sF1 .* sF2
%
% Syntax
%   sF = sF1.*sF2
%   sF = a.*sF1
%   sF = sF1.*a
%
% Input
%   sF1, sF2 - S2FunHarmonic
%
% Output
%   sF - S2FunHarmonic
%

sF = times@S2FunHarmonic(sF1,sF2);

% try to preserve symmetry
if isa(sF1,'S2FunHarmonicSym') && isa(sF2,'S2FunHarmonicSym')
  sym = disjoint(sF1.CS,sF2.CS);
  sF = S2FunHarmonicSym(sF.fhat,sym);
end

author: mtex-toolbox

repo: mtex

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


text:
function test_suite = test_surficial_io_niml_dset
% tests for AFNI NIML input/output
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_niml_dset_dataset_io()
    if cosmo_skip_test_if_no_external('afni')
        return;
    end

    ds=get_niml_dataset();
    fn=cosmo_make_temp_filename('_tmp','.niml.dset');
    cleaner=onCleanup(@()delete(fn));

    formats={'ascii','binary','binary.lsbfirst','binary.msbfirst'};

    n_formats=numel(formats);

    for k=1:n_formats
        format=formats{k};
        niml_str=get_niml_dset_string(format);

        fid=fopen(fn,'w');
        fwrite(fid,niml_str);
        fclose(fid);
        ds2=cosmo_surface_dataset(fn);
        assert_dataset_equal(ds,ds2);

        cosmo_map2surface(ds,fn,'encoding',format);
        ds3=cosmo_surface_dataset(fn);
        assert_dataset_equal(ds,ds3);

        % try mapping using functionality in AFNI Matlab library
        s=struct();
        s.data=ds.samples';
        s.node_indices=ds.a.fdim.values{1}(ds.fa.node_indices)-1;
        s.stats=ds.sa.stats;
        s.labels=ds.sa.labels;

        afni_niml_writesimple(s, fn, format);
        ds4=cosmo_surface_dataset(fn);
        assert_dataset_equal(ds,ds4);

        % ensure that mapping to struct works
        s2=cosmo_map2surface(ds,'-niml_dset','encoding',format);
        assertEqual(s2,s);

    end

function assert_dataset_equal(x,y)
    assertElementsAlmostEqual(x.samples,y.samples,'absolute',1e-4);
    assertEqual(sort(fieldnames(x)),sort(fieldnames(y)));
    assertEqual(x.fa,y.fa);
    assertEqual(x.sa,y.sa);
    assertEqual(x.a,y.a);

function ds=get_niml_dataset()
    ds=cosmo_synthetic_dataset('ntargets',2,'nchunks',1,...
                        'type','surface');
    ds.a.fdim.values{1}=1+[1 4 7 8 5 3];
    ds.sa=struct();
    ds.sa.stats={'Ftest(3,4)';'Zscore()'};
    % TODO: use these complicated labels
    % ds.sa.labels={'label1<>"&''';'label2'};
    ds.sa.labels={'label1';'label2'};

function string=get_niml_dset_string(format)
    string=[...
            get_niml_dset_header() ...
            get_niml_dset_data(format) ...
            get_niml_dset_footer()
            ];

function string=get_niml_dset_header()
    string=bytes_printf([...
            '<AFNI_dataset\n' ...
            '  dset_type="Node_Bucket"\n' ...
            '  self_idcode="XYZ_QJCWHYMMSUIUXCMBUZJASSJT"\n' ...
            '  filename="data.niml.dset"\n' ...
            '  label="data.niml.dset"\n' ...
            '  ni_form="ni_group" >\n' ...
            ]);

function string=get_niml_dset_footer()
    string=bytes_printf([...
            '<AFNI_atr\n' ...
            '  atr_name="COLMS_RANGE"\n' ...
            '  ni_type="String"\n' ...
            '  ni_dimen="1" >\n' ...
            '"-3.6849477 2.0316862 1 0;-1.3264908 2.3386572 2 4"'...
                        '</AFNI_atr>\n' ...
            '<AFNI_atr\n' ...
            '  atr_name="COLMS_LABS"\n' ...
            '  ni_type="String"\n' ...
            '  ni_dimen="1" >\n' ...
            '"label1;label2"</AFNI_atr>\n' ...
            '<AFNI_atr\n' ...
            '  atr_name="COLMS_TYPE"\n' ...
            '  ni_type="String"\n' ...
            '  ni_dimen="1" >\n' ...
            '"Generic_Float;Generic_Float"</AFNI_atr>\n' ...
            '<AFNI_atr\n' ...
            '  atr_name="COLMS_STATSYM"\n' ...
            '  ni_type="String"\n' ...
            '  ni_dimen="1" >\n' ...
            '"Ftest(3,4);Zscore()"</AFNI_atr>\n' ...
            '</AFNI_dataset>\n' ...
            ]);
    % TODO: use these complicated labels:
    % '"label1<>"&';label2"</AFNI_atr>\n' ...

function bytes=get_niml_dset_data(format)

    ds=get_niml_dataset();
    data=ds.samples;
    nodes_base0=ds.a.fdim.values{1}(ds.fa.node_indices)-1;

    [nrows,ncols]=size(data);

    row_sep=sprintf('\n');
    col_sep=' ';

    switch lower(format)
        case 'ascii'
            data_pat=[repmat(['%.6f' col_sep],1,nrows) row_sep];
            bucket_data=uint8(sprintf(data_pat,data));

            node_pat=['%d' row_sep];
            node_data=uint8(sprintf(node_pat,nodes_base0));
            ni_form='';

        case {'binary','binary.lsbfirst','binary.msbfirst'};
            data_single=single(data(:)');
            node_int=int32(nodes_base0(:)');

            [unused,unused,endian]=computer();
            computer_format=sprintf('binary.%ssbfirst',lower(endian));

            do_swap=~strcmp(format,'binary') && ...
                            ~strcmp(format, computer_format);

            if do_swap
                if cosmo_wtf('is_octave')
                    % Octave 3.8 does not support byteswap for single
                    data_single_int=typecast(data_single,'int32');
                    data_single_int_swap=swapbytes(data_single_int);
                    data_single=typecast(data_single_int_swap,'single');
                else
                    data_single=swapbytes(data_single);
                end
                node_int=swapbytes(node_int);
            end

            if strcmp(format,'binary')
                ni_form_format=computer_format;
            else
                ni_form_format=format;
            end

            bucket_data=typecast(data_single,'uint8');
            node_data=typecast(node_int,'uint8');

            ni_form=sprintf('  ni_form="%s"\n',ni_form_format);

        otherwise
            error('unsupported format %s', format);
    end

    bytes=[...
            bytes_printf([...
                '<SPARSE_DATA\n' ...
                '  data_type="Node_Bucket_data"\n' ...
                '  ni_type="2*float"\n' ...
                ni_form ...
                '  ni_dimen="6" >' ])...
            bucket_data ...
            bytes_printf([...
                '</SPARSE_DATA>\n' ...
                '<INDEX_LIST\n' ...
                '  data_type="Node_Bucket_node_indices"\n' ...
                '  sorted_node_def="No"\n' ...
                '  COLMS_RANGE="1 8 0 3"\n' ...
                '  COLMS_LABS="Node Indices"\n' ...
                '  COLMS_TYPE="Node_Index"\n' ...
                '  ni_type="int"\n' ...
                ni_form ...
                '  ni_dimen="6" >']) ...
            node_data ...
            bytes_printf([...
                '</INDEX_LIST>\n' ])...
            ];

function b=bytes_printf(varargin)
    b=uint8(sprintf(varargin{:}));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [tests,pass,perf]=test_mrgann(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
test_string={'wqrs(''mitdb/100'');mrgann(''mitdb/100'',''atr'',''wqrs'',''testAnn'')'};

clean_up={['delete([pwd filesep ''mitdb'' filesep ''*'']);' ...
          'rmdir([pwd filesep ''mitdb''],''s'');';]};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 11-  System response to discrete time sinusoidal input

syms n
w0=2;
Hw0=(3*j*w0+2)/(w0+1);

mag=abs(Hw0);
phas=angle(Hw0);

y=3*mag*sin(w0*n+1+ phas);

n=0:50;
y=subs(y,n);
stem(n,y);
legend('y[n]');

author: Sable

repo: mcbench-benchmarks


text:
function test_ft_sourcewrite

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_sourcewrite

gridsize = 38*48*41;
nsamples = 20;
source = [];
source.dim = [38,48,41];
source.pos = randn(gridsize,3);
source.time = 1:nsamples;
source.pow = randn(gridsize,nsamples);
source.inside = 1:gridsize;

cfg = [];
cfg.parameter = 'pow';
cfg.filename = tempname;
cfg.filetype = 'nifti';
ft_sourcewrite(cfg, source)
delete([cfg.filename '.nii'])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Isoparametric Formulation Implementation
% clear memory
clear all
close all
clc
% E: modulus of elasticity
% A: area of cross section
% L: length of bar
E=8; L=4;
u_exact=@(x) (56-8*(x-2)-24*heaviside(x-5))/2/x;
hold on;
ezplot(u_exact,[2 6])
title('Exact Solution v.s. FEM','interpreter','latex');
xlabel('x','interpreter','latex');
ylabel('Axial stress, $\it{\sigma}_{x}$','interpreter','latex','FontSize',12);

for i=1:4   %For NEL = 1, 2, 4, 8
fprintf( '\nNumber of elements:%d\n\n',2^(i-1) );
% numberElements: number of elements
numberElements=2^(i-1);
% numberNodes: number of nodes
numberNodes=2*numberElements+1;
A=zeros(1,numberElements);
% generation of coordinates and connectivities
NNOD=3;
nodeCoordinates=linspace(2,L+2,numberNodes);
%Generate element length vector
for i=1:numberElements
    Le(i)=L/numberElements;
    elementNodes(i,:)=[(i-1)*2+1 (i-1)*2+2 (i-1)*2+3];
    A(i)=(nodeCoordinates(2*i-1)+Le(i)/2)*2;
end
% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  detJacobian=Le(e)/2;
  invJacobian=1/detJacobian;
  ngp = 3;
  [w,xi]=gauss1d(ngp);
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(end)));
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL3(xi(ip));
      B=naturalDerivatives*invJacobian;
      stiffness(elementDof,elementDof)=...
      stiffness(elementDof,elementDof)+ B'*B*w(ip)*detJacobian*E*A(e);
      force(elementDof)=force(elementDof)+...
          8*shape'*detJacobian*w(ip);
  end
  if(nodeCoordinates(elementDof(end))==5)
     x=(5-xc)/detJacobian;
     [s,n]=shapeFunctionL3(x);
     force(elementDof)= force(elementDof)+...
         24*s';
  end
  if(nodeCoordinates(elementDof(1))<5&&...
          nodeCoordinates(elementDof(3))>5)
     x=(5-xc)/detJacobian;
     [s,n]=shapeFunctionL3(x);
     force(elementDof)= force(elementDof)+...
         24*s';
  end
end
% boundary conditions and solution
% prescribed dofs
prescribedDof=[1];
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactionsPretty(displacements,stiffness, ...
    numberNodes,prescribedDof,force)
fprintf('Axial stress\n')
fprintf('element\t\taxial stress\n')
%Stress and strain recovery
ngp = 3;
elementNodeCoor=zeros(numberElements*NNOD,1);
elementNodeStr=zeros(numberElements*NNOD,1);
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:);
  detJacobian=Le(e)/2;
  invJacobian=1/detJacobian;
  xi=[-1 0 1];
  for ip=1:NNOD;
      [shape,naturalDerivatives]=shapeFunctionL3(xi(ip));
      B=naturalDerivatives*invJacobian;
      elementNodeCoor(ip+(e-1)*NNOD,1)=nodeCoordinates(elementDof(ip));
      elementNodeStr(ip+(e-1)*NNOD,1)=E*B*displacements(elementDof,1);
      fprintf('%2.0f\t%2.0fth node\t%10.4e\n', e, ip,elementNodeStr(ip+(e-1)*NNOD,1))
  end
end

%post process
switch numberElements
    case 1
        str='r*--';
    case 2
        str='g*--';
    case 4
        str='k*--';
    case 8
        str='m*--';
end

plot(elementNodeCoor,elementNodeStr,str)
hold on;
end
legend('Exact solution','NEL=1','NEL=2','NEL=4','NEL=8','interpreter','latex');

repo: Aero-matlab


text:
% demo_ECG_sleep_staging_slpdb
% Note: Raw slpdb data should be converted to Matlab format by wfdb2mat function before analysis.

fid=fopen('RECORDS');
cc=textscan(fid,'%s');
cc=cc{1};
fclose(fid);

class=1;
% class = 1, four types classification
% class = 4, two types classification
%   sleep_stage:  = 1: Wake,
%                   2: REM sleep,
%                   3: NREM light sleep,
%                   4: NREM deep sleep, when class = 1
%   sleep_stage:  = 1: NREM sleep,
%                   2: Wake+REM, when class = 4

% par
for i=1:length(cc)
    do_ECG_sleep_staging_slpdb_sub(cc{i},class)
end

function do_ECG_sleep_staging_slpdb_sub(record,class)

record
[tm,signal,Fs,siginfo]=rdmat([record 'm']);
ECGdata=signal(:,1);
fs=250;
% class=1;

sleep_stage = ECG_sleep_staging_prediction(ECGdata,fs,class);

% read stage type from annotation
[ann,type,subtype,chan,num,comments]=read_ann(record,'st');
ann_epoch_len=30; % annotation epoch is 30s
seg_length=300; % segment length 300s (5 min)
stage=[];
stage(1:round(ann(end)/Fs/ann_epoch_len)+1)=NaN;

for j=1:length(ann)
    time_j=round(ann(j)/Fs/ann_epoch_len)+1;
    switch (comments{j,1}(1))
        case 'W'
            stage(time_j)=1;
        case 'R'
            stage(time_j)=2;
        case '1'
            stage(time_j)=3;
        case '2'
            stage(time_j)=3;
        case '3'
            stage(time_j)=4;
        case '4'
            stage(time_j)=4;
    end
end
if class==4
    stage(find(stage==1))=2;
    stage(find(stage>2))=1;
end
save(['ECG_sleep_staging_result_class_' num2str(class) '_' record],'sleep_stage','stage');
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
IN_DIR = 'images/stacked';
OUT_DIR = 'images/detected';
mkdir(OUT_DIR);
N = 88;

for i = 1 : N
    fn = sprintf('%s/%d.jpg', IN_DIR, i);
    I = imread(fn);
    detectEllipses(I);
    saveas(gcf, sprintf('%s/%d.png', OUT_DIR, i));
    close;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [convRVTOut, rvtOut, verbose] = tapas_physio_create_rvt_regressors(...
    ons_secs, sqpar, model_rvt, verbose)
% computes respiratory response function regressor and respiratory volume per time
%
%    [convRVTOut, rvtOut, verbose] = tapas_physio_create_rvt_regressors(...
%                                   ons_secs, sqpar, model_rvt, verbose)
% References:
%
%   Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008.
%       The respiration response function: The temporal dynamics of
%       fMRI signal fluctuations related to changes in respiration.
%       NeuroImage 40, 644-654.
%
%   Harrison, S.J., Bianchi, S., Heinzle, J., Stephan, K.E., Iglesias, S.,
%   Kasper L., 2021.
%   A Hilbert-based method for processing respiratory timeseries.
%   NeuroImage, 117787. https://doi.org/10.1016/j.neuroimage.2021.117787

% IN
%   ons_secs.
%       ons_secs            ons_secs structure with variable `fr`
%                           (filtered respiratory signal time series)
%       sqpar               scan timing information (sequence parameters)
%                           slice onsets etc.
%       model_rvt           rvt modeling parameter structure. e.g.
%                           model_rvt.method
%                           'hilbert' (default, [Harrison2021]) or
%                           'peaks' [Birn2006]
%
% OUT
%   convRVTOut          [nScans, nDelays, nSampleSlices]
%                       respiratory response function regressor after
%                       convolution for specified delays and downsampled
%                       to given slices.
% EXAMPLE
%   [convHRV, hr] = tapas_physio_create_hrv_regressor(physio_out.ons_secs, physio_out.sqpar);
%
%   See also tapas_physio_rvt_hilbert tapas_physio_rvt_peaks tapas_physio_rrf

% Author: Lars Kasper
% Created: 2014-01-20
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering,
%               University of Zurich and ETH Zurich.
%
% This file is part of the physIO toolbox, which is released under the
% terms of the GNU General Public Licence (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.

if nargin < 3
    physio = tapas_physio_new;
    model_rvt = physio.model.rvt;
end

delays = model_rvt.delays;

if nargin < 4
    verbose.level = [];
    verbose.fig_handles = [];
end

slicenum = 1:sqpar.Nslices;

% Calculate RVT
sample_points  = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum);
switch lower(model_rvt.method)
    case 'peaks'
        [rvt, ~, ~, verbose] = tapas_physio_rvt_peaks(ons_secs.fr, ons_secs.t, sample_points, verbose);
    case 'hilbert'
        [rvt, verbose] = tapas_physio_rvt_hilbert(ons_secs.fr, ons_secs.t, sample_points, verbose);
    otherwise
        error('Unrecognised value for ''rvt.method'' (%s)!', model_rvt.method)
end
rvt = rvt / max(abs(rvt)); % normalize for reasonable range of regressor

if verbose.level >=2
    verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params();
    set(gcf, 'Name', 'Model: Convolution Respiration RVT X RRF');
    subplot(2,2,1)
    plot(sample_points,rvt, 'g');xlabel('time (seconds)');
    title('Respiratory volume per time');
    ylabel('a.u.');
end

% Generate RRF
dt = sqpar.TR / sqpar.Nslices;
t = 0:dt:60;  % seconds
rrf = tapas_physio_rrf(t);
rrf = rrf / max(abs(rrf));

if verbose.level >= 2
    subplot(2,2,2)
    plot(t, rrf,'g'); xlabel('time (seconds)');
    title('Respiratory response function');
end

% Convolve and rescale for display purposes
convRVT = tapas_physio_conv(rvt, rrf, 'causal');
convRVT = convRVT / max(abs(convRVT));

if verbose.level >= 2
    subplot(2,2,3)
    plot(sample_points, convRVT,'g');xlabel('time (seconds)');
    title('Resp vol time X resp response function');
end

% Create shifted regressors convolved time series, which is equivalent to
% delayed response functions according to Wikipedia (convolution)
%
% "Translation invariance[edit]
% The convolution commutes with translations, meaning that
%
% \tau_x ({f}*g) = (\tau_x f)*g = {f}*(\tau_x g)\,
% where \tau_x is the translation of the function f by x defined by
% (\tau_x f)(y) = f(y-x).

% remove mean and linear trend to fulfill periodicity condition for
% shifting
convRVT = detrend(convRVT);

% TODO: what happens at the end/beginning of shifted convolutions?
nDelays = numel(delays);
nShiftSamples = ceil(delays/dt);

% resample to slices needed
nSampleSlices = numel(sqpar.onset_slice);
nScans = numel(sample_points(sqpar.onset_slice:sqpar.Nslices:end));

rvtOut = zeros(nScans,nSampleSlices);
convRVTOut = zeros(nScans,nDelays,nSampleSlices);
samplePointsOut = zeros(nScans,nSampleSlices);

for iDelay = 1:nDelays
    convRVTShifted = circshift(convRVT, nShiftSamples(iDelay));
    for iSlice = 1:nSampleSlices
        onset_slice = sqpar.onset_slice(iSlice);
        rvtOut(:,iSlice) = rvt(onset_slice:sqpar.Nslices:end)';
        convRVTOut(:,iDelay,iSlice) = convRVTShifted(onset_slice:sqpar.Nslices:end);
        samplePointsOut(:,iSlice) = sample_points(onset_slice:sqpar.Nslices:end);
    end
end

if verbose.level >= 2
    subplot(2,2,4)
    [tmp, iShiftMin] = min(abs(delays));
    hp{1} = plot(samplePointsOut, rvtOut,'k--');hold all;
    hp{2} = plot(samplePointsOut, squeeze(convRVTOut(:,iShiftMin,:)),'g');
    xlabel('time (seconds)');
    title('RVT regessor');
    legend([hp{1}(1), hp{2}(1)], 'respiratory volume / time (a. u.)', ...
        'respiratory response regressor');
end

end

author: translationalneuromodeling

repo: tapas


text:
function [A,b,x] = spikes(n,t_max)
%SPIKES Test problem with a "spiky" solution.
%
% [A,b,x] = spikes(n,t_max)
%
% Artificially generated discrete ill-posed problem.
%
% The solution x consists of a unit step at t = .5, and a pulse train
% of spikes of decrasing magnitude at t = .5, 1.5, 2.5, ...
%
% The parameter t_max is optional; its default value is 5.
% It controls the length of the pulse train.

% Per Christian Hansen, IMM, 04/21/97.

% Initialization.
if (nargin == 1), t_max = 5; end
del = t_max/n;

% Compute the matrix A.
[t,sigma] = meshgrid(del:del:t_max,del:del:t_max);
A = sigma./(2*sqrt(pi*t.^3)).*exp(-(sigma.^2)./(4*t));

% Compute the right-hand side b and the solution x.
if (nargout > 1)
  heights = 2*ones(t_max,1); heights(1) = 25;
  heights(2) = 9; heights(3) = 5; heights(4) = 4; heights(5) = 3;
  x = zeros(n,1); n_h = 1;
  peak = 0.5/t_max; peak_dist = 1/t_max;
  if (peak < 1)
    n_peak = round(peak*n); x(n_peak) = heights(n_h);
    x(n_peak+1:n) = ones(n-n_peak,1);
    peak = peak + peak_dist; n_h = n_h + 1;
  end
  while (peak < 1)
    x(round(peak*n)) = heights(n_h);
    peak = peak + peak_dist; n_h = n_h + 1;
  end
  b = A*x;
end

repo: REST


text:
function allCombos=genAllCapacitatedCombos(numPerBin,numBins,totalItems)
%%GENALLCAPACITATEDCOMBOS Given totalItems items, generate all combinations
%               combinations of placing the items into numBins unlabeled
%               bins where each bin can hold numPerBin items. If
%               numPerBin=1, then this is just generating standard
%               combinations of numBins items taken from a total of
%               totalItems. The capacitated combinations are equivalent to
%               putting totalItems labeled balls into numBins unlabeled
%               bins such that each bin holds numPerBin items and the
%               ordering of the items in the bins does not matter.
%
%INPUTS: numPerBin The scalar integer capacity of the bins; numPerBin>=1.
%          numBins The scalar integer number of bins; numBins>=1.
%       totalItems The scalar integer total number of items to place in the
%                  bins; totalitems>=numBins*numPerBin. If omitted or an
%                  empty matrix is passed, then the default of
%                  totalItems=numBins*numPerBin is used.
%
%OUTPUTS: allCombos The numPerBinXnumBinsXnumCapCombos set of capacitated
%                   combinations. allCombos(:,:,i) is the ith capacitated
%                   combination. The rows are the items in each bin and the
%                   columns are the bins. Items are numbered from 1 to
%                   totalItems.
%
%If we assume that totalItems=n=numPerBin*numBins, then we can generate the
%combinations by fixing the first element of each bin to the smallest index
%unassigned to any previous bin in the current combination being
%constructed. Thus, the first bin will always contain item 1. The second
%bin will always contain the smallest unassigned item that is not in the
%first bin and so on. Thus, we have to go through standard
%(n-(k-1)*numPerBin-1) choose (numPerBin-1) combinations of the remaining
%indices for the kth bin conditioned on the previous bins. This function
%implements a non-recursive algorithm for generating such combinations.
%
%In the general case, if totalItems>n, then we have to additionally go
%through all possible combinations of n items chosen from totalItems. This
%means generating all combinations as if totalItems=n and then remapping
%the indices for each subset of n items from totalItems.
%
%A total of numCapacitatedCombos(numPerBin,numBins,totalItems) capacitated
%combinations will be produced. Note that the number of capacitated
%combinations can grow rapidly and one might prefer visiting them
%sequentially using getNextCapacitatedCombo rather than generating all of
%them at once.
%
%EXAMPLE:
%Here we put six items into two bins of capacity 3:
% allCombos=genAllCapacitatedCombos(3,2,6)
%One wil get 10 combinations:
% allCombos(:,:,1)=[1, 4;
%                   2, 5;
%                   3, 6];
% allCombos(:,:,2)=[1, 3;
%                   2, 5;
%                   4, 6];
% allCombos(:,:,3)=[1, 3;
%                   2, 4;
%                   5, 6];
% allCombos(:,:,4)=[1, 3;
%                   2, 4;
%                   6, 5];
% allCombos(:,:,5)=[1, 2;
%                   3, 5;
%                   4, 6];
% allCombos(:,:,6)=[1, 2;
%                   3, 4;
%                   5, 6];
% allCombos(:,:,7)=[1, 2;
%                   3, 4;
%                   6, 5];
% allCombos(:,:,8)=[1, 2;
%                   4, 3;
%                   5, 6];
% allCombos(:,:,9)=[1, 2;
%                   4, 3;
%                   6, 5];
% allCombos(:,:,10)=[1, 2;
%                    5, 3;
%                    6, 4];
%
%July 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The total number of items to assign in combinations across the bins.
n=numBins*numPerBin;

if(nargin<3||isempty(totalItems))
    totalItems=n;
end

if(n>totalItems)
    error('totalItems cannot be less than numBins*numPerBin')
end

%Special case for standard combinations.
if(numPerBin==1)
    allCombos=genAllCombinations(totalItems,numBins,1,1);
    numCopies=size(allCombos,2);
    allCombos=reshape(allCombos,[numPerBin,numBins,numCopies]);
    return
end

%The total number of items to assign in combinations.
n=numBins*numPerBin;

%The total number of capacitated combinations if n=totalItems.
numCombosNoExtra=numCapacitatedCombos(numPerBin,numBins);

%If n~=totalItems, then this is the number of variants of each capacitated
%combination that will have to be generated.
numCopies=binomial(totalItems,n);

%To store the results.
allCombos=zeros(numPerBin,numBins,numCopies,numCombosNoExtra);

numPerBinM1=numPerBin-1;

curIdxCombos=zeros(numPerBin,numBins);
curCombos=zeros(numPerBin-1,numBins);

%This holds the indices that can be assigned in each bin and in subsequent
%bins.
idxLists=zeros(n,numBins);
idxLists(:,1)=1:n;

curSet=1;
curBin=1;
backtracking=false;
while(curBin>0)
    %This is the number of elements that can be assigned in the current
    %bin including making a hard assignment on the first index.
    numInIdx=n-(curBin-1)*numPerBin;

    if(backtracking==false)
        %We are just entering this bin.
        %The first element in each bin is fixed to the first index in the
        %idxList.
        curIdxCombos(1,curBin)=idxLists(1,curBin);

        %In this level, we have to go through all combinations of the
        %remaining n-curBin*(numPerBin-1); unassigned elements of
        %idxLists(:,curBin) in this bin. Each time an assignment is made in
        %this bin, those elements are removed from idxLists going into the
        %next higher level so that we never get repeats.
        curIdxCombos(2:numPerBin,curBin)=idxLists(2:numPerBin,curBin);
        %Current combination of the indices in idxLists in this bin.
        curCombos(:,curBin)=1:numPerBinM1;

        if(curBin==numBins)
            %If we are in the final bin, then save the combo with all
            %variants of indices and move to the previous bin.

            I=1:n;
            for curCopy=1:numCopies
                allCombos(:,:,curCopy,curSet)=I(curIdxCombos);

                I=getNextCombo(I,totalItems,1);
            end

            curSet=curSet+1;

            backtracking=true;
            curBin=curBin-1;
            continue;
        else
            %Otherwise, set idxLists for the next level and move to the
            %next bin.numInIdxNext
            numInIdxNext=numInIdx-numPerBin;
            idxLists(1:numInIdxNext,curBin+1)=idxLists((numPerBin+1):numInIdx,curBin);

            curBin=curBin+1;
        end
    else
        %If here, we are backtracking.

        %We subtract 1 because we have made a hard assignment on the first
        %index.
        nextCombo=getNextCombo(curCombos(:,curBin),numInIdx-1,1);
        if(isempty(nextCombo))
            %Continue backtracking.
            curBin=curBin-1;
            continue;
        end
        %Otherwise, record the current combination.
        curCombos(:,curBin)=nextCombo;

        %Record the indices for the current combination.
        nextCombo=nextCombo+1;
        curIdxCombos(2:numPerBin,curBin)=idxLists(nextCombo,curBin);

        %Now, set the values in idxLists at the next bin level to those
        %that are in idxLists but not in nextCombo. We note that the
        %indices are given in increasing order and that the first index is
        %skipped because we know it is assigned.
        curAssignedIdx=2;
        curNextIdx=1;
        for curListIdx=2:numInIdx
            if(idxLists(curListIdx,curBin)==curIdxCombos(curAssignedIdx,curBin))
                curAssignedIdx=curAssignedIdx+1;
                if(curAssignedIdx>numPerBin)
                    %Then all of the rest of the indices are assigned.
                    numInIdxNext=numInIdx-numPerBin;
                    idxLists(curNextIdx:numInIdxNext,curBin+1)=idxLists((curListIdx+1):numInIdx,curBin);
                    break;
                end
            else
                idxLists(curNextIdx,curBin+1)=idxLists(curListIdx,curBin);
                curNextIdx=curNextIdx+1;
            end
        end

        %Go to the next bin.
        curBin=curBin+1;
        backtracking=false;
    end
end

allCombos=allCombos(:,:,:);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [varargout]=sv3(varargin)

% function [hf]=sv3(M,v,optionStruct)
% ------------------------------------------------------------------------
% sv3 (slice view 3D) is a 3D slice viewer function. The 3D image data M is
% rendered using 3 mutally orthogonal slices. The option input v is the
% voxel size which scales the image appropriately in the 3 directions. The
% 3rd optional input vixOptStruct which can be used to make custom
% visualization settings. The default structure containts the following:
%
% optionStructDefault.colormap=gray(250); %colormap
% optionStructDefault.clim=[min(M(~isnan(M))) max(M(~isnan(M)))]; %color limits
% optionStructDefault.fontColor='w'; %font color
% optionStructDefault.fontSize=20; %font size
% optionStructDefault.edgeColor='none';
% optionStructDefault.figStruct=figStruct; %figure options (see cFigure)
% optionStructDefault.sliceIndices=round(size(M)/2); %Default mid-slices
% optionStructDefault.alphaLevel=1;
% optionStructDefault.origin=[0 0 0];
% optionStructDefault.updateFrequency=10;
%
% See also: sliceViewer, sv2, imx
%
% Change log:
% 2018/06/06 Added initial slice indices as option to input structure
% 2018/06/06 Added basic description at the top of this function
% 2019/08/09 Changed to use uicontrol slider rather than java slider due
% to future removal of javacomponent
% 2021/11/03 Fixed issue with specification of image origin and shifting
% slider
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        M=varargin{1};
        v=ones(1,3);
        optionStruct=[];
    case 2
        M=varargin{1};
        v=varargin{2};
        optionStruct=[];
    case 3
        M=varargin{1};
        v=varargin{2};
        optionStruct=varargin{3};
end

%Expand voxel size if it is scalar
if numel(v)==1
    v=v.*ones(1,3);
end

M=double(M); %Conver the image to a double

figStruct.Name='GIBBON: Slice viewer'; %Figure name
figStruct.Color='k'; %Figure background color
figStruct.ColorDef='black'; %Setting colordefinitions to black

optionStructDefault.colormap=gray(250); %colormap
optionStructDefault.clim=[min(M(~isnan(M))) max(M(~isnan(M)))]; %color limits
optionStructDefault.fontColor='w'; %font color
optionStructDefault.fontSize=20; %font size
optionStructDefault.edgeColor='none';
optionStructDefault.figStruct=figStruct; %figure options (see cFigure)
optionStructDefault.sliceIndices=round(size(M)/2); %Default mid-slices
optionStructDefault.alphaLevel=1;
optionStructDefault.origin=[0 0 0];
optionStructDefault.updateFrequency=10;

[optionStruct]=structComplete(optionStruct,optionStructDefault,1);

M=double(M);

%%
% Plot settings
scrollBarWidth=20; %Scrollbar width
fontColor=optionStruct.fontColor;
fontSize=optionStruct.fontSize;
cMap=optionStruct.colormap;
cLim=optionStruct.clim;
figStruct=optionStruct.figStruct;
alphaLevel=optionStruct.alphaLevel;
originLoc=optionStruct.origin;
updateFrequency=optionStruct.updateFrequency;

if all(isnan(cLim)) | all(isinf(cLim))
    cLim=[-1 1];
end

if diff(cLim)<eps
   cLim=cLim+[-1 1];
end

%%

%Defining row, column and slice indicices for slice patching
sliceIndexI=optionStruct.sliceIndices(1); %(close to) middle row
sliceIndexJ=optionStruct.sliceIndices(2); %(close to) middle column
sliceIndexK=optionStruct.sliceIndices(3); %(close to) middle slice

%%

[ax,ay,az]=im2cart([size(M,1)+1 0],[size(M,2)+1 0],[size(M,3)+1 0],v);

navString=['I: ',num2str(sliceIndexI),', J:  ',num2str(sliceIndexJ),', K: ',num2str(sliceIndexK)];

hf=cFigure(figStruct);
ht=gtitle(navString);
axis equal; axis tight; view(3);  axis vis3d; axis([ax(2) ax(1) ay(2) ay(1) az(2) az(1)]); grid on; box on; hold on;
colormap(cMap); colorbar;
caxis(cLim);
set(gca,'fontSize',fontSize);
drawnow;

%%

%Initialize sliders
hSlider_I= uicontrol(hf,'Style','slider','Position',[0,0,scrollBarWidth,round(hf.Position(4))]);
set(hSlider_I,'Value',sliceIndexI,'Min',1,'Max',size(M,1),'SliderStep',[1/(size(M,1)-1) 1/(size(M,1)-1)]);
hSlider_I.Callback={@plotSlice,{hf,hSlider_I,1}};
addlistener(hSlider_I,'ContinuousValueChange',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_I,1}));
addlistener(hSlider_I,'Value','PostSet',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_I,1}));

hSlider_J= uicontrol(hf,'Style','slider','Position',[1*scrollBarWidth,0,scrollBarWidth,round(hf.Position(4))]);
set(hSlider_J,'Value',sliceIndexJ,'Min',1,'Max',size(M,2),'SliderStep',[1/(size(M,2)-1) 1/(size(M,2)-1)]);
hSlider_J.Callback={@plotSlice,{hf,hSlider_J,2}};
addlistener(hSlider_J,'ContinuousValueChange',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_J,2}));
addlistener(hSlider_J,'Value','PostSet',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_J,2}));

hSlider_K= uicontrol(hf,'Style','slider','Position',[2*scrollBarWidth,0,scrollBarWidth,round(hf.Position(4))]);
set(hSlider_K,'Value',sliceIndexK,'Min',1,'Max',size(M,3),'SliderStep',[1/(size(M,3)-1) 1/(size(M,3)-1)]);
hSlider_K.Callback={@plotSlice,{hf,hSlider_K,3}};
addlistener(hSlider_K,'ContinuousValueChange',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_K,3}));
addlistener(hSlider_K,'Value','PostSet',@(hObject, event) plotSlice(hObject,event,{hf,hSlider_K,3}));

%% Set resize function

set(hf,'ResizeFcn',{@setScrollSizeFunc,{hf,scrollBarWidth,hSlider_I,hSlider_J,hSlider_K}});

%%

t=clock;
hf.UserData.sv3.time=t;
hf.UserData.sv3.Name=figStruct.Name;
hf.UserData.sv3.M=M;
hf.UserData.sv3.v=v;
hf.UserData.sv3.patchTypes={'si','sj','sk'};
hf.UserData.sv3.sliceIndices=[sliceIndexI sliceIndexJ sliceIndexK];
hf.UserData.sv3.fontColor=fontColor;
hf.UserData.sv3.fontSize=fontSize;
hf.UserData.sv3.edgeColor=optionStruct.edgeColor;
hf.UserData.sv3.hp=nan(3,1);
hf.UserData.sv3.ht=ht;
hf.UserData.sv3.M_plot=M;
hf.UserData.sv3.alphaLevel=alphaLevel;
hf.UserData.sv3.origin=originLoc;
hf.UserData.sv3.sliderHandles=[hSlider_I,hSlider_J,hSlider_K];
hf.UserData.sv3.updateFrequency=updateFrequency;

%%
set(hSlider_I,'Value',sliceIndexI);
set(hSlider_J,'Value',sliceIndexJ);
set(hSlider_K,'Value',sliceIndexK);

%Initialize view
pause(1/updateFrequency); %Wait so plot will update
plotSlice([],[],{hf,hSlider_I,1});
hf.UserData.sv3.time=t; %Reset clock so this happens now
plotSlice([],[],{hf,hSlider_J,2});
hf.UserData.sv3.time=t; %Reset clock so this happens now
plotSlice([],[],{hf,hSlider_K,3});

drawnow;

%%
varargout{1}=hf;

end

function plotSlice(~,~,inputCell)

hf=inputCell{1};

jSlider=inputCell{2};
dirOpt=inputCell{3};
sliceIndex=round(get(jSlider,'Value'));
hf.UserData.sv3.sliceIndices(dirOpt)=sliceIndex;
sliceIndices=hf.UserData.sv3.sliceIndices;

dt=1/hf.UserData.sv3.updateFrequency;
t=hf.UserData.sv3.time;
t2=clock;
dtt=etime(t2,t); %Elapsed time

if dtt>dt %If ready to update
    hf.UserData.sv3.time=t2;

    M=hf.UserData.sv3.M_plot;
    v=hf.UserData.sv3.v;
    patchType=hf.UserData.sv3.patchTypes{dirOpt};

    logicPatch=false(size(M));
    switch dirOpt
        case 1
            logicPatch(sliceIndex,:,:)=1;
        case 2
            logicPatch(:,sliceIndex,:)=1;
        case 3
            logicPatch(:,:,sliceIndex)=1;
    end

    figure(hf); %TEMP FIX for bug in MATLAB 2018

    if isnan(hf.UserData.sv3.hp(dirOpt))
        [F,V,C]=im2patch(M,logicPatch,patchType);
        [V(:,1),V(:,2),V(:,3)]=im2cart(V(:,2),V(:,1),V(:,3),v);
        V=V+hf.UserData.sv3.origin(ones(size(V,1),1),:);
        hf.UserData.sv3.hp(dirOpt)= gpatch(F,V,C,hf.UserData.sv3.edgeColor,hf.UserData.sv3.alphaLevel);
    else
        V=get(hf.UserData.sv3.hp(dirOpt),'Vertices');
        switch dirOpt
            case 1
                V(:,2)=(sliceIndex-0.5).*v(1)+hf.UserData.sv3.origin(2);
            case 2
                V(:,1)=(sliceIndex-0.5).*v(2)+hf.UserData.sv3.origin(1);
            case 3
                V(:,3)=(sliceIndex-0.5).*v(3)+hf.UserData.sv3.origin(3);
        end
        set(hf.UserData.sv3.hp(dirOpt),'CData',M(logicPatch)); %Set color data
        set(hf.UserData.sv3.hp(dirOpt),'Vertices',V); %Set vertices
    end

    navString=['I: ',num2str(sliceIndices(1)),', J:  ',num2str(sliceIndices(2)),', K: ',num2str(sliceIndices(3))];

    set(hf.UserData.sv3.ht,'string',navString);

    hf.Name=[hf.UserData.sv3.Name,' ',navString];
end
end

function setScrollSizeFunc(~,~,inputCell)
hf=inputCell{1};
w=inputCell{2};

for q=3:numel(inputCell)
    hSlider=inputCell{q};
    posData=[w*(q-3),0,w,round(hf.Position(4))];
    set(hSlider,'Position',posData);
end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ mColumnVectorImage ] = ImageToColumnsSliding( mInputImage, vBlockSize )
% ----------------------------------------------------------------------------------------------- %
% [ mColumnImage ] = ImageToColumns( mInputImage, blockRadius )
%   Creates an column image from the sliding neighborhood in mInpuImage
% Input:
%   - mInputImage           -   Input image.
%                               Structure: Image Matrix (Single Channel)
%                               Type: 'Single' / 'Double'.
%                               Range: [0, 1].
%   - vBlockSize            -   Block Size.
%                               Structure: 2D Vector.
%                               Type: 'Single' / 'Double'.
%                               Range: [0, 1].
% Output:
%   - mColumnVectorImage    -   Column Vector Image.
%                               Structure: Image Matrix (Single Channel)
%                               Type: 'Single' / 'Double'.
%                               Range: [0, 1].
% Remarks:
%   1.  Prefixes:
%       -   'm' - Matrix.
%       -   'v' - Vector.
%   2.  Converts each sliding `vBlockSize(1)` by `vBlockSize(2) block of
%       `mInputImage` into a column of `mColumnVectorImage` with no zero
%       padding.
%   3.  Shouldn't be used for images larger than 400x400 and blocks of
%       51x51.
% TODO:
%   1.  I
%   Release Notes:
%   -   1.0.000     20/03/2015  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

[numRows, numCol] = size(mInputImage);
blockNumRows = vBlockSize(1);
blockNumCols = vBlockSize(2);

% Create Hankel-like indexing sub matrix.
nc = numRows - blockNumRows + 1;
nn = numCol - blockNumCols + 1;

vColumnIdx = [(0:(blockNumRows - 1))]';
vRowIdx = [1:nc];

t = vColumnIdx(:, ones(nc, 1)) + vRowIdx(ones(blockNumRows, 1), :);    % Hankel Subscripts
tt = zeros(blockNumRows * blockNumCols, nc);
rows = 1:blockNumRows;
for ii = 0:(blockNumCols - 1)
    tt(((ii * blockNumRows) + rows), :) = t + (numRows * ii);
end
mColumnVectorImageIdx = zeros((blockNumRows * blockNumCols), (nc * nn));
cols = 1:nc;
for jj = 0:(nn - 1)
    mColumnVectorImageIdx(:, ((jj * nc) + cols)) = tt + (numRows * jj);
end

mColumnVectorImage = mInputImage(mColumnVectorImageIdx);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function cpf(fig_no,target_name,format)
% This function cpf.m (change picture format)
% takes a picture in one format
% and saves it in another format.
%
% Syntax: load_pics(fig_no,target_name,format)
% e.g.  : load_pics(1,'holiday_on_Bornholm.bmp','epsc2');
%
% fig no      : Figure number where the figure is displayed
% target_name : Filename of new file. target_name = ''
%               gives target_name=input_filename.
% format      : One of the output formats availeble in MATLAB,
%               type "help fileformats" or
%                    "help print"

if nargin ~= 3
 error('Please specify: fig no., target name and output format');
end

% Pick file
[filename, pathname,filterspec] = uigetfile('*.*');

if strcmp(target_name,'')
    target_name = filename(1:end-4);
end

% Load file as image
[A,map]=imread(filename);

% Display image
figure(fig_no)
clf
image(A);
colormap(map);
axis equal off

% Save in new format
saveas(fig_no,target_name,format);


text:
function pass = test_get(pref)
% Test GET.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) cos(x.*y.*z));
[core, cols, rows, tubes] = tucker(f);

pass(1) = norm([-1 1 -1 1 -1 1] - f.domain) < tol;

pass(2) = norm(core(:) - f.core(:)) < tol;

pass(3) = norm(cols - f.cols) < tol;

pass(4) = norm(rows - f.rows) < tol;

pass(5) = norm(tubes - f.tubes) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function plotAmps(vw, scans)
%
% plotAmps(vw)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in the current ROI.
%
% gmb   5/25/98
% jw    6/7/08 added std to user data

% Compute means across scans, for all pixels in the
% currently selected ROI.
if exist('scans', 'var')
    [meanAmps,meanPhs,seAmps,meanStd] = vectorMeans(vw, scans);
else
    [meanAmps,meanPhs,seAmps,meanStd] = vectorMeans(vw);
end
selectGraphWin

% Header
ROIname     = vw.ROIs(vw.selectedROI).name;
headerStr   = ['Mean of amplitudes, ROI ',ROIname];
set(gcf,'Name',headerStr);

% Plot the bar graph
clf
fontSize = 14;
meanAmps
meanStd

mybar(meanAmps,meanStd);
xlabel('Scan','FontSize',fontSize);
ylabel('Mean Amplitude','FontSize',fontSize);
ylim =get(gca,'YLim');
set(gca,'YLim',ylim*1.1);
set(gca,'FontSize',fontSize);

% Save the data in gca('UserData')
data.y = meanAmps;
data.s = meanStd;
set(gca,'UserData',data);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ncs] = read_neuralynx_ncs(filename, begrecord, endrecord)

% READ_NEURALYNX_NCS reads a single continuous channel file
%
% Use as
%   [ncs] = read_neuralynx_ncs(filename)
%   [ncs] = read_neuralynx_ncs(filename, begrecord, endrecord)

% Copyright (C) 2005-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_neuralynx_ncs.m 2885 2011-02-16 09:41:58Z roboos $

if nargin<2
  begrecord = 1;
end
if nargin<3
  endrecord = inf;
end

% the file starts with a 16*1024 bytes header in ascii, followed by a number of records
hdr = neuralynx_getheader(filename);
fid = fopen(filename, 'rb', 'ieee-le');

% determine the length of the file
fseek(fid, 0, 'eof');
headersize = 16384;
recordsize = 1044;
NRecords   = floor((ftell(fid) - headersize)/recordsize);

if NRecords>0
  if (ispc), fclose(fid); end
  % read the timestamp from the first and last record
  hdr.FirstTimeStamp = neuralynx_timestamp(filename, 1);
  hdr.LastTimeStamp  = neuralynx_timestamp(filename, inf);
  if (ispc), fid = fopen(filename, 'rb', 'ieee-le'); end
else
  hdr.FirstTimeStamp = nan;
  hdr.LastTimeStamp  = nan;
end

if begrecord==0 && endrecord==0
  % only read the header
elseif begrecord<1
  error('cannot read before the first record');
elseif begrecord>NRecords
  error('cannot read beyond the last record')
elseif endrecord>NRecords
  endrecord = NRecords;
end

if begrecord>=1 && endrecord>=begrecord
  % rewind to the first record to be read
  status = fseek(fid, headersize + (begrecord-1)*recordsize, 'bof');
  if status~=0
    error('cannot jump to the requested record');
  end

  numrecord    = (endrecord-begrecord+1);
  TimeStamp    = zeros(1,numrecord,'uint64');
  ChanNumber   = zeros(1,numrecord);
  SampFreq     = zeros(1,numrecord);
  NumValidSamp = zeros(1,numrecord);
  Samp         = zeros(512,numrecord);  % this allows easy reshaping into a 1xNsamples vector

  for k=1:numrecord
    % read a single continuous data record
    TimeStamp(k)    = fread(fid,   1, 'uint64=>uint64');
    ChanNumber(k)   = fread(fid,   1, 'int32');
    SampFreq(k)     = fread(fid,   1, 'int32');
    NumValidSamp(k) = fread(fid,   1, 'int32');
    Samp(:,k)       = fread(fid, 512, 'int16');
    % mark the invalid samples
    Samp((NumValidSamp+1):end,k) = nan;
  end

  % store the record data in the output structure
  ncs.TimeStamp    = TimeStamp;
  ncs.ChanNumber   = ChanNumber;
  ncs.SampFreq     = SampFreq;
  ncs.NumValidSamp = NumValidSamp;
  % apply the scaling factor from ADBitVolts and convert to uV
  ncs.dat          = Samp * hdr.ADBitVolts * 1e6;

end
fclose(fid);

% store the header info in the output structure
ncs.NRecords = NRecords;
ncs.hdr      = hdr;

author: PatternRecognition

repo: OpenBMI


text:
function [err,time,solver,eqn] = femPoisson3(node,elem,pde,bdFlag,option,varargin)
%% FEMPOISSON3 solve Poisson equation by various finite element methods
%
%   FEMPOISSON3 computes approximations to the Poisson equation on a
%   sequence of meshes obtained by uniform refinement of a input mesh.
%
% See also Poisson, crack, Lshape
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

%% Check input arguments
if nargin >=1 && ischar(node)
    option.elemType = node;
    clear node
end
if ~exist('node','var') || ~exist('elem','var')
    [node,elem] = cubemesh([0,1,0,1,0,1],0.25); % default mesh is a cube
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
    pde = sincosdata3;                          % default data
end
if ~exist('bdFlag','var')
    bdFlag = setboundary3(node,elem,'Dirichlet');
end

%% Parameters
option = femoption(option);
maxIt = option.maxIt;   maxN = option.maxN; L0 = option.L0;
elemType = option.elemType; refType = option.refType;

%% Generate an initial mesh
for k = 1:L0
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect3(node,elem,bdFlag);
    end
end

%% Initialize err
errL2 = zeros(maxIt,1);   errH1 = zeros(maxIt,1);
erruIuh = zeros(maxIt,1); errMax = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    switch elemType
        case 'P1'     % piecewise linear function P1 element
            [u,Du,eqn,info] = Poisson3(node,elem,pde,bdFlag,option);
        case 'CR'     % piecewise linear function CR element
            [u,Du,eqn,info] = Poisson3CR(node,elem,pde,bdFlag,option);
        case 'P2'     % piecewise quadratic function
            [u,Du,eqn,info] = Poisson3P2(node,elem,pde,bdFlag,option);
        case 'WG'     % weak Galerkin element
            [u,Du,eqn,info] = Poisson3WG(node,elem,pde,bdFlag,option);
    end
    % compute error
    tic;
    if isfield(pde,'Du')
        if ~isempty(Du)
            errH1(k) = getH1error3(node,elem,pde.Du,Du);
        else
            errH1(k) = getH1error3(node,elem,pde.Du,u);
        end
    end
    if isfield(pde,'exactu')
        errL2(k) = getL2error3(node,elem,pde.exactu,u);
        % interpolation
        switch elemType
            case 'P1'
                uI = Lagrangeinterpolate(pde.exactu,node,elem);
            case 'CR'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'CR',eqn.face);
            case 'P2'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'P2',eqn.edge);
            case 'WG'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'WG',eqn.face);
        end
        erruIuh(k) = sqrt((u-uI)'*eqn.A*(u-uI));
        errMax(k) = max(abs(u-uI));
    end
    errTime(k) = toc;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    if option.printlevel>1
        fprintf('Time to compute the error %4.2g s \n H1 err %4.2g    L2err %4.2g \n',...
            errTime(k),errH1(k), errL2(k));
    end
    % record solver information
    itStep(k) = info.itStep;
    stopErr(k) = info.stopErr;
    flag(k) = info.flag;
    % plot
    N(k) = length(u);
    if option.plotflag && N(k) < 2e3 % show mesh and solution for small size
       figure(1);  showresult3(node,elem,u);
    end
    if N(k) > maxN
        break;
    end
    % refine mesh
    tic;
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect3(node,elem,bdFlag);
    end
    meshTime(k) = toc;
end

%% Plot convergence rates
if option.rateflag
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.55,0.4]);
    subplot(1,2,1)
    showrate2(N(1:k),errH1(1:k),1,'-*','||Du-Du_h||',...
              N(1:k),errL2(1:k),1,'k-+','||u-u_h||');
    subplot(1,2,2)
    showrate2(N(1:k),erruIuh(1:k),1,'m-+','||Du_I-Du_h||',...
              N(1:k),errMax(1:k),1,'r-*','||u_I-u_h||_{\infty}');
end

%% Output
err = struct('N',N,'H1',errH1(1:k),'L2',errL2(1:k),...
             'uIuhH1',erruIuh(1:k),'uIuhMax',errMax(1:k));
time = struct('N',N,'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assmble',assembleTime(1:k),'mesh',meshTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error
ts = zeros(k,3); ts = char(ts);
display(' #Dof   ||u-u_h||     ||Du-Du_h||   ||DuI-Du_h||  ||uI-u_h||_{max}');
display([num2str(err.N) ts num2str(err.L2,'%0.5e') ts num2str(err.H1,'%0.5e')...
         ts num2str(err.uIuhH1,'%0.5e') ts num2str(err.uIuhMax,'%0.5e')]);


text:
function [sys,x0,str,ts]=ADRC_1(t,x,u,flag,r0,h,B01,B02,r)

switch flag
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2
        sys=mdlUpdate(x,u,r0,h,B01,B02);
    case 3
        sys=mdlOutputs(x,r);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=4;
    sizes.NumOutputs=2;
    sizes.NumInputs=3;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdate(x,u,r0,h,B01,B02)
    e1=x(1)-u(1);
    fh=fhan(e1,x(2),r0,h);
    sys(1)=x(1)+h*x(2);
    sys(2)=x(2)+h*fh;
    e2=x(3)-u(2);
    sys(3)=x(3)+h*(x(4)-B01*e2+u(3));
    sys(4)=x(4)+h*(-B02*e2);

function sys=mdlOutputs(x,r)
    e3=x(1)-x(3);
    sys(1)=r*e3-x(4);
    sys(2)=x(1);
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;

function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function r = rank(f, tol)
%RANK      Rank of a SEPARABLEAPPROX.
%   RANK(F) produces an estimate of the rank of the approximant F.
%
%   RANK(F, TOL) is the number of singular values of F greater than TOL*N, where
%   N is the first singular value of F.
%
% See also LENGTH.

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

% Empty check:
if ( isempty(f) )
    r = [];
    return
end

% HARD tolerance.
if ( nargin == 1 )
    tol = 0;
end

% Compute the singular values of f.
s = svd( f );

% Check for zero function.
if ( max(s) == 0  )
    r = 0;
else
    % r = no. of s.v. above relative tol.
    r = find(s/s(1) > tol, 1, 'last');
end

end


text:
function len = edgeLength3d(edge, varargin)
%EDGELENGTH3D Return the length of a 3D edge.
%
%   L = edgeLength3D(EDGE);
%   Returns the length of a 3D edge, with following representation:
%   [x1 y1 z1 x2 y2 z2].
%
%   Example
%     p1 = [1 1 1];
%     p2 = [3 4 5];
%     edge = createEdge3d(p1, p2);
%     edgeLength3d(edge)
%     ans =
%         5.3852
%
%   See also
%     edges3d, createEdge3d, drawEdge3d

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-08-29,    using Matlab 9.4.0.813654 (R2018a)
% Copyright 2018 INRA - Cepia Software Platform.

if nargin == 1
    dp = edge(:, 4:6) - edge(:, 1:3);
else
    dp = varargin{1} - edge;
end

len = hypot(hypot(dp(:,1), dp(:,2)), dp(:,3));

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%-------------------------------------------------------------------------
% Coupling Kinetic and Fluid System of Euler Equations
% To solve Shock Tube Problem
%
% Kinetic equation use SBBGK and Hydrodynamic equation use Roe Euler
%
% Based on:
%  [1] Pierre Degond, Giacomo Dimarco and Luc Mieussens
%      A moving interface method for dynamic kinetic-fluid coupling
%      Journal of Computation Physics 227(2007)1176-1208
%
% By Manuel Diaz and Yun-Da Tsai
% 007@IAM 25.01.2013
%-------------------------------------------------------------------------

clear all; clc; close all;

%% Global Variables
global CFL r_time theta dt dtdx nx
global w k nv
global gamma etpfix

%% Controling Parameters
name        ='SBBGK1d'; % Simulation Name
CFL         = 0.05;     % CFL condition
r_time      = 1/10000;  % Relaxation time
tEnd        = 0.04;      % End time
theta       = 0;        % {-1} BE, {0} MB, {1} FD.
quad        = 2;        % for NC = 1 , GH = 2
method      = 1;        % for TVD = 1, WENO3 = 2, WENO5 = 3
IC_case     = 1;        % IC: {1}Sod's, {2}LE, {3}RE, {4}DS, {5}SS, {6}Cavitation
plot_figs   = 1;        % 0: no, 1: yes please!
write_ans   = 0;        % 0: no, 1: yes please!
gamma       = 2;      % Ratio of specific heats
flxtype     = 2;        % {1} Roe, {2} LF, {3} LLF, {4} Upwind <-non-conservative!
etpfix      = 0.90;     % {#} Harten's sonic entropy fix value, {0} no entropy fix

%% Space Discretization
nx  = 100;                      % number of cells
x   = linspace(0,1,nx);         % Physical domain -x
dx  = max(x(2:end)-x(1:end-1)); % delta x

%% Load Initial Condition
[z0,ux0,t0,p0,rho0,E0] = SSBGK_IC1d(x,IC_case);

%% Load Initial Cut Function
 xa = 0.5;   % buffer left boundary
 xb = 0.5125;   % buffer right boundary
 h0  = cutfunc(x,xa,xb);  % Physical Cut Function

%  h0 = ones(size(x));
%  h0 = zeros(size(x));

%% Discretization of the Velocity Space
% Microscopic Velocity Discretization (using Discrete Ordinate Method)
% that is to make coincide discrete values of microscopic velocities with
% values as the value points for using a quadrature method, so that we can
% integrate the velocity probability distribution to recover our
% macroscopics properties.
switch quad

    case{1} % Newton Cotes Quadrature:
    V  = [-20,20];  % range: a to b
    nv = 200;       % nodes desired (may not the actual value)
    [c,w,k] = cotes_xw(V(1),V(2),nv,5); % Using Netwon Cotes Degree 5

    case{2} % Gauss Hermite Quadrature:
    nv = 100;          % nodes desired (the actual value)
    [c,w] = GaussHermite(nv); % for integrating range: -inf to inf
    k = 1;            % quadrature constant.
    w = w.*exp(c.^2); % weighting function of the Gauss-Hermite quadrature

    otherwise
        error('Order must be between 1 and 2');
end

%% Applying DOM
% The actual nv value will be computed using 'lenght' vector function:
nv = length(c);
% Remap velocity points
    c = repmat(c,1,nx);     w = repmat(w,1,nx);
% Remap classical IC
    [rho0,ux0,p0] = apply_DOM(rho0,ux0,p0,nv);
% Remap Semiclassical IC
    [z0,t0,E0] = apply_DOM(z0,t0,E0,nv);
% Remap h coeficient
  h0 = repmat(h0,nv,1);

%% Semi0-classical Equilibrium Distribution Function
M0 = f_equilibrium_1d(z0,ux0,c,t0,theta);

%% Load initial Conditions and Spliting of information
M = M0; rho = rho0; ux = ux0; t = t0; p = p0; z = z0; h = h0; E = E0;

%% Split ICs in R and L
rhor = h.*rho0;     rhol = (1-h).*rho0;
uxr = h.*ux0;       uxl = (1-h).*ux0;
pr = h.*p0;         pl = (1-h).*p0;

% [zr,~,tr,~] = macroproperties1d(rhor,rhor.*uxr,pr+0.5*rhor.*uxr.^2,nx,nv,theta);
% [zl,~,tl,~] = macroproperties1d(rhol,rhol.*uxl,pl+0.5*rhol.*uxl.^2,nx,nv,theta);
% Ml = f_equilibrium_1d(zl,uxl,c,tl,theta) ;
% fr = f_equilibrium_1d(zr,uxr,c,tr,theta);
 fr=h.*M;
 Ml=(1-h).*M;
%% Main Loop
% Compute next time step
 dt = dx*CFL/max(c(:,1));

% dt =1/20000;
time  = 0:dt:tEnd;
dtdx = dt/dx;
% Ml(isnan(Ml)) = 0;     fr(isnan(fr)) = 0;
f = fr + Ml; % computed here for ploting purposes

count = 1; % iteration counter
tic;
for tsteps = time

    CFL=dtdx*max(c(:,1));
% Plot IC
   if plot_figs == 1; surf(f); end;

% Compute vector 'q'
   q = [rhol(1,:) ; rhol(1,:).*uxl(1,:) ; pl(1,:)+0.5*rhol(1,:).*uxl(1,:).^2];

% Update Physical cut function 'h'
     h_next = h; % this means: fixed buffer assumption

%     Evaluate Modified Boltzmann BGK
    [fr_next] = ModSBBGK(h,h_next,M,fr,Ml,c,flxtype);

% Evaluate Modificed Roe Euler Solver
    [rhol,rhoul,El] = ModEuler(h,h_next,q,fr_next,c);

%     plot partial result
%     if plot_figs == 1;
%         subplot(1,3,1); plot(x,rhol,'o'); title('Density');
%         subplot(1,3,2); plot(x,rhoul./rhol,'o'); title('Velocity');
%         subplot(1,3,3); plot(x,El,'o'); title('Energy');
%     end

%     macroscopic properties
    [zl,uxl,tl,pl] = macroproperties1d(rhol,rhoul,El,nx,nv,gamma,theta);

    % Apply DOM
    [tl,zl,uxl] = apply_DOM(tl,zl,uxl,nv);

    % Compute Ml_next
     Ml_next = f_equilibrium_1d(zl,uxl,c,tl,theta) ;
%  ****************************************************

    % New time step info: sum left and righ values with 'NaN' filter sum
    % function.
    Ml_next(isnan(Ml_next)) = 0;     fr_next(isnan(fr_next)) = 0;
    % Total f
     f  = fr_next + Ml_next;

     % Update New macroquantity
   [rho,rhou,E] = macromoments1d(k,w,f,c);

    %Apply Neumann BC's in total variables
    f(:,1)  = f(:,2);                 f(:,end)  = f(:,end-1);
    rho(:,1)= rho(:,2);            rho(:,end)= rho(:,end-1);
    rhou(:,1) = rhou(:,2);       rhou(:,end) = rhou(:,end-1);
    E(:,1)  = E(:,2);                 E(:,end)  = E(:,end-1);

    [rho,rhou,E] = apply_DOM(rho,rhou,E,nv);

    % Recover Semiclassical Conditions for next time step
     [z,ux,t,p] = macroproperties1d(rho,rhou,E,nx,nv,gamma,theta);

%  update New equilibrium
     M=f_equilibrium_1d(z,ux,c,t,theta);
     fr=h.*M;

%  preparing new Ml
     rhol = (1-h).*rho;
     uxl = (1-h).*ux;
     pl = (1-h).*p;
     [zl,~,tl,~] = macroproperties1d(rhol,rhol.*uxl,pl+0.5*rhol.*uxl.^2,nx,nv,gamma,theta);
% update New Ml
     Ml=   f_equilibrium_1d(zl,uxl,c,tl,theta) ;
     Ml(isnan(Ml)) = 0;

% update counter
    count = count+1;

% update plot
    drawnow
end
toc;

% write/plot final output
%plot(x,r_total,'o');


text:
function [z,ASAcontrol] = convol(x,y,shift1,shift2,last)
%CONVOL Convolution sum
%   Z = CONVOL(X,Y) is the convolution sum of two vectors X and Y.
%
%   CONVOL(X) is the auto-convolution sum of a vector X.
%
%   CONVOL(X,Y,SHIFT1,SHIFT2) or CONVOL(X,SHIFT1,SHIFT2) returns only the
%   elements from SHIFT1 up to SHIFT2.
%
%   CONVOL is an ARMASA main function.
%
%   See also: CONVOLREV, DECONVOL, CONV.

%Header
%=============================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
switch nargin
case 1
   if isa(x,'struct'), ASAcontrol=x; x=[];
   else, ASAcontrol=[];
   end
   y=[]; shift1=[]; shift2=[];
case 2
   if isa(y,'struct'), ASAcontrol=y; y=[];
   else, ASAcontrol=[];
   end
   shift1=[]; shift2=[];
case 3
   if isa(shift1,'struct'), ASAcontrol=shift1; shift1=[]; shift2=[];
   else, ASAcontrol=[]; shift2=shift1; shift1=y; y=[];
   end
case 4
   if isa(shift2,'struct'), ASAcontrol=shift2; shift2=shift1; shift1=y; y=[];
   else, ASAcontrol=[];
   end
case 5
   if isa(last,'struct'), ASAcontrol=last;
   else, error(ASAerr(39))
   end
otherwise
   error(ASAerr(1,mfilename))
end

if isequal(nargin,1) & ~isempty(ASAcontrol)
      %ASAcontrol is the only input argument
   ASAcontrol.error_chk = 0;
   ASAcontrol.run = 0;
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 6 12 17 20];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 4 15 0 0];

%Checks
%------

if ~any(strcmp(fieldnames(ASAcontrol),'error_chk')) | ASAcontrol.error_chk
      %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(x)
      error(ASAerr(11,'x'))
   end
   if ~isempty(y) & ~isnum(y)
      error(ASAerr(11,'y'))
   end
   if ~isempty(shift1) & (~isnum(shift1) | ...
         ~isintscalar(shift1) | shift1<0)
      error(ASAerr(17,'shift1'))
   end
   if ~isempty(shift2) & (~isnum(shift2) | ...
         ~isintscalar(shift2) | shift2<0)
      error(ASAerr(17,'shift2'))
   end

   %Input argument value checks
   if ~isavector(x)
      error([ASAerr(14) ASAerr(15,'x')])
   else
      l_x = length(x);
      l_y = length(x);
   end
   if ~isempty(y)
      if ~isavector(y)
         error([ASAerr(14) ASAerr(15,'y')])
      end
      l_y = length(y);
   end
   if ~(isempty(shift1) & isempty(shift2)) & ...
         (shift1 > shift2 | shift2 > l_x+l_y-1 | ...
         shift1==0)
      error(ASAerr(20,{'shift1','shift2'}))
   end
end

if ~any(strcmp(fieldnames(ASAcontrol),'version_chk')) | ASAcontrol.version_chk
      %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);
end

if ~any(strcmp(fieldnames(ASAcontrol),'run')) | ASAcontrol.run
      %Run the computational kernel
   ASAcontrol.run = 1;

%Main
%========================================================================

s_x = size(x);
l_x = prod(s_x);
if isempty(y)
   s_y = s_x;
   l_y = l_x;
else
   s_y = size(y);
   l_y = prod(s_y);
end
if l_x>l_y
   l_long = l_x;
   l_short = l_y;
else
   l_long = l_y;
   l_short = l_x;
end
if isempty(shift1)
   shift1 = 1;
   shift2 = l_long+l_short-1;
end
if s_x(1)==1 & s_y(1)==1
   z = zeros(1,shift2-shift1+1);
else
   z = zeros(shift2-shift1+1,1);
end

%Determination of the fastest approach
if l_long<81
   mode=2;
else
   if isempty(y)
      l_fft = nxtppow(l_long+1+max(shift2-l_long,l_long-shift1));
      t_fft = 14.6*l_fft*log(l_fft)+6*l_fft;
      thresh = min(l_long,shift2);
      t_direct = thresh^2-shift1^2+thresh+shift1;
      if shift2>thresh
         t_direct = t_direct+4*l_long*shift2-shift2^2-3*l_long^2;
      end
      if t_fft<t_direct
         mode = 3;
      else
         mode = 1;
      end
   else
      t_direct = 0;
      start = shift1;
      if shift1<l_short
         thresh = min(l_short,shift2);
         t_direct = thresh^2-shift1^2;
         start = l_short;
      end
      if start>=l_short
         thresh = min(l_long,shift2);
         t_direct = t_direct+2*l_short*max(0,(thresh-l_short));
      end
      if shift2>l_long
         t_direct = t_direct+2*l_short*(shift2-l_long)-(shift2-l_long)^2;
      end
      if shift1<l_short
         l_eff = shift2;
      else
         l_eff = shift2-shift1+l_short;
      end
      t_filt = 0.6*l_eff*l_short+1.2*l_short^2;
      l_fft = nxtppow(l_long+l_short);
      t_fft = 14.6*l_fft*log(l_fft)+6*l_fft;
      [mode,index] = sort([t_direct t_filt t_fft]);
      mode = index(1);
   end
end

if mode==1 %Direct approach,
           %minimized number of flops
  if isempty(y)
     y = x;
  end
  if l_x>=l_y
      if s_x(2)>1
         long = x;
      else
         long = x';
      end
      if s_y(2)>1;
         short = y(l_y:-1:1)';
      else
         short = y(l_y:-1:1);
      end
      first = shift1;
      last = shift2;
      k = 1;
      k_increm = 1;
   else
      if s_x(2)>1
         short = x';
      else
         short = x;
      end
      if s_y(2)>1;
         long = y(l_y:-1:1);
      else
         long = y(l_y:-1:1)';
      end
      first = l_x+l_y-shift2;
      last = l_x+l_y-shift1;
      k = last-first+1;
      k_increm = -1;
   end
   start = first;
   run_in = l_short-first;
   if run_in>0
      index_temp = run_in+1;
      i_stop = min(last-first,run_in-1);
      for i = 0:i_stop
         z(k) = long(1:first+i)*short(index_temp-i:l_short);
         k = k+k_increm;
      end
      start = l_short;
   end
   run_center = start-l_short;
   if run_center>=0
      if last<=l_long
         stop = last;
      else
         stop = l_long;
      end
      index_temp = run_center+1;
      i_stop = stop-start;
      for i=0:i_stop
         z(k) = long(index_temp+i:start+i)*short;
         k = k+k_increm;
      end
   end
   run_out = last-l_long;
   if run_out>0
      index_temp = l_long-l_short+1;
      i_start = max(1,first-l_long);
      for i = i_start:run_out
         z(k) = long(index_temp+i:l_long)*short(1:l_short-i);
         k = k+k_increm;
      end
   end
elseif mode==2 %Filter approach,
               %non-optimal in the sense of flops but sometimes
               %faster because 'filter' is a built-in function
   if isempty(y)
      y = x;
   end
   if l_x>l_y
      if s_x(1)==1 & s_y(1)==1
         long = zeros(1,l_x+l_y-1);
         long(1:l_x) = x;
      else
         long = zeros(l_x+l_y-1,1);
         long(1:l_x) = x;
      end
      short = y;
   else
      if s_x(1)==1 & s_y(1)==1
         long = zeros(1,l_x+l_y-1);
         long(1:l_y) = y;
      else
         long = zeros(l_x+l_y-1,1);
         long(1:l_y) = y;
      end
      short = x;
   end
   first = shift1;
   last = shift2;
   if first-l_short>0
      start = first-l_short+1;
      first = l_short;
   else
      start = 1;
   end
   stop = last;
   z = filter(short,1,long(start:stop));
   z = z(first:end);
elseif mode==3 %FFT approach,
               %fast for long sequences
   if ~isempty(y) %cross-convolution
      first = shift1;
      last = shift2;
      run_in = l_short-first;
      if run_in<0
         start = 1-run_in;
         first = l_short;
      else
         start = 1;
      end
      if last<l_long
         stop = last;
      else
         stop = l_long;
      end
      first = first+1;
      last = last-start+2;
      l_long = stop-start+1;
      sig_mtx = zeros(l_fft,2);
      if l_x>l_y
         x = x(stop:-1:start);
         sig_mtx(l_fft-l_long+1:l_fft,1) = x(:);
         sig_mtx(1:l_y,2) = y(:);
      else
         y = y(stop:-1:start);
         sig_mtx(1:l_x,2) = x(:);
         sig_mtx(l_fft-l_long+1:l_fft,1,1) = y(:);
      end
      transf_mtx = fft(sig_mtx);
      z = real(fft(transf_mtx(:,1).*conj(transf_mtx(:,2))));
      if s_x(1)==1 & s_y(1)==1
         z = z(first:last)'/l_fft;
      else
         z = z(first:last)/l_fft;
      end
   else %auto-convolution
      l_sig = l_x;
      first = shift1+1;
      last = shift2+1;
      sig_mtx = zeros(l_fft,2);
      sig_mtx(1:l_x,2) = x(:);
      x = x(l_x:-1:1);
      sig_mtx(l_fft-l_x+1:l_fft,1) = x(:);
      transf_mtx = fft(sig_mtx);
      z = real(fft(transf_mtx(:,1).*conj(transf_mtx(:,2))));
      if s_x(1)==1 & s_y(1)==1
         z = z(first:last)'/l_fft;
      else
         z = z(first:last)/l_fft;
      end
   end
end

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   z = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% [2000 12 4 15  0  0]   W. Wunderink           wwunderink01@freeler.nl
% [2000 12 6 12 17 20]   W. Wunderink           wwunderink01@freeler.nl


text:
function probVal=trivarNormCDF(b,mu,R)
%%TRIVARNORMCDF Evaluate the cumulative density function of the trivariate
%              normal distribution with a specified mean and covariance
%              matrix. This evaluates Pr{x1<b(1), x2<b(2), x3<b(3)) where
%              the random vector is [x1;x2;x3].
%
%INPUTS: b A 2X1 vector [b1;b2] such that b1 is the upper bound of the
%           first variable and b2 is the upper bound on the second
%           variable.
%        mu The 2X1 mean of the distribution. If this is omitted or an
%           empty matrix is passed, then the default of [0;0] is used.
%         R The 2X2 covariance matrix of the distribution. If this is
%           omitted or an empty matrix is passed, then R=eye(2,2) is used.
%
%OUTPUTS: probVal The value of the trivariate normal CDF.
%
%This function implements the second Plankett method that is described in
%[1]. Specifically, Equation 14 is used.
%
%REFERENCES:
%[1] A. Genz, "Numerical computation of rectangular bivariate and
%    trivariate normal and t probabilites," Statistics and Computing, vol.
%    14, pp. 251-260, 2004.
%
%July 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public releas

if(nargin<3||isempty(R))
    R=eye(3,3);
end

if(nargin<2||isempty(mu))
    mu=[0;0;0];
end

%The exteme for one bound.
if(any(b)==-Inf)
    probVal=0;
    return;
end

%If any of the bounds is Inf, then then problem reduces to a 2D integral.
if(b(1)==Inf)
    R=R(2:3,2:3);
    mu=mu(2:3);
    b=b(2:3);
    probVal=bivarNormCDF(b,mu,R);
    return
elseif(b(2)==Inf)
     R=R([1,3],[1,3]);
     mu=mu([1;3]);
     b=b([1;3]);
     probVal=bivarNormCDF(b,mu,R);
     return;
elseif(b(3)==Inf)
    R=R(1:2,1:2);
    mu=mu(1:2);
    b=b(1:2);
    probVal=bivarNormCDF(b,mu,R);
    return
end

%Center the distribution.
b=b-mu;

%Scale the R matrix (and the associated b values) to make the diagonals of
%R all 1.
S=diag(1./sqrt(diag(R)));
b=S*b;
R=S*R*S';

%The correlation values.
pVec=[R(2,1);R(3,1);R(3,2)];

%We want the variables of integration to be permuted to minimize
%max(abs(p21),abs(p31)). This means that out of the above values, p21, p31,
%and p32, we want to change the ordering so that p21 and p31 are the values
%with the smallest magnitudes.
[~,idx]=sort(abs(pVec),'ascend');
pVec=pVec(idx);
b=b(idx);

b1=b(1);
b2=b(2);
b3=b(3);
p21=pVec(1);
p31=pVec(2);
p32=pVec(3);

R2D=[1,  p32;
     p32,1];
%The first term in Equation 14 in 1.
term1=GaussianD.CDF(b1)*bivarNormCDF([b2;b3],[],R2D);
f=@(t)costFun(t,b1,b2,b3,p21,p31,p32);

RelTol=1e-18;
AbsTol=1e-18;
term2=integral1DAdaptive(f,[0;1],[],[],[],RelTol,AbsTol);
probVal=term1+term2;

%In case finite precision limitations make the value invalid.
probVal=min(max(probVal,0),1);

end

function val=costFun(t,b1,b2,b3,p21,p31,p32)

t2=t.^2;
denomTerm2=1-p31^2*t2-p21^2*t2-p32^2+2*t2*p31*p21*p32;
denomTerm3=1-p21^2*t2-p31^2*t2-p32^2+2*t2*p31*p21*p32;
u2Hat=(b2*(1-p31^2*t2)-b1*t*(p21-p31*p32)-b3*(p32-p31*p21*t2))./sqrt((1-p31^2.*t2).*denomTerm2);
u3Hat=(b3*(1-p21^2*t2)-b1*t*(p31-p21*p32)-b2*(p32-p31*p21*t2))./sqrt((1-p21^2.*t2).*denomTerm3);

r=p31*t;
f2=(b1^2+b3^2-2*r*b1*b3)./(1-r.^2);
f2(isnan(f2))=Inf;
r=p21*t;
f3=(b1^2+b2^2-2*r*b1*b2)./(1-r.^2);
f3(isnan(f3))=Inf;
%The above NaN stuff should make it work properly if, for example, b1 and
%b2 are near realmax or realmin.

Phiu3Hat=GaussianD.CDF(u3Hat);
Phiu2Hat=GaussianD.CDF(u2Hat);

val=Phiu3Hat.*p21.*exp(-f3/2)./sqrt(1-p21^2.*t2)+Phiu2Hat.*p31.*exp(-f2/2)./sqrt(1-p31.^2.*t2);

end

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


text:
function anss = rombquad(f,a,b,m)
%ROMBQUAD(f,a,b,m) evaluates the integral of f from a to b.
% f is an inline function or function handle.
% a is the lower limit of integration.
% b is the upper limit of integration.
% m - is the exponent in the relative error.  For example m = 10 implies an
% estimated relative error of 1e-10.
%
% Example:
%
%           f = @(x) x.^2; % Vectoized function handle.
%           rombquad(f,0,2); % The integral from 0 to 2.
%
% See also quad, quadl, quadv and gaussquad (on the FEX)
%
% Author:  Matt Fig
% Contact:  popkenai@yahoo.com

if nargin<3
   error ('Not enough arguments. See help.')
elseif nargin<4 || floor(m)~=m || m <1
    m = 10;
end

if isinf(a) || isinf(b)
   error('Infinite limits not allowed.')
end

h = 2.^((1:m)-1)*(b-a)/2^(m-1);             % These are the intervals used.
k1 = 2.^((m-2):-1:-1)*2+1;                  % Index into the intervals.
f1 = feval(f,a:h(1):b);                     % Function evaluations.
R = zeros(1,m);                             % Pre-allocation.
% Define the starting vector:
for ii = 1:m
	R(ii) = 0.5*h(ii)*(f1(1)+2*...
            sum(f1(k1(end-ii+1):k1(end-ii+1)-1:end-1)) + f1(end));
end
% Interpolations:
for jj = 2:m
    jpower = (4^(jj-1)-1);
    for kk = 1:(m-jj+1)
        R(kk) = R(kk)+(R(kk)-R(kk+1))/jpower;
    end
end
anss = R(1);