text:
function []=anim8_DIC_3D_pairs_faceMeasure_direction(DIC3DAllPairsResults,faceMeasureString,varargin)
%% function for plotting 3D-DIC results of face measures as color + direction as arrows in STEP3.
% plotting 3D surfaces from camera pairs, animation changing
% with time, and the faces colored according to faceMeasureString
% this function is called in plotMultiDICPairResults
%
% Options:
% anim8_DIC_3D_pairs_faceMeasure(DIC_3Dallpairs_results,faceMeasureString)
% anim8_DIC_3D_pairs_faceMeasure(DIC_3Dallpairs_results,faceMeasureString,optStruct)
%
% Inputs:
% * DIC3DAllPairsResults
% * faceMeasureString: can be any of the following:
%   'Epc1','Epc2','epc1','epc2','Lamda1','Lamda2', or a combination of two  in a cell array to plot side by side
% * optStruct: optional structure for plotting options which may include any of the following fields:
%   - smoothLogic: logical variable for smoothing (true)/not smoothing (false) the face measure
%   - FaceAlpha: transparacy of the faces (scalar between 0 and 1, where zero is transparent and 1 is opaque)
%   - colorBarLimits: a 2x1 scalar vector for the colobar limits. if not set, it's automatic
%   - dataLimits: a 2x1 scalar vector for the data limits of the face measure. if a face measure is outside these limits, it is set to NaN. if not set no face is set to NaN
%   - colorMap
%   - zDirection: 1 for z up and -1 for z down
%   - lineColor: line color for the mesh. can be for example 'b','k','none',etc...
%   - supTitleString=faceMeasureString;

%% Assign plot options
Narg=numel(varargin);
switch Narg
    case 1
        optStruct=varargin{1};
    case 0
        optStruct=struct;
    otherwise
        ('wrong number of input arguments');
end

% complete the struct fields
if ~isfield(optStruct,'smoothLogic')
    optStruct.smoothLogic=0;
end
if ~isfield(optStruct,'FaceAlpha')
    optStruct.FaceAlpha=1;
end
if ~isfield(optStruct,'dataLimits')
    optStruct.dataLimits=[-inf inf];
end
if ~isfield(optStruct,'zDirection') % 1 or -1
    optStruct.zDirection=1;
end
if ~isfield(optStruct,'lineColor') % 'none' or 'k'
    optStruct.lineColor='none';
end
if ~isfield(optStruct,'maxCorrCoeff')
    optStruct.maxCorrCoeff=[];
end

%%
nFrames=numel(DIC3DAllPairsResults{1}.Points3D);

quiverScaleFactor=20;

DIC_3Djoined_results=joinPairs(DIC3DAllPairsResults);
[xl,yl,zl]=axesLimits(DIC_3Djoined_results.Points3D);

%% Assign the right face measure into FC

if iscell(faceMeasureString)
    faceMeasureCell=faceMeasureString;
    nStrains=numel(faceMeasureString);
    switch nStrains
        case 1
            FC=cell(nFrames,nStrains);
            D=cell(nFrames,nStrains);
            Ds=cell(nFrames,nStrains);
            Vc=cell(nFrames,nStrains);
        case 2
            FC=cell(nFrames,nStrains);
            D=cell(nFrames,nStrains);
            Ds=cell(nFrames,nStrains);
            Vc=cell(nFrames,nStrains);
        otherwise
            error('wrong number of cells in faceMeasureString cell array (must be 1 or 2)');
    end
elseif ischar(faceMeasureString)
    nStrains=1;
    faceMeasureCell=cell(1);
    faceMeasureCell{1}=faceMeasureString;
    FC=cell(nFrames,1);
    D=cell(nFrames,1);
    Ds=cell(nFrames,1);
    Vc=cell(nFrames,1);
else
    error('wrong face measure (second input variable)');
end

if ~isfield(optStruct,'supTitleString') % 'none' or 'k'
    optStruct.supTitleString=faceMeasureCell;
end

directionStringCell=cell(nStrains,1);
for is=1:nStrains
    switch faceMeasureCell{is}
        case 'Epc1'
            directionStringCell{is}='Epc1vecCur';
        case 'Epc2'
            directionStringCell{is}='Epc2vecCur';
        case 'epc1'
            directionStringCell{is}='epc1vec';
        case 'epc2'
            directionStringCell{is}='epc2vec';
        case 'Lamda1'
            directionStringCell{is}='Epc1vecCur';
        case 'Lamda2'
            directionStringCell{is}='Epc2vecCur';
        otherwise
            error('unexpected face measure string. plots not created');
    end
    for it=1:nFrames
        FC{it,is}=DIC_3Djoined_results.Deform.(faceMeasureCell{is}){it}; % face color (strain)
        D{it,is}=DIC_3Djoined_results.Deform.(directionStringCell{is}){it}; % direction (unit vector)
        Vc{it,is}=DIC_3Djoined_results.FaceCentroids{it};
        if ~isempty(optStruct.maxCorrCoeff)
            corrNow=DIC_3Djoined_results.FaceCorrComb{it};
            FC{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
            D{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
            Vc{it,is}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
        end
        switch faceMeasureCell{is}
            case {'Epc1','Epc2','epc1','epc2'}
                Ds{it,is}=quiverScaleFactor*FC{it,is}.*D{it,is}; % direction with magnitude (scaled vector)
            case {'Lamda1','Lamda2'}
                Ds{it,is}=quiverScaleFactor*(FC{it,is}-1).*D{it,is}; % direction with magnitude (scaled vector)
        end

    end
end
FCmat = cell2mat(FC);
if ~isfield(optStruct,'colorBarLimits')
    switch faceMeasureCell{is}
        case {'Epc1','Epc2','epc1','epc2'}
            optStruct.colorBarLimits=[-prctile(abs(FCmat(:)),100) prctile(abs(FCmat(:)),100)];
        case {'Lamda1','Lamda2'}
            optStruct.colorBarLimits=[1-prctile(abs(FCmat(:)-1),100) 1+prctile(abs(FCmat(:)-1),100)];
    end
end
colorBarLogic=1;
if ~isfield(optStruct,'colorMap')
    optStruct.colorMap=.8*coldwarm;
end

%% Plot

animStruct=struct;

hf=cFigure;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

for is=1:nStrains
    subplot(1,nStrains,is);

    axisGeom;
    ax=gca;
    ax.CameraUpVector=[0 0 optStruct.zDirection];
    colormap(optStruct.colorMap);
    if colorBarLogic
        colorbar;
        caxis(optStruct.colorBarLimits);
    end

    title(optStruct.supTitleString{is});
    axis off
    % camlight headlight

    it=1;
    Fnow=DIC_3Djoined_results.Faces;
    Pnow=DIC_3Djoined_results.Points3D{it};
    Vnow=Vc{it,is};
    FCnow=FC{it,is};
    Dnow=Ds{it,is};
    if optStruct.smoothLogic
        [FCnow]=triSmoothFaceMeasure(FCnow,Fnow,Pnow,[],[]);
    end
    FCnow(FCnow<optStruct.dataLimits(1))=NaN;
    FCnow(FCnow>optStruct.dataLimits(2))=NaN;

    hp(is)=gpatch(Fnow,Pnow,FCnow,optStruct.lineColor,optStruct.FaceAlpha); hold on
    hq(is)=quiver3(Vnow(:,1),Vnow(:,2),Vnow(:,3),Dnow(:,1),Dnow(:,2),Dnow(:,3),0,'Color',.2*[1 1 1],'ShowArrowHead','off','AutoScale','on'); hold on;
%
    h_ax=gca;
    h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

end

%% fill in the animstruct

animStruct.Time=1:nFrames;
animStruct.Handles=cell(1,nFrames);
animStruct.Props=cell(1,nFrames);
animStruct.Set=cell(1,nFrames);

for it=1:nFrames
    animStruct.Handles{it}=[];
    animStruct.Props{it}=cell(1,8*nStrains);
    animStruct.Set{it}=cell(1,8*nStrains);
    Pnow=DIC_3Djoined_results.Points3D{it};

    for is=1:nStrains
        FCnow=FC{it,is};
        if optStruct.smoothLogic
            [FCnow]=triSmoothFaceMeasure(FCnow,Fnow,Pnow,[],[]);
        end
        FCnow(FCnow<optStruct.dataLimits(1))=NaN;
        FCnow(FCnow>optStruct.dataLimits(2))=NaN;

        Fnow=DIC_3Djoined_results.Faces;
        Pnow=DIC_3Djoined_results.Points3D{it};
        Vnow=Vc{it,is};
        Dnow=Ds{it,is};

        animStruct.Handles{it}=[animStruct.Handles{it} hp(is) hp(is) hq(is) hq(is) hq(is) hq(is) hq(is) hq(is)]; %Handles of objects to animate (add one every pair)

        animStruct.Props{it}{1+8*(is-1)}='CData';
        animStruct.Props{it}{2+8*(is-1)}='Vertices'; %Properties of objects to animate
        animStruct.Props{it}{3+8*(is-1)}='XData'; %Properties of objects to animate
        animStruct.Props{it}{4+8*(is-1)}='YData'; %Properties of objects to animate
        animStruct.Props{it}{5+8*(is-1)}='ZData'; %Properties of objects to animate
        animStruct.Props{it}{6+8*(is-1)}='UData'; %Properties of objects to animate
        animStruct.Props{it}{7+8*(is-1)}='VData'; %Properties of objects to animate
        animStruct.Props{it}{8+8*(is-1)}='WData'; %Properties of objects to animate

        animStruct.Set{it}{1+8*(is-1)}=FCnow;
        animStruct.Set{it}{2+8*(is-1)}=Pnow; %Property values for to set in order to animate
        animStruct.Set{it}{3+8*(is-1)}=Vnow(:,1); %Property values for to set in order to animate
        animStruct.Set{it}{4+8*(is-1)}=Vnow(:,2); %Property values for to set in order to animate
        animStruct.Set{it}{5+8*(is-1)}=Vnow(:,3); %Property values for to set in order to animate
        animStruct.Set{it}{6+8*(is-1)}=Dnow(:,1); %Property values for to set in order to animate
        animStruct.Set{it}{7+8*(is-1)}=Dnow(:,2); %Property values for to set in order to animate
        animStruct.Set{it}{8+8*(is-1)}=Dnow(:,3); %Property values for to set in order to animate

        h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

    end
end

anim8(hf,animStruct);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function tree(objectorder,heights)

%TREE creates a tree in which the leaves represent
%   objects.  The vertical coordinate of the junction
%   of two branches is the dissimilarity between the
%   corresponding clusters (maximal 30 objects allowed).
%
% The algorithm is fully described in:
%   Kaufman, L. and Rousseeuw, P.J. (1990),
%   "Finding groups in data: An introduction to cluster analysis",
%   Wiley-Interscience: New York (Series in Applied Probability and
%   Statistics), ISBN 0-471-87876-6.
%
% Required input arguments:
%   objectorder :  order of objects
%   heights     : diameter of cluster before dividing it
%                 (=length of banner)
%
% I/O:
%   tree(objectorder,heights)
%
% Example (subtracted from the referenced book)
%   load agricul.mat
%   result = diana(agricul,[4 4],0,0,1);
%   tree(result.objectorder,result.heights)
%
% The output of TREE is a figure containing the
%   agglomerative (agnes) or divise (diana) tree.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Wai Yan Kong (May 2006)
% Last Revision: 28/09/2006

clf reset
whitebg([1 1 1]);

if (nargin<2)
    error('Two input arguments required')
elseif (nargin>2)
    error('Too many input arguments')
end

if(size(objectorder,2)~=size(heights,2)+1)
    error('Missing values in objectorder or heights')
end

Heights=H(heights);

number=size(objectorder,2);

if (number>30)
    error('Only 30 objects allowed')
end

%midden=[];
middle(1)=0;
Maxi=0;
Mini=0;
PrevMID=middle(1);
lengt=[]; %length
high=[]; %height
indices=[];
L=[];

[maxim,index]=max(Heights);
if(1<=index-1)
    [Prevmax,Previndex]=max(Heights(1:(index-1)));
else
    Prevmax=maxim;
    Previndex=index;
end
if(index+1<=number-1)
    [Postmax,Postindex]=max(Heights((index+1):(number-1)));
    Postindex=Postindex+index+1-1;
else
    Postmax=maxim;
    Postindex=index;
end

if(Postindex+1<=number-1)
    [Post2max,Post2index]=max(Heights(Postindex+1:number-1));
    Post2index=Post2index+Postindex+1-1;
else
    Post2max=Postmax;
    Post2index=Postindex;
end

if(Postindex+1<=Post2index-1)
    [Betweenmax,Betweenindex]=max(Heights(index+1:Postindex-1));
    Betweenindex=Betweenindex+index+1-1;
else
    Betweenmax=Postmax;
    Betweenindex=Postindex;
end

L=cat(2,L,maxim);

high=cat(2,high,[Post2max,Postmax,Betweenmax,maxim,Prevmax]);
indices=cat(2,indices,[Post2index,Postindex,Betweenindex,index,Previndex]);

M=0;
extra=number/2;
lengt(1)=Postindex-Previndex+extra+2;
PrevLEN=lengt(1);
rectangle('Position',[middle(1)-(lengt(1)/2),maxim,lengt(1),0.0001]);

NbanFirst=1;
k=1;
branch=0;
over=0;
Element=0;
Special=0;
Spec=0;
Sp=0;
S=0;
ww=0;
if(maxim<20)
    extrawaystick=0.8;
    extrawaytext=1.5;
else
    extrawaytext=0;
    extrawaystick=0;
end

right=0;
direct11=0;
direct111=0;
righttree=0;
if(index==1)
    direct1=1;
else
    direct1=0;
end

Last=size(lengt,2);
LastM=size(middle,2);
LastL=size(lengt,2);

while(k<=number)
    w=0;
    LastH=size(high,2);
    if((high(LastH)~=high(LastH-1)) & Element==0)
        while(branch==0)
            S=0;
            w=w+1;
            ww=ww+1;

            if(over==0)
                middle=cat(2,middle,middle(LastM)-(lengt(LastL)/2));
            end

            if(righttree==1)
                middle=cat(2,middle,middle(LastM)-(lengt(LastL)/2));
                righttree=0;
            end

            if(Sp==1)
                LastM=size(middle,2);
                LastL=size(lengt,2);
                middle(LastM)=middle(LastM)+lengt(LastL)/2;
                Sp=0;
                S=1;
            end
            LastM=size(middle,2);

            Output=maxim-Prevmax;
            rectangle('Position',[middle(LastM),Prevmax,0.0001,Output]);

            if(Betweenmax==Postmax | over==0)
                Post2max=Postmax;
                Post2index=Postindex;
            end

            if((Betweenmax==Postmax | over==0) & (right~=1 | over==0))
                Postmax=maxim;
                Postindex=index;
            end

            maxim=Prevmax;
            index=Previndex;

            if(NbanFirst<=(index-1))
                [Prevmax,Previndex]=max(Heights(NbanFirst:(index-1)));
                Previndex=Previndex+NbanFirst-1;
            else
                Prevmax=maxim;
                Previndex=index;
            end

            if(index+1<=Postindex-1)
                [Betweenmax,Betweenindex]=max(Heights(index+1:Postindex-1));
                Betweenindex=Betweenindex+index+1-1;
            else
                Betweenmax=Postmax;
                Betweenindex=Postindex;
            end

            if(right==1 & over~=0)
                high=cat(2,high,[Prevmax Prevmax]);
                indices=cat(2,indices,[Previndex Previndex]);
            else
                high=cat(2,high,Prevmax);
                indices=cat(2,indices,Previndex);
            end

            if(Postindex-Previndex<=0)
                lengt=cat(2,lengt,1);
            else
                if(Postindex-Previndex-direct111>0)
                    if(ww==1)
                        lengt=cat(2,lengt,Postindex-Previndex-direct111+2);
                    else
                        lengt=cat(2,lengt,Postindex-Previndex-direct111);
                    end
                else
                    if(ww==1)
                        lengt=cat(2,lengt,Postindex-Previndex+2);
                    else
                        lengt=cat(2,lengt,Postindex-Previndex);
                    end
                end
            end
            L=cat(2,L,maxim);

            if(direct111~=0)
                direct111=0;
            end

            LastH=size(high,2);
            LastL=size(lengt,2);
            LastM=size(middle,2);

            rectangle('Position',[middle(LastM)-(lengt(LastL)/2),maxim,lengt(LastL),0.0001]);

            if(over~=0)
                over=0;
            end

            if(right~=0)
                right=0;
            end

            if(NbanFirst==index)
                branch=1;
            end
        end

        rectangle('Position',[middle(LastM)-(lengt(LastL)/2),maxim-1+extrawaystick,0.0001,1-extrawaystick]);
        if(objectorder(k)>10)

            text(middle(LastM)-(lengt(LastL)/2)-0.25,maxim-2.0+extrawaytext,num2str(double(objectorder(k))));
            T=middle(LastM)-(lengt(LastL)/2-0.25);
        else

            text(middle(LastM)-(lengt(LastL)/2)-0.1,maxim-2.0+extrawaytext,num2str(double(objectorder(k))));
            T=middle(LastM)-(lengt(LastL)/2-0.1);
        end

        Maxi=max(Maxi,T);
        Mini=min(Mini,T);

        if(Betweenmax~=Postmax)
            high=high(1:LastH-2);
            indices=indices(1:LastH-2);

            Prevmax=Betweenmax;
            Previndex=Betweenindex;

            high=cat(2,high,Prevmax);
            indices=cat(2,indices,Previndex);
            LastH=size(high,2);

            LastL=size(lengt,2);
            middle(LastM)=middle(LastM)+lengt(LastL)/2;
            over=1;
        end
        direct1=0;
        Special=0;
        Spec=0;
    else
        if(Special==1)
            Special=0;
            Spec=1;
        end

        if(Sp==1)
            Sp=0;
        end

        if(S==1)
            S=2;
        end

        direct111=0;
        LastH=size(high,2);
        LastM=size(middle,2);
        LastL=size(lengt,2);

        if(direct1==0)

            rectangle('Position',[middle(LastM)+(lengt(LastL)/2),high(LastH)-1+extrawaystick,0.0001,1-extrawaystick]);

            text(middle(LastM)+lengt(LastL)/2-0.1,high(LastH)-2.0+extrawaytext,num2str(double(objectorder(k))));

            T=middle(LastM)+lengt(LastL)/2-0.1;
            Maxi=max(Maxi,T);
            Mini=min(Mini,T);
        else

            rectangle('Position',[middle(LastM)-(lengt(LastL)/2),high(LastH)-1+extrawaystick,0.0001,1-extrawaystick]);

            text(middle(LastM)-lengt(LastL)/2-0.1,high(LastH)-2.0+extrawaytext,num2str(double(objectorder(k))));

            T=middle(LastM)-lengt(LastL)/2-0.1;
            Maxi=max(Maxi,T);
            Mini=min(Mini,T);
            direct1=0;
            direct11=1;
        end

        if(high(LastH-1)==high(LastH))
            high=high(1:LastH-2);
            indices=indices(1:LastH-2);
        else
            high=high(1:LastH-1);
            indices=indices(1:LastH-1);
        end
        LastH=size(high,2);

        maxim=high(LastH);
        index=indices(LastH);

        if(LastH-1>=1)
            Postmax=high(LastH-1);
            Postindex=indices(LastH-1);
        end

        if(LastH-2>=1)
            Post2max=high(LastH-2);
            Post2index=indices(LastH-2);
        end

        if(index+1<=Postindex-1)
            [Betweenmax,Betweenindex]=max(Heights(index+1:Postindex-1));
            Betweenindex=Betweenindex+index+1-1;
        else
            Betweenmax=Postmax;
            Betweenindex=Postindex;
        end

        if(Betweenmax~=Postmax)
            Prevmax=Betweenmax;
            Previndex=Betweenindex;
        else
            Prevmax=maxim;
            Previndex=index;
        end

        if(high(LastH)>high(LastH-1))
            righttree=1;
        end

        high(LastH)=Prevmax;
        indices(LastH)=Previndex;

        LastM=size(middle,2);
        LastL=size(lengt,2);
        Last=size(L,2);

        if(LastM-1>=1)
            middle=middle(1:LastM-1);
        end

        if(LastL-1>=1)
            lengt=lengt(1:LastL-1);
        end

        if(Last-1>=1)
            L=L(1:Last-1);
        end

        LastM=size(middle,2);
        LastL=size(lengt,2);
        Last=size(L,2);

        if(Last>0)
            while(maxim>L(Last))
                lengt=lengt(1:Last-1);
                L=L(1:Last-1);
                Last=size(L,2);
            end
        end
        LastL=size(lengt,2);

        middle(LastM)=middle(LastM)+(lengt(LastL)/2);

        if(Betweenmax~=Postmax)
            Element=0;
        else
            Element=1;
            middle(LastM)=middle(LastM)-(lengt(LastL)/2);
        end

        if(righttree==1 | direct11==1)
            if(Spec==1)
                Sp=1;
            end

            PrevMAX=L(1);
            [maxim,index]=max(Heights(NbanFirst+1:number-1));
            index=index+NbanFirst+1-1;

            if(NbanFirst+1<=index-1)
                [Prevmax,Previndex]=max(Heights(NbanFirst+1:index-1));
                Previndex=Previndex+NbanFirst+1-1;
            else
                Prevmax=maxim;
                Previndex=index;
            end

            if(index+1<=number-1)
                [Postmax,Postindex]=max(Heights(index+1:number-1));
                Postindex=Postindex+index+1-1;
            else
                Postmax=maxim;
                Postindex=index;
            end

            if(index+1<=Postindex-1)
                [Betweenmax,Betweenindex]=max(Heights(index+1:Postindex-1));
                Betweenindex=Betweenindex+index+1-1;
            else
                Betweenmax=Postmax;
                Betweenindex=Postindex;
            end

            if(Postindex+1<=number-1)
                [Post2max,Post2index]=max(Heights(Postindex+1:number-1));
                Post2index=Post2index+Postindex+1-1;
            else
                Post2max=Postmax;
                Post2index=Postindex;
            end

            rectangle('Position',[PrevMID+(PrevLEN/2),maxim,0.0001,PrevMAX-maxim]);

            LastM=size(middle,2);
            LastMID=PrevMID+(PrevLEN/2);
            PrevMID=LastMID;
            middle=[];
            middle(1)=LastMID;
            LastM=size(middle,2);

            high=[];
            indices=[];
            high=cat(2,high,[Post2max,Postmax,Betweenmax,maxim,Prevmax]);

            indices=cat(2,indices,[Post2index,Postindex,Betweenindex,index,Previndex]);
            LastH=size(high,2);

            L=[];
            L=cat(2,L,maxim);
            Last=size(L,2);
            lengt=[];
            if(Postindex>Previndex)
                if(M<extra)
                    M=M+1;
                end
                lengt(1)=Postindex-Previndex+(extra-M);

                rectangle('Position',[LastMID-lengt(1)/2,maxim,lengt(1),0.0001]);
                direct1=1;
            elseif(Postindex<=Previndex & Betweenmax==Postmax)
                lengt(1)=1;
                rectangle('Position',[LastMID-1/2,maxim,1,0.0001]);
                %%%%%%%

                rectangle('Position',[LastMID-1/2,maxim-1+extrawaystick,0.0001,1-extrawaystick]);
                if(objectorder(k+1)>10)

                    text(LastMID-1/2-0.25,maxim-2.0+extrawaytext,num2str(double(objectorder(k+1))));
                    Mini=min(Mini,LastMID-1/2-0.25);
                else

                    text(LastMID-1/2-0.1,maxim-2.0+extrawaytext,num2str(double(objectorder(k+1))));
                    Mini=min(Mini,LastMID-1/2-0.1);
                end

                rectangle('Position',[LastMID+1/2,maxim-1+extrawaystick,0.0001,1-extrawaystick]);

                text(LastMID+1/2-0.1,maxim-2.0+extrawaytext,num2str(double(objectorder(k+2))));
                Maxi=max(Maxi,LastMID+1/2-0.1);

                k=k+2;
            end
            PrevLEN=lengt(1);
            LastL=size(lengt,2);

            if(Element==1)
                Element=0;
            end

            if(direct11==1)
                if(NbanFirst+1<index)
                    middle=cat(2,middle,middle(LastM)-lengt(LastL)/2);
                    LastM=size(middle,2);
                    direct111=1;
                    over=1;
                end
                direct11=0;
            end

            Special=1;
            PrevMAX=maxim;
        end
        if(S==2 && k+1==number)

            rectangle('Position',[LastMID+lengt(1)/2,PrevMAX-1+extrawaystick,0.0001,1-extrawaystick]);

            text(LastMID+lengt(1)/2-0.1,PrevMAX-2.0+extrawaytext,num2str(double(objectorder(k+1))));
            Maxi=max(Maxi,LastMID+lengt(1)/2-0.1);
            Mini=min(Mini,LastMID+lengt(1)/2-0.1);
            k=k+1;
        end
        over=1;
    end
    branch=0;

    k=k+1;
    NbanFirst=NbanFirst+1;

end
axis([Mini-1,Maxi+1,min(Heights)-5,max(Heights)+1]);
XT=[];
set(gca,'XTick',XT);
set(gca,'XTickLabel',[]);

%---
function res=H(vector)

lengt=size(vector,2);
Nvector=vector;
i=1;
for i=1:lengt
    for j=i+1:lengt
        if (vector(i)==vector(j))
            Nvector(j)=vector(j)+i*0.0001;
            i=i+1;
        end
    end
end
res=Nvector;

author: marianux

repo: ecg-kit

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


text:
function omega_en_n = transport_rate(lat, Vn, Ve, h)
% transport_rate: calculates the transport rate in the navigation frame.
%
% INPUT
%	lat, 1x1 latitude (rad).
%	Vn, 1x1 North velocity (m/s).
%   Ve, 1x1 East velocity (m/s).
%   h, altitude (m).
%
% OUTPUT
%	omega_en_n, 3x3 skew-symmetric transport rate matrix (rad/s).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   Paul D. Groves. Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems. Second Edition.
% Eq. 5.44, page 177.
%
% Version: 002
% Date:    2022/01/25
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

h = abs(h);

[RM,RN] = radius(lat);

om_en_n(1,1) =   Ve / (RN + h);              % North
om_en_n(2,1) = -(Vn / (RM + h));             % East
om_en_n(3,1) = -(Ve * tan(lat) / (RN + h));  % Down

omega_en_n = skewm(om_en_n);

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function D = all_pairs_distances(V,U)
  % ALL_PAIRS_DISTANCES compute distances between each point i in V and point j
  % in U
  %
  % This is obsolete: use pdist2 instead
  %
  % D = all_pairs_distances(V,U)
  %
  % Inputs:
  %   V  #V by dim list of points
  %   U  #U by dim list of points
  % Outputs:
  %   D  #V by #U matrix of distances, where D(i,j) gives the distance between
  %     V(i,:) and U(j,:)
  %

  warning('obsolete. Call `pdist2` directly instead');

  assert(size(V,2) == size(U,2));

  %% Super-slow for loops method
  %% Elapsed time is 128.969203 seconds.
  %D = zeros(size(V,1),size(U,1));
  %for i = 1:size(V,1)
  %  for j = 1:size(U,1)
  %    D(i,j) = sqrt(sum((V(i,:)-U(j,:)).^2,2));
  %  end
  %end

  %% Slow, especially if hitting memory paging
  %% Elapsed time is 2.519262 seconds.
  %D = ...
  %  squeeze(sqrt(sum((repmat(V,[1 1 size(U,1)])- ...
  %  permute(repmat(U,[1 1 size(V,1)]),[3 2 1])).^2,2)));

  %% Faster single for loop, single repmat method, handles memory paging better
  %% Elapsed time is 1.160419 seconds.
  %D = zeros(size(V,1),size(U,1));
  %for i = 1:size(V,1)
  %  D(i,:) = ...
  %    sqrt(sum((repmat(V(i,:),[size(U,1) 1]) - U).^2,2))';
  %end

  if size(V,1) > 1000 || size(U,1) > 1000
    % Fast single for loop, single bsxfun method, but matlab can handle both for
    % loops in bsxfun. This seems to be a nice balance between hitting a memory
    % wall for big input and still being fast for medium and slow input
    % Elapsed time is 0.863138 seconds.
    D = zeros(size(V,1),size(U,1));
    for i = 1:size(V,1)
      D(i,:) = ...
        sqrt(sum(bsxfun(@minus,V(i,:),U).^2,2));
    end
  else
    % Fastest
    % Elapsed time is 0.653082 seconds.
    D = permute(sqrt(sum(bsxfun(@minus,V,permute(U,[3 2 1])).^2,2)),[1 3 2]);
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [g,L,info] = nonsepgabpars_from_window(g,a,M,lt,L,callfun)
%NONSEPGABPARS_FROM_WINDOW  Compute g and L from window
%   Usage: [g,g.info,L] = gabpars_from_window(f,g,a,M,lt,L);
%
%   Use this function if you know a window and a lattice
%   for the NONSEPDGT. The function will calculate a transform length L and
%   evaluate the window g into numerical form.
%
%   If the transform length is unknown (as it usually is unless explicitly
%   specified by the user), set L to be [] in the input to this function.

if nargin<6
  stacknames=dbstack;
  callfun=stacknames(2).name;
end;

if isempty(L)
  if isnumeric(g)
    L=length(g);
  else
    L=dgtlength(1,a,M,lt);
  end;
else
  Lcheck=dgtlength(L,a,M,lt);
  if Lcheck~=L
    error('%s: Invalid transform size L',upper(mfilename));
  end;
end;

[g,info] = comp_window(g,a,M,L,lt,'NONSEPGABDUAL');

if (info.isfir)
  if info.istight
    g=g/sqrt(2);
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


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

%%
clear; close all; clc;

%% Syntax
% |[Vg]=sweepCurveSmooth(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function X = plus(X,Y)
%PLUS Merges two LMI objects to one LMI

if isa(X,'constraint')
    X = lmi(X);
elseif isa(X,'sdpvar')
    X = lmi(X);
end

if isa(Y,'constraint')
    Y = lmi(Y);
elseif isa(Y,'sdpvar')
    Y = lmi(Y);
end

% Support set+[]
if isempty(X)
    X = Y;
    return
elseif isempty(Y)
    return
end

if ~((isa(X,'lmi')) & (isa(Y,'lmi')))
    error('Both arguments must be constraints')
end

nX = length(X.LMIid);
nY = length(Y.LMIid);
if nX==0
    X = Y;
    return
end
if nY == 0
    return;
end

xBlock = isa(X.clauses{1},'cell');
yBlock = isa(Y.clauses{1},'cell');
if yBlock && xBlock && length(X.clauses{1})>1  && length(Y.clauses{1})>1
    % Both objects are long blocks of constraints. Join these on high level
    for i = 1:length(Y.clauses)
        X.clauses{end+1} = Y.clauses{i};
    end
elseif ~xBlock && ~yBlock
    % Both have been flattened

    % Maybe we should build cells of cells
    if length(X.clauses)> 99
        temp = X.clauses;
        X.clauses = [];
        X.clauses{1} = temp;
        X.clauses{2} = {};
        jx = length(X.clauses);
        for i = 1:length(Y.clauses)
            X.clauses{jx}{end+1} = Y.clauses{i};
        end
    else
        for i = 1:length(Y.clauses)
            X.clauses{end+1} = Y.clauses{i};
        end
    end
else
    % This is the standard case. X has been populated with a bunch of
    % constraints (growing list) while Y is yet another constraint to be
    % added to that list.
    if ~xBlock
        % Lift to a single block
        temp = X.clauses;
        X.clauses = [];
        X.clauses{1} = temp;
    end
    % New block in X? This is where performance comes from. Handling cells
    % of cells is way faster in MATLAB, than a long cell (quadratic running
    % time, appears to be nasty copying going on when adding new elements)
    if length(X.clauses{end}) >= 100
        X.clauses{end+1} = {};
    end
    % Special case for performance
    if isa(Y.clauses{1},'cell') && length(Y.clauses)==1
        j = length(X.clauses);
        for i = 1:length(Y.clauses{1})
            X.clauses{j}{end+1} = Y.clauses{1}{i};
        end
    else
        Y = flatten(Y);
        j = length(X.clauses);
        for i = 1:length(Y.clauses)
            X.clauses{j}{end+1} = Y.clauses{i};
        end
    end
end
aux = X.LMIid;
X.LMIid = [X.LMIid Y.LMIid];

% VERY FAST UNIQUE BECAUSE THIS IS CALLED A LOT OF TIMES....
if ~(max(aux) < min(Y.LMIid))
    i = sort(X.LMIid);
    i = i(diff([i NaN])~=0);
    if length(i)<nX+nY
        % Flatten first. This typically doesn't happen, so we accept this.
        X.clauses = [X.clauses{:}];
        [i,j] = unique(X.LMIid);
        X = subsref(X,struct('type','()','subs',{{j}}));
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function AIC=aic(nParameters, nObservations, Residuals)
%Computes Akaike information criterion for a model
%
%
%
% Implements http://en.wikipedia.org/wiki/Akaike_information_criterion

% HISTORY:
% ER wrote it 01/2010

AIC =2*(nParameters)+nObservations*(log(sum(Residuals.^2)/nObservations));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%This program decomposes an image into blocks.

%Parameters
% inImg         -   Input Gray Image
% blkSize       -   Window size for block processing
% out           -   Output 4 dimensional matrix with blocks.
%
%Author : Athi Narayanan S
%Student, M.E, EST,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%s_athi1983@yahoo.co.in
%http://sites.google.com/site/athisnarayanan/

function out = MatDec(inImg, blkSize)

[m,n]=size(inImg);

r3=m/blkSize;
c3=n/blkSize;
q4=0;
q1=0;

for i=1:r3
    for j=1:c3
        for s=1:blkSize
            for k=1:blkSize
                p3=s+q4;
                q2=k+q1;
                out(s,k,i,j)=inImg(p3,q2);
            end
        end
        q1=q1+blkSize;
    end
    q4=q4+blkSize;q1=0;
end

author: Sable

repo: mcbench-benchmarks


text:
function area = polygon_area_2d_2 ( n, v )

%*****************************************************************************80
%
%% POLYGON_AREA_2D_2 computes the area of a polygon in 2D.
%
%  Discussion:
%
%    The area is the sum of the areas of the triangles formed by
%    node N with consecutive pairs of nodes.
%
%    If the vertices are given in counterclockwise order, the area
%    will be positive.  If the vertices are given in clockwise order,
%    the area will be negative.
%
%    Thanks to Martin Pineault for noticing that an earlier version
%    of this routine would not correctly compute the area of a nonconvex
%    polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(2,N), the coordinates of the vertices.
%
%    Output, real AREA, the area of the polygon.
%
  area = 0.0;

  for i = 1 : n - 2

    t(1:2,1:3) = [ v(1:2,i)'; v(1:2,i+1)'; v(1:2,n)' ]';

    area_triangle = triangle_area_2d ( t );

    area = area + area_triangle;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ft_standalone(varargin)

% FT_STANDALONE is the entry function of the compiled FieldTrip application.
% The compiled application can be used to execute FieldTrip data analysis
% scripts.
%
% This function can be started on the interactive MATLAB command line as
%   ft_standalone script.m
%   ft_standalone script1.m script2.m ...
%   ft_standalone jobfile.mat
% or after compilation on the Linux/macOS command line as
%   fieldtrip.sh <MATLABROOT> script.m
%   fieldtrip.sh <MATLABROOT> script1.m script2.m ...
%   fieldtrip.sh <MATLABROOT> jobfile.mat
%
% It is possible to pass additional options on the MATLAB command line like
% this on the MATLAB command line
%   ft_standalone --option value scriptname.m
% or on the Linux/macOS command line
%   fieldtrip.sh <MATLABROOT> --option value scriptname.m
% The options and their values are automaticallly made available as local
% variables in the script execution environment.
%
% See also FT_COMPILE_STANDALONE

% Copyright (C) 2011-2018 by Robert Oostenveld
%
% This file is part of FieldTrip.
%
% megconnectome 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.
%
% megconnectome 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 megconnectome.  If not, see <http://www.gnu.org/licenses/>.

% this function uses assignin/evalin. The alternative is to use assign/eval
% but then the script execution might collide with local variables inside this
% megconnectome function
workspace = 'caller';

if nargin==0
  % show the help for the compiled application
  fprintf('You should start the compiled FieldTrip application like this:\n');
  fprintf('  fieldtrip <MATLABROOT>  script.m\n');
  fprintf('  fieldtrip <MATLABROOT>  script1.m script2.m ...\n');
  fprintf('  fieldtrip <MATLABROOT>  jobfile.mat\n');
  fprintf('\n');
  return
end

% separate the --options from the filenames
[options, varargin] = hcp_getopt(varargin{:});

if any(strcmp(options(1:2:end), 'version'))
  % only show the provenance information, do not run anything
  varargin = {};
end

% get the general provenance information (e.g. user and hostname)
prov = hcp_provenance;

for i=1:length(varargin)
  fname = varargin{i};

  % get the specific provenance information for the script
  % prov.script{i}.filename = fname;
  % prov.script{i}.md5sum   = hcp_md5sum(fname);

  % The following section with "which" does not work in the compiled application
  %
  %   if isempty(which(fname))
  %     error('The file %s cannot be found\n', fname);
  %   end
  %   % ensure that the script name includes the full path and extension
  %   fname = which(fname);

  [p, f, ext] = fileparts(fname);

  switch ext
    case '.m'
      if ~exist(fname, 'file')
        error('The script %s cannot be found\n',fname);
      end

      fid = fopen(fname);
      if fid == -1, error('Cannot open %s',fname); end
      S = fscanf(fid,'%c');
      fclose(fid);

      % store the script content in the provenance
      % prov.script{i}.content  = sprintf('\n%s\n', S);

      % capture all screen output in a diary file
      diary off
      diaryfile = tempname;
      diary(diaryfile);

      try
        % ensure that subsequent scripts do not interfere with each other
        % keep the ft_default global variable
        global ft_default
        localcopy = ft_default;
        evalin(workspace,'clear all');
        global ft_default
        ft_default = localcopy;

        % make the options available as local variables
        for j=1:2:length(options)
          key = options{j};
          val = options{j+1};

          if ismember(val(1), {'[', '{'})
            % the value is a string that represents an array
            evalin(workspace, sprintf('%s = %s;', key, val));
          elseif ismember(val(1), {'1', '2', '3', '4', '5', '6', '7', '8', '9', '0', '-', '+'}) && all(ismember(cellstr(val(:)), {'1', '2', '3', '4', '5', '6', '7', '8', '9', '0', '-', '+', '.', 'e'}))
            % the value is a string that represents a number
            evalin(workspace, sprintf('%s = %s;', key, val));
          else
            assignin(workspace, key, val);
          end
        end

        % add some information about the script to the global ft_default variable
        % this will be picked up by hcp_write_provenance
        % ft_default.prov.script = prov.script{i};

        % evaluate this script
        evalin(workspace,S);

        % remove the provenance information about the script
        % ft_default.prov.script = [];

      catch err
        % remove the provenance information about the script
        % ft_default.prov.script = [];

        fprintf('Execution failed: %s\n', fname);
        rethrow(err);
      end

    case '.mat'
      % The job input consists of a file like this
      %   /Users/robert/tmp/robert_mbp_p36978_b1_j001_input.mat
      % which after execution results in an output file
      %   /Users/robert/tmp/robert_mbp_p36978_b1_j001_output.mat

      if length(fname)<10
        error('The jobfile %s is incorrect\n',fname);
      elseif ~strcmp(fname(end-9:end), '_input.mat')
        error('The jobfile %s is incorrect\n',fname);
      end

      jobid = fname(1:end-10);
      [p, jobid] = fileparts(jobid);
      cd(p);

      % read the job input arguments from jobid_input.mat, delete the input file and execute the job
      qsubexec(jobid);

      if isempty(getenv('PBS_ENVIRONMENT'))
        % this suggests that the job is running outside of the torque environment, e.g. for local testing
        % execution of qsubexec is expected to result in a stdout and stderr log file from torque
        logout       = fullfile(p, sprintf('%s.o', jobid));
        logerr       = fullfile(p, sprintf('%s.e', jobid));
        fclose(fopen(logout, 'w'));
        fclose(fopen(logerr, 'w'));
      end

      % read the job results from jobid_output.mat, delete the output file
      % FIXME what to do with the output variable?
      argout = qsubget(jobid);
      disp(argout)

    otherwise
      error('Unknown input format for %s, should be *.m or *.mat\n',fname);
  end % switch type of input argument

  % force close all figures that were created
  close all

  diary off
  fid = fopen(diaryfile);
  if fid == -1, error('Cannot open %s',diaryfile); end
  S = fscanf(fid,'%c');
  fclose(fid);
  delete(diaryfile);

  % add the screen output to the provenance structure
  % prov.script{i}.stdout = sprintf('\n%s\n', S);

end % for each of the input arguments

% display the provenance structure on screen as an XML stream
% disp(fileexchange_struct2xml(struct('megconnectome', prov)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function assign(key, val)
assignin('caller', key, val);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [CRLB,J]=getEstimatorMinMSEBound(RInv,stateJacob,statJacob,biasVec,biasJacob,numMeasDims,checkPosDef)
%%GETESTIMATORMINMSEBOUND This evaluates a generalization of the Cramer-Rao
%       lower bound (CRLB) specifically for the case where measurements are
%       corrupted with multivariate Gaussian noise. The measurement model
%       is z=h(x)+w, where x is the (deterministic) state, h is the
%       (possible nonlinear) measurement function and w is zero-mean
%       Gaussian noise with covariance matrix R. h can vary between
%       measurements that are being fused. This offers a lower bound on the
%       mean-squared error of a biased multivariate statistic. The
%       statistic is just a transformation of the state, so zStat=f(x).
%       Biased means that if T(x) is the estimator, then the expecttation
%       E{T(x)-f(x)}=b(x), where b(x) is the bias. If one doesn't know the
%       bias and its gradient, omitting those terms provides a lower bound
%       on the covariance matrix of an unbiased statistic. If the statistic
%       Jacobian is omitted, then the bound is directly on the state x.
%       The bias term can be useful when one already has an estimator with
%       a known or approximated bias and they want to approximate its
%       accuracy with the MSE matrix bound rather than analytically
%       evaluating its true MSE.
%
%INPUTS: RInv The zDimXzDimXnumMeas set of inverse covariance matrices
%          associated with the multivariate Gaussian noise corrupting each
%          of the numMeas measurements. If the dimensionality of the
%          measurements varies, then zDim is the maximum
%          dimensionality of any measurement and then numMeasDims is
%          required so that if sel=1:numMeasDims(k), then RInv(sel,sel,k)
%          is the submatrix used for the kth measurement.
% stateJacob The zDimXxDimXnumMeas Jacobian matrices of derivatives of the
%          measurement function h taken with respect to the elements of
%          the target state for every measurement. If a single zDimXxDim
%          matrix is passed, then it is assumed that this is the same for
%          all numMeas measurement. If the dimensionality of the
%          measurements varies, then zDim is the maximum
%          dimensionality of any measurement and then numMeasDims is
%          required so that if sel=1:numMeasDims(k), then
%          stateJacob(sel,:,k) is the submatrix used for the kth
%          measurement.
% statJacob The statDimXxDim matrix of derivatives of the the statistic
%          function f(x) taken with respect to the state x. If this is
%          omitted or an empty matrix is passed, then an xDimXxDim identity
%          matrix is used, which is the same as h(x)=x.
%  biasVec A statDimX1 vector, if provided. This is the bias of the assumed
%          estimator of the statistics. If this value and the next value
%          are omitted or empty matrices are provided, then the CRLB for an
%          unbiased statistic is computed.
% biasJacob A statDimXxDim matrix of derivatives of the estimator bias b(x)
%          taken with respect to the state x. If this is provided, then
%          biasVec must be provided. If omitted or an empty matrix is
%          passed, then this is assumed to be 0.
% numMeasDims If all measurement shave the same dimensionality , then this
%          input should be omitted or an empty matrix passed. Otherwise,
%          this is a length numMeas vector that specified how many
%          dimensions each measurement has.
% checkPosDef If this is true, then a check is performed as to whether the
%          Fisher information matrix in the algorithm is full rank.
%          If it is not, then the CRLB matrix will be returned empty and no
%          attempt will be made to invert it. The default if omitted or an
%          empty matrix is passed is false.
%
%OUTPUTS: CRLB The statDimXstatDim CRLB matrix or, if the FIM was singular,
%              and checkPosDef is true, an empty matrix.
%            J The Fisher infromation matrix that went into computing CRLB.
%
%A simple derivation of the standard multivariate CRLB is derived in [1].
%If in the derivation there, one replaces the difference T_k(x)-x_k (the
%difference of the best estimator of the kth component of x and the state)
%with T_k(x)-f_k(x), where f_k(x) is the kth component of a desired
%statistic and thus T_k(x) becomes an estimator of a component of a
%statistic, then one will derive the CRLB of a statistic instead of just of
%the state. In the scalar case, this was already in Rao's original paper in
%[2].
%
%To consider a biased statistic, one then replaces the assumption in the
%paper of a zero expected value: E{T_k(x)-x_k}=0 with
%E{T_k(x)-f_k(x)}=b_k(x). The rest of the derivation in [1] is pretty much
%the same, making sure that one is evaluating Cov{Z} and not just E{Z*Z'}.
%Completing the derivation just using E{Z*Z'} and not Cov{Z} will lead one
%to have a lower bound that is smaller by a bias*bias' term.
%
%Note that is is possible for a biased estimator to have a lower MSE than
%an unbiased estimator. One such example is given in [3].
%
%EXAMPLE 1:
%This is a simple example, where z=x+w with w having the identity matrix as
%the covariance matrix. We take x to be 2X1. We want to estimate
%f(x)=x'*x. The assumed biased estimator is T(z)=z'*z+sum(z)/100. (Note
%that even T(z)=z'*z is biased, but we choose to make the bias non-
%constant). The bias in 2D can be found to be b(x)=2+sum(x)/100.
%The analytic mean squared error in 2D is
%MSE=(80002+40001*x(1)^2+2*x(1)*(400+x(2))+x(2)*(800+40001*x(2)))/10000
%We plot the MSE and the CRLB for a fixed x(1), while varying x(2). One can
%see a reasonably good level of agreement, though the bound is not tight.
% MSE=@(x)(80002+40001*x(1)^2+2*x(1)*(400+x(2))+x(2)*(800+40001*x(2)))/10000;
% biasVec=@(x)(2+sum(x)/100);
% biasJacob=@(x)[1/100,1/100];
% statJacob=@(x)2*x';
% stateJacob=eye(2,2);
% RInv=eye(2,2);
% numPts=200;
% x1=-1/2;
% x2=linspace(-5,5,numPts);
% MSEVal=zeros(1,numPts);
% CRLB=zeros(1,numPts);
% for curPt=1:numPts
%     xCur=[x1;x2(curPt)];
%     MSEVal(curPt)=MSE(xCur);
%     CRLB(curPt)=getEstimatorMinMSEBound(RInv,stateJacob,statJacob(xCur),biasVec(xCur),biasJacob(xCur));
% end
% figure(1)
% clf
% hold on
% plot(x2,MSEVal,'-b','linewidth',2)
% plot(x2,CRLB,'--k','linewidth',2)
% legend('Actual MSE','CRLB','location','north')
%
%REFERENCES:
%[1] A. Papoulis and S. U. Pillai, Probability, Random Variables and
%    Stochastic Processes, 4th ed. Boston: McGraw Hill, 2002.
%[2] C. R. Rao, "Information and the accuray attainable in the estimation
%    of statistical parameters," Bulletin of the Calcutta Mathematical
%    Society, vol. 37, no. 3, pp. 81-91, 1945.
%[3] P. Stoica and R. L. Moses, "On biased estimators and the unbiased
%    Cramér-Rao lower bound," Signal Processing, vol. 21, no. 4, pp. 349-
%    350, Dec. 1990.
%
%November 2022 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numMeas=size(RInv,3);
xDim=size(stateJacob,2);
statDim=size(statJacob,1);

if(nargin<7||isempty(checkPosDef))
    checkPosDef=false;
end

if(nargin<6)
    numMeasDims=[];
end

if(nargin<5||isempty(biasJacob))
    biasJacob=zeros(statDim,xDim);
end

if(nargin<4||isempty(biasVec))
    biasVec=zeros(statDim,1);
end

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

if(numMeas>1&&size(stateJacob,3)==1)
    stateJacob=repmat(stateJacob,[1,1,numMeas]);
end

%First get the FIM
J=zeros(xDim,xDim);
if(isempty(numMeasDims))
    %If all measurements are zDim in size
    for k=1:numMeas
        J=J+stateJacob(:,:,k)'*RInv(:,:,k)*stateJacob(:,:,k);
    end
else
    for k=1:numMeas
        zDimCur=numMeasDims(k);
        sel=1:zDimCur;
        J=J+stateJacob(sel,:,k)'*RInv(sel,sel,k)*stateJacob(sel,:,k);
    end
end

if(checkPosDef)
    %Check whether or not J is positive definite.
    if(matrixRank(J)<xDim)
        CRLB=[];
        return;
    end
end

CRLB=(statJacob+biasJacob)*inv(J)*(statJacob+biasJacob)'+biasVec*biasVec';

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 size = calculateSize( fileName )
%CALCULATESIZE Calculates the size of a MATLAB file for use in the game Cody
%   CALCULATESIZE(FILENAME) determines the number of nodes in the parse tree
%   associated with the file FILENAME.
%
%   For more information, visit http://www.mathworks.com/matlabcentral/cody
%
%   Copyright 1984-2012 The MathWorks, Inc.

    % MTREE is an undocumented function. It will change in the future
    % versions of MATLAB.
    t = mtree(fileName,'-file');
    size = length(t.nodesize);
end


text:
function s = ymdf_to_s_english ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_S_ENGLISH writes an English YMDF date into a string.
%
%  Format:
%
%    BC OS YYYY/MM/DD.FF
%    AD OS YYYY/MM/DD.FF
%    AD NS YYYY/MM/DD.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, string S, a representation of the date.
%
  if ( y < 0 )
    s = sprintf ( 'BC OS %d/%02d/%05.2f', -y, m, d + f );
  elseif ( y < 1752 || ( y == 1752 && m < 9 ) || ( y == 1752 && m == 9 && d < 3 ) )
    s = sprintf ( 'AD OS %d/%02d/%05.2f', y, m, d + f );
  else
    s = sprintf ( 'AD NS %d/%02d/%05.2f', y, m, d + f );
  end

  return
end


text:
% 2D mobile robot example.
% by Olzhas Adiyatov
% 08/28/2013

map = struct('name', 'bench_june1.mat', 'start_point', [-12.5 -5.5], 'goal_point', [7 -3.65]);
max_iter =  10e3;
max_nodes = 3e3;
is_benchmark = false;
rand_seed = 40;
variant = 'FNSimple2D';
result = rrt_star_fn(map, max_iter, max_nodes, is_benchmark, rand_seed, variant);

author: olzhas


text:
% MAIN.m  --  Lesson 8 -- Create GIF
%
% This script performs an simulation and animation of a cart-pole as it
% moves passively along a horizontal track.
%
% In this lesson I've removed the "real time" animation and replaced it
% with a script that generates a MP4 video file instead.
%
% NOTE: The video writer object will only write mp4 video files to Windows
% and Mac platforms (no linux for now)
%

clc; clear; clear global;

%%%% Initial State
z0 = [
    0.0;   %horizontal position
    (pi/180)*80;  %pendulum angle (wrt gravity)
    0.3;   %horizontal velocity
    0.0];  %pendulum angular rate

%%%% Physical Parameters  (big mass and inertia for "slow" physics)
p.m1 = 18.0;  % (kg) Cart mass
p.m2 = 12.0;  % (kg) pole mass
p.g = 9.81;  % (m/s^2) gravity
p.l = 5.0;   % (m) pendulum (pole) length

%%%% Time vector
t = linspace(0,10,250);  %Simulation time stamps

%%%% Function Handle
dynFun = @(t,z)( cartPoleDynamics(z, p) );

%%%% Simulate the system!
options = odeset(...
    'RelTol',1e-8, ...
    'AbsTol',1e-8);
[~, z] = ode45(dynFun, t, z0, options);   %  <-- This is the key line!
z = z';

%%%% Plots:
figure(1); clf; hold on;
plotCartPole(t,z);  %Moved plotting to its own function

%%%% Animation:

% Convert states to cartesian positions:
pos = cartPolePosition(z,p);
x1 = pos(1,:);
y1 = pos(2,:);
x2 = pos(3,:);
y2 = pos(4,:);

% Plotting parameters:
p.w = 0.6*p.l;  %Width of the cart
p.h = 0.4*p.l;  %Height of the cart
p.r = 0.1*p.l;  % Radius of the pendulum bob

% Compute the extents of the drawing, keeping everything in view
padding = 0.2*p.l;  %Free space around edges
xLow = min(min(x1 - 0.5*p.w,  x2 - p.r)) - padding;
xUpp = max(max(x1 + 0.5*p.w,  x2 + p.r)) + padding;
yLow = min(min(y1 - 0.5*p.h,  y2 - p.r)) - padding;
yUpp = max(max(y1 + 0.5*p.w,  y2 + p.r)) + padding;
extents = [xLow,xUpp,yLow,yUpp];

% Create and clear a figure:
figHandle = figure(2); clf; hold on;

% Compute the verticies of a star, just for fun;
star = getStarVerticies(7,0.5);  % 7 verticies, spoke ratio of 0.5
p.star = 0.6*p.r*star;  %Rescale the star;

% Set up MP4 and log data
vidObj = VideoWriter('cartPoleAnimation.mp4','MPEG-4');
vidObj.FrameRate = 25;
vidObj.Quality = 100;
open(vidObj);

frameRate = vidObj.FrameRate;
nFrame = floor(frameRate*t(end));
frameDelay = 1/frameRate;
time = 0;
for i=1:nFrame
    % Compute the position of the system at the current real world time
    posDraw = interp1(t',pos',time')';

    % Redraw the image
    drawCartPole(time,posDraw,extents,p);

    % Write data to video file
    writeVideo(vidObj,getframe(figHandle));

    % time step system to next frame:
    time = time + frameDelay;
end

close(vidObj);

author: MatthewPeterKelly

repo: dscTutorials


text:
function [region_parts, orig_sp, varargout] = spagglom(I_rgb, opts, varargin)
% The main function for method described in "Generating object segmentation
% proposals using global and local search".
%
% INPUTS
%
%   'I_rgb' is an rgb image of type uint8
%   'opts' is an options cell array, initialized by running script
%       spagglom_options.m
%   Additionally, precalculated superpixelation for the image and/or pixel
%       features can be supplied - please see the code for details. Script
%       test_recalls.m also contains few examples.
%
% OUTPUTS
%
%   'region_parts' is a cell array, each item of which is a set of natural
%       numbers
%   'orig_sp' is a cell array. Each element represents a superpixel,
%       containing various information. The indexing corresponds to
%       region_parts. For example region_parts{1} = [4, 8] means that the
%       first region is made up of superpixels 4 and 8, which can be
%       accessed by orig_sp{4} and orig_sp{8}.
%   Additionally, histograms for each generated region proposal are
%       returned if third output variable is specified. See the very bottom
%       of this function for details.
%
%
% Requires VLFeat. Version 9.16 was used for developing.

%% Parse varargin
I_seg = [];
words = [];
histograms = [];
argl = length(varargin);
if mod(argl,2) ~= 0 % string-value pairs
    error('Give pairs of extra arguments.');
end
for argi = 0:((argl/2) - 1)
    arg_str = varargin(2*argi + 1);
    arg_val = varargin(2*argi + 2);
    switch arg_str{1}
        case 'I_seg' % Integer map of image superpixelation
            I_seg = arg_val{1};
        case 'words' % Pixelwise features
            words = arg_val{1};
        case 'histograms' % Histograms for each superpixel. Using this is not as straightforward as the above two
            histograms = arg_val{1};
        otherwise
            error('Argument error.');
    end
end

[h, w, ~] = size(I_rgb); % image size

opts.sp_reindexed = 0;

region_parts{1} = []; % suppresses errors if no region_parts or graphcut_parts are generated below
graphcut_parts{1} = [];

%% Initial segmentation

% Calculate the initial segmentation
[sp, K] = initial_seg(I_rgb, opts, I_seg);

orig_sp = sp; % make a copy of sp. This will be used by the scoring routine.

%% Features
if isempty(histograms) % by default, this is true
    % Compute pixel-wise features
    if isempty(words)
        [words, k] = compute_features(I_rgb, I_rgb, 'rgb', opts);
    else % use user 'words'
        k = [];
        for r = 1:length(words)
            k(r) = max(words{r}(:)); % histogram sizes
        end
    end

    % Compute histograms from features and add them to superpixels
    sp = compute_histograms(sp, words, k, h, w);

else % use user supplied histograms
    assert(length(histograms) == length(sp)); % one set of histograms for each superpixel
    hist_num = length(histograms{1});
    % Important: You must setup 'opts.features' in spagglom_options.m to
    % match what type of histograms you are supplying, in correct order.

    for r = 1:length(sp) % for each superpixel
        for fn = 1:hist_num % for each histogram
            sp{r}.hist{fn} = histograms{r}{fn}; % copy the histogram values
        end
    end
end

%% Agglomeration of superpixels

% Calculate similarity scores between all neighboring superpixel pairs
scores = similarity_scores(sp, K, opts);

%%
%% Run the non-branching part of the algorithm. "Refining the superpixelation". <<< PHASE 1 >>>
opts.phase = 1; % Critically important line

%[sp, K, scores, region_parts_phase1] = spagglom_sub(sp, K, scores, 0, opts);
[sp, K, ~, ~] = spagglom_sub(sp, K, scores, 0, opts); % returns the current state of 'sp' and 'K'

%% Update indexing of superpixels

if 1 % update indexing for speedup. After this region_parts_primary cannot be used!
    [sp, K, scores] = update_sp_indexing(sp, K, h, w, opts); % makes rest of the algorithm (and scoring) faster by decreasing amount of parts and sp
    orig_sp = sp; % Replace the orig_sp
    opts.sp_reindexed = 1; % a note that this step was performed
else
    % Not updating the indexing allows for returning the
    % region_parts_phase1 above, because they can be represented using the
    % orig_sp defined earlier.
end

%%
%% Run the local and global search <<< PHASE 2 >>>
opts.phase = 2; % Critically important line
sp = setup_on_edge(sp, h, w); % Setup on_edge variables. Required before graphcut_regions()

graphcut_parts{1} = graphcut_regions(sp, K, opts);

initial_regions_created = 0; % default 0
[~, ~, ~, region_parts{1}] = spagglom_sub(sp, K, scores, initial_regions_created, opts); % CREATE INITIAL REGIONS HERE (YES 0 <-> 1 NO)
initial_regions_created = 1; % may prevent unnecessary duplicates if the above local search is repeated below

% Example of how to run in another color space
% [sp_nrgb, scores_nrgb] = change_color_space('nrgb', I_rgb, sp, K, im_num, h, w, [], opts);
% [~, ~, ~, region_parts{2}] = spagglom_sub(sp_nrgb, K, scores_nrgb, initial_regions_created, opts);
% graphcut_parts{2} = graphcut_regions(sp_nrgb, K, opts);

%% Combine all found regions
% Above you can collect region_parts{1}, region_parts{2}, region_parts{3}, ... as
% many as you want. You can only add region_parts_phase1 if you did NOT
% update superpixel indexing above.
%region_parts = horzcat(region_parts_phase1, region_parts{:}, graphcut_parts{:});
region_parts = horzcat(region_parts{:}, graphcut_parts{:});

%% Keep only unique entries
regions_bin = false(length(region_parts), length(orig_sp)); % Represent regions as binary masks of the superpixelation orig_sp
for ku = 1:length(region_parts)
    regions_bin(ku, region_parts{ku}) = 1;
end

% Use unique() to remove duplicates
[~, ue, ~] = unique(regions_bin, 'rows'); % This could perhaps be optimized further by summing rows of regions_bin and running unique for groups having the same sum (the same amount of parts)
region_parts = region_parts(ue); % The return value

%% Output region histograms if desired
if nargout == 3
    varargout{1} = cell(1,length(region_parts));
    for ku = 1:length(region_parts)
        varargout{1}{ku} = merge_histograms(orig_sp, region_parts{ku});
    end
    %varargout{1}{10}{2} % Example: gives the second histogram (the rgb
    % histogram in default case) of 10th region
end

author: Cloud-CV

repo: object-proposals


text:
function [x,y,w] = gqcircsegm ( n, omega, r )

%
%% GQCIRCSEGM computes the nodes and weights of a product gaussian formula
% on a circle segment of a disk centered at the origin
% with angles in [0,omega]

% uses the routines:
%
% r_jacobi.m, gauss.m from
% www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
%
% trigauss.m
% http://www.math.unipd.it/~marcov/mysoft/trigauss.m
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    Gaspare Da Fies, Marco Vianello,
%    University of Padova
%
%  Parameters:
%
%    Input, integer n: algebraic degree of exactness.
%
%    Input, real omega: half-length of the angular interval, 0<omega<=pi
%
%    Input, r: radius of the disk
%
%    Output, xyw: (ceil((n+1)/2) x ceil((n+2)/2) x 3 array of (xnodes,ynodes,weights)

%
%  Compute the trigonometric gaussian formula on the arc.
%
  tw = trigauss ( n + 2, 0, omega );
%
%  Compute the algebraic gaussian formula on [-1,1];
%
  ab = r_jacobi ( ceil((n+1)/2), 0, 0 );
  xw = gauss ( ceil((n+1)/2), ab );
%
%  Create the grid.
%
  [ t, theta ] = meshgrid ( xw(:,1), tw(:,1) );
  [ w1, w2 ]   = meshgrid ( xw(:,2), tw(:,2) );
%
%  Nodal cartesian coordinates and weights.
%
  s = sin ( theta(:) );
  x = (r * cos ( theta(:) ))';
  y = (r * t(:) .* s)';
  w = (r^2 * s.^2 .* w1(:) .* w2(:))';

  return
end


text:
function g = SigmaNoisy(Reference,Target,Clean)

%ASSUMES THAT TARGET IMAGE IS NOT REFERENCE IMAGE OR CLEAN IMAGE

sumcn=0;

[m1,n1] = size(Reference);
[m2,n2] = size(Target);
[m3,n3] = size(Clean);

r=[m1;m2;m3;];
c=[n1/3;n2/3;n3/3;];

rows=min(r);
cols=min(c);
for i=1:3

    x1=(Reference(1:rows,1:cols,i));
    y=(Target(1:rows,1:cols,i));
    x2=(Clean(1:rows,1:cols,i));

    avg = (x1 + x2)/2;
    diffNoisy=y-avg;
    diffFull=x1-x2;

    varx1 = var(diffFull(:))/4;

    %The sigma values of Target image RGB channels
    vary= max(var(diffNoisy(:)) - varx1,0);

    g(:,i) = sqrt(vary);
    sumcn= sumcn+vary;

end
sigmaNoisy = sqrt(sumcn/3);

g(:,4)=sigmaNoisy;

author: lbasek

repo: image-denoising-benchmark


text:
function [ u, nullty, ifault ] = cholesky ( a, n, nn )

%*****************************************************************************80
%
%% CHOLESKY computes the Cholesky factorization of a PDS matrix.
%
%  Discussion:
%
%    For a positive definite symmetric matrix A, the Cholesky factor U
%    is an upper triangular matrix such that A = U' * U.
%
%    This routine was originally named "CHOL", but that conflicts with
%    the name of a built in MATLAB routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by Michael Healy.
%    Modifications by AJ Miller.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Michael Healy,
%    Algorithm AS 6:
%    Triangular decomposition of a symmetric matrix,
%    Applied Statistics,
%    Volume 17, Number 2, 1968, pages 195-197.
%
%  Parameters:
%
%    Input, real A((N*(N+1))/2), a positive definite matrix
%    stored by rows in lower triangular form as a one dimensional array,
%    in the sequence
%    A(1,1),
%    A(2,1), A(2,2),
%    A(3,1), A(3,2), A(3,3), and so on.
%
%    Input, integer N, the order of A.
%
%    Input, integer NN, the dimension of the array used to store A,
%    which should be at least (N*(N+1))/2.
%
%    Output, real U((N*(N+1))/2), an upper triangular matrix,
%    stored by columns, which is the Cholesky factor of A.
%
%    Output, integer NULLTY, the rank deficiency of A.  If NULLTY is zero,
%    the matrix is judged to have full rank.
%
%    Output, integer IFAULT, an error indicator.
%    0, no error was detected;
%    1, if N < 1;
%    2, if A is not positive semi-definite.
%    3, if NN < (N*(N+1))/2.
%
%  Local Parameters:
%
%    Local, real ETA, should be set equal to the smallest positive
%    value such that 1.0 + ETA is calculated as being greater than 1.0 in the
%    accuracy being used.
%

%
%  Initialize output.
%
  u = [];
  nullty = 0;
  ifault = 0;

  eta = 1.0E-09;

  if ( n <= 0 )
    ifault = 1;
    return
  end

  if ( nn < floor ( ( n * ( n + 1 ) ) / 2 ) )
    ifault = 3;
    return
  end

  u(1:nn) = 0.0;
  j = 1;
  k = 0;
  ii = 0;
%
%  Factorize column by column, ICOL = column number.
%
  for icol = 1 : n

    ii = ii + icol;
    x = eta * eta * a(ii);
    l = 0;
    kk = 0;
%
%  IROW = row number within column ICOL.
%
    for irow = 1 : icol

      kk = kk + irow;
      k = k + 1;
      w = a(k);
      m = j;

      for i = 1 : irow - 1
        l = l + 1;
        w = w - u(l) * u(m);
        m = m + 1;
      end

      l = l + 1;

      if ( irow == icol )
        break
      end

      if ( u(l) ~= 0.0 )

        u(k) = w / u(l);

      else

        u(k) = 0.0;

        if ( abs ( x * a(k) ) < w * w )
          ifault = 2;
          return
        end

      end

    end
%
%  End of row, estimate relative accuracy of diagonal element.
%
    if ( abs ( w ) <= abs ( eta * a(k) ) )

      u(k) = 0.0;
      nullty = nullty + 1;

    else

      if ( w < 0.0 )
        ifault = 2;
        return
      end

      u(k) = sqrt ( w );

    end

    j = j + icol;

  end

  return
end


text:
function inds = fernsInds( data, fids, thrs )
% Compute indices for each input by each fern.
%
% USAGE
%  inds = fernsInds( data, fids, thrs )
%
% INPUTS
%  data     - [NxF] N length F binary feature vectors
%  fids     - [MxS] feature ids for each fern for each depth
%  thrs     - [MxS] threshold corresponding to each fid
%
% OUTPUTS
%  inds     - [NxM] computed indices for each input by each fern
%
% EXAMPLE
%
% See also fernsClfTrain, fernsClfApply
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.50
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

inds = fernsInds1( data, fids, thrs );

%%% OLD MATLAB CODE -- NOW IN MEX
% [M,S]=size(fids); N=size(data,1);
% inds = zeros(N,M,'uint32');
% for n=1:N
%   for m=1:M
%     for s=1:S
%       inds(n,m)=inds(n,m)*2;
%       if( data(n,fids(m,s))<thrs(m,s) )
%         inds(n,m)=inds(n,m)+1;
%       end
%     end
%   end
% end
% inds=inds+1;
% end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function varargout = imsmooth(varargin)
% VL_IMSMOOTH  Smooth image
%   J = VL_IMSMOOTH(I,SIGMA) convolves the image I by an isotropic
%   Gaussian kernel of standard deviation SIGMA.  I must be an array
%   of doubles. IF the array is three dimensional, the third dimension
%   is assumed to span different channels (e.g. R,G,B). In this case,
%   each channel is convolved independently.
%
%   VL_IMSMOOTH() accepts the following options:
%
%   Kernel::
%     Selects between GAUSSIAN and TRIANGULAR kernels. The triangular
%     kernel support has 2*SIGMA-1 sampels. Kernels have unit mass.
%
%   Padding::
%     Selects between ZERO or CONTINUITY padding method to handle the
%     image boundaries. ZERO extends the input image with zeroes
%     around the border, and CONTINUITY extends the image with
%     constant pixels.
%
%   Step::
%     Sets the subsampling step. A subsampling step of STEP pixels
%     causes J(1:STEPS:end, 1:STEPS:end, :) to be computed. This is
%     useful to downsample the image.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_imsmooth(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [varargout]=febStruct2febFile(FEBioInputStruct)

% function [varargout]=febStruct2febFile(FEBioInputStruct)
% ------------------------------------------------------------------------
%
% This function uses the input FEBioInputStruct to generate an FEBio .feb file.
%
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/05/27: Updated for GIBBON
% 2015/05/09: Added biphasic capabilities
% 2015/05/10: Added traction load capabilities
% 2017/08/18: Adding support to allow incomplete feb file creation (only
% write fields provided).
% 2018/01/17: Added support for febiospec 2.5 and made it default
% 2018/01/18: Decided to make 2.5 available in new seperate system
%------------------------------------------------------------------------

%%

warning('febStruct2febFile is depricated and will be removed in future releases. Update your codes to use febioStruct2xml, see HELP_febioStruct2xml (furthermore febio_spec version 2.5 is recommended)');

%%

dispStartTitleGibbonCode('Writing FEBio XML object');

%% Set default display setting if missing
if ~isfield(FEBioInputStruct,'disp_opt')
    FEBioInputStruct.disp_opt=0;
end

%% Initialize docNode object
domNode = com.mathworks.xml.XMLUtils.createDocument('febio_spec'); %Create the overall febio_spec field

febioSpecVersion='2.0';

%Set febio_spec
febio_spec = domNode.getDocumentElement;
if ~isfield(FEBioInputStruct,'febio_spec')
    FEBioInputStruct.febio_spec.version=febioSpecVersion;
elseif ~isfield(FEBioInputStruct.febio_spec,'version')
    FEBioInputStruct.febio_spec.version=febioSpecVersion;
end
febio_spec.setAttribute('version',FEBioInputStruct.febio_spec.version); %Adding version attribute

disp(['Using febio_spec: ',FEBioInputStruct.febio_spec.version]);

%% Add comment if present
if isfield(FEBioInputStruct,'commentField')
    commentString = FEBioInputStruct.commentField;
else %Default comment
    commentString = ['Created using GIBBON, ',datestr(now)];
end
commentNode = domNode.createComment(commentString);
febio_spec.appendChild(commentNode);

%% DEFINING MODULE LEVEL
if ~isfield(FEBioInputStruct,'Module')
    FEBioInputStruct.Module.Type='solid'; %Use solid as default module
end
domNode=addModuleLevel_FEB(domNode,FEBioInputStruct);

%% DEFINE CONTROL SECTION
if isfield(FEBioInputStruct,'Control')
    domNode=addControlLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINING GLOBALS LEVEL
if isfield(FEBioInputStruct,'Globals')
    domNode=addGlobalsLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINING MATERIAL LEVEL
if isfield(FEBioInputStruct,'Materials')
    domNode=addMaterialLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINING GEOMETRY LEVEL

if isfield(FEBioInputStruct,'Geometry')
    writeMethod=1; %Use to switch between text mode (1) or XML model (2)
    switch writeMethod
        case 1 % TEXT FILE PARSING (faster for large models)
            domNode=addGeometryLevel_TXT(domNode,FEBioInputStruct);
        case 2 %XML PARSING
            domNode=addGeometryLevel_FEB(domNode,FEBioInputStruct);
    end
end

%% DEFINE BOUNDARY CONDITIONS LEVEL AND LOADS LEVEL
if isfield(FEBioInputStruct,'Boundary')
    [domNode]=addBoundaryLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE LOAD LEVEL
if isfield(FEBioInputStruct,'Loads')
    [domNode]=addLoadsLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE CONTACT
if isfield(FEBioInputStruct,'Contact')
    [domNode]=addContactLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE CONSTRAINTS LEVEL
if isfield(FEBioInputStruct,'Constraints')
    [domNode]=addConstraintsLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE LOADDATA LEVEL
if isfield(FEBioInputStruct,'LoadData')
    [domNode]=addLoadDataLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE STEP LEVEL
if isfield(FEBioInputStruct,'Step')
    [domNode]=addStepLevel_FEB(domNode,FEBioInputStruct);
end

%% DEFINE OUTPUT LEVEL
domNode=addOutputLevel_FEB(domNode,FEBioInputStruct);

%% CREATE OUTPUT OR EXPORT XML FILE

switch nargout
    case 0
        disp('Writing .feb file');
        xmlwrite_xerces(FEBioInputStruct.run_filename,domNode); %Custom XML write function
%         if isfield(FEBioInputStruct,'topCommentLine')
%             exportFEB_XML(FEBioInputStruct.run_filename,domNode,FEBioInputStruct.topCommentLine); % Saving XML file
%         else
%             exportFEB_XML(FEBioInputStruct.run_filename,domNode); % Saving XML file
%         end
    case 1
        varargout{1}=domNode;
end
dispDoneGibbonCode;

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/>.


text:
function vEst=RROnlyStaticVelEst(rr,xTx,xRx,zTar)
%%RRONLYSTATICVELEST Perform unweighted least-squares estimation of a
%           target velocity vector given only range rate measurements from
%           different bistatic channels (or from multiple receivers if the
%           target is the transmitter). This function can work in 2D and
%           3D space and it will produce a least-squares estimate when more
%           than the minimum number of range rate measurements needed is
%           given. This uses a non-relativistic model for the range rates
%           and ignores possible atmospheric effects.
%
%INPUTS: rr A numMeasX1 or 1XnumMeas vector of range rates. numMeas>=3 for
%           the velocity vector to be observable in 3D and numMeas>=2 in
%           2D.
%       xTx A 6XnumMeas matrix (in 3D or a 4XnumMeas matrix in 2D) of
%           stacked transmitter position and velocity vectors. If the
%           target is the transmitter, then an empty matrix should be
%           passed. xTx(:,n)=[x;y;z;xDot;yDot;zDot] is the state (3
%           position and 3 velocity components) corresponding to the nth
%           range rate measurement. If all of the transmitters are in the
%           same spot, then a single 6X1 or 4X1 vector can be passed.
%       xRx A 6XnumMeas matrix (in 3D or a 4XnumMeas matrix in 2D) of
%           stacked transmitter position and velocity vectors.
%           states1(:,n)=[x;y;z;xDot;yDot;zDot] is the state (3 position
%           and 3 velocity components) corresponding to the nth range rate
%           measurement. If all of the receivers are in the same spot, then
%           a single 6X1 or 4X1 vector can be passed.
%      zTar The 3X1 (in 3D) or 2X1 (in 2D) Cartesian position of the
%           target.
%
%OUTPUTS: vEst The 3X1 (in 3D) or 2X1 (in 2D) least squares Cartesian
%              velocity estimate of the target.
%
%This implements the least squares velocity vector estimation procedure of
%Equation 41 in Section IV E of [1]. The case of the transmitter being the
%target is handled specially (to get rid fo the singulaity).
%
%EXAMPLE 1:
%This is just a simple example showing that for three range rate values,
%one can get back the orignal velocity value in 3D.
% zTar=[0;40e3;40e3];
% vTar=[400;-200;100];
% xTx1=[100;10e3;3e3;50;50;-50];
% xTx2=[0;0;0;0;0;-20];
% xTx3=[10e3;10e3;3e3;100;-100;100];
%
% xRx1=[-10e3;0;3e3;100;100;100];
% xRx2=[0;10e3;30;-80;-200;-20];
% xRx3=xRx2;
%
% states1=[xTx1,xTx2,xTx3];
% states2=[xRx1,xRx2,xRx3];
%
% xTar=[zTar;vTar];
%
% useHalfRange=false;
% rr=zeros(3,1);
% rr(1)=getRangeRate(xTar,useHalfRange,xTx1,xRx1);
% rr(2)=getRangeRate(xTar,useHalfRange,xTx2,xRx2);
% rr(3)=getRangeRate(xTar,useHalfRange,xTx3,xRx3);
%
% vEst=RROnlyStaticVelEst(rr,states1,states2,zTar)
%One will see that vEst=vzTar.
%
%EXAMPLE 2:
%This is the same as example 1, except the demonstration is now occurring
%in 2D.
% zTar=[0;40e3];
% vTar=[400;-200];
% xTx1=[100;10e3;50;50];
% xTx2=[0;0;0;0];
% xTx3=[10e3;10e3;100;-100];
%
% xRx1=[-10e3;0;100;100];
% xRx2=[0;10e3;-80;-200];
% xRx3=xRx2;
%
% states1=[xTx1,xTx2,xTx3];
% states2=[xRx1,xRx2,xRx3];
%
% xTar=[zTar;vTar];
%
% useHalfRange=false;
% rr=zeros(3,1);
% rr(1)=getRangeRate(xTar,useHalfRange,xTx1,xRx1);
% rr(2)=getRangeRate(xTar,useHalfRange,xTx2,xRx2);
% rr(3)=getRangeRate(xTar,useHalfRange,xTx3,xRx3);
%
% vEst=RROnlyStaticVelEst(rr,states1,states2,zTar)
%One will see that vEst=vzTar.
%
%EXAMPLE 3:
%In this example, the range rate is one-way from the target (e.g. the
%target is an emitter). The relative error of the velocity estimate with
%the truth (noiseless scenario) is within finite precision limits.
% zTar=randn(3,1);
% vTar=randn(3,1);
% numMeas=5;
% xRx=randn(6,numMeas);
% xTar=[zTar;vTar];
% useHalfRange=false;
% rr=zeros(numMeas,1);
% for k=1:numMeas
%     rr(k)=getRangeRate(xTar,useHalfRange,xTar,xRx(:,k));
% end
% vEst=RROnlyStaticVelEst(rr,[],xRx,zTar);
% RelErr=max(abs((vTar-vEst)./vTar))
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numMeas=length(rr);
rr=rr(:);

if(size(xTx,2)==1)
    xTx=repmat(xTx,1,numMeas);
end

if(size(xRx,2)==1)
    xRx=repmat(xRx,1,numMeas);
end

posDim=length(zTar);

zRx=xRx(1:posDim,:);
vRx=xRx((posDim+1):(2*posDim),:);

h=bsxfun(@minus,zTar,zRx);
hNorm=sqrt(sum(h.*h,1));
h=bsxfun(@rdivide,h,hNorm);

if(~isempty(xTx))
    %The target is not the transmitter.
    zTx=xTx(1:posDim,:);
    vTx=xTx((posDim+1):(2*posDim),:);
    hi=bsxfun(@minus,zTar,zTx);
    hiNorm=sqrt(sum(hi.*hi,1));
    hi=bsxfun(@rdivide,hi,hiNorm);

    rDotB=rr+sum(h.*vRx,1)'+sum(hi.*vTx,1)';
    Hv=bsxfun(@plus,h',hi');
else
    %The target is the transmitter.
    rDotB=rr+sum(h.*vRx,1)';
    Hv=h';
end

vEst=lsqminnorm(Hv,rDotB);

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.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% EUROPEAN OPTION PRICE COMPARISON (RUN SCRIPT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Compare Methods For European Options Under Hestons Model
% Author:      Justin Kirkby
%
% Methods: 1) Kahl-Jackel-Lord Approach (Fourier, Heston Model)
%          2) PROJ (Kirkby, 2015), European Levy/Heston
%          3) CONV (Lord, Fang, Bervoets, Oosterlee, 2008), European Levy/Heston
%          4) Carr-Madan (2008), European Levy/Heston Pricer
%          5) Regime Switching Fourier PROJ (Cui, Kirkby, Nguyen, 2017) - Stoch Vol / RS Pricer
%          6) Time-Changed Markov Chain (Cui, Kirkby, Nguyen, 2019), assumes rho = 0
%          7) Monte Carlo, Using Lord et al (2010) Low-Bias Schemes
%               ... More to come
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);
addpath('../../PROJ/LEVY/European_Options')
addpath('../../PROJ/LEVY/RN_CHF')
addpath('../../PROJ/LEVY/Helper_Functions')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
call = 1;    %For call use 1 (2 for put)
S_0  = 100;  %Initial price
W    = 100;  %Strike            %NOTE: no error handling in place for extreme values of W (increase grid if strike falls outside)
r    = .00;  %Interest rate (NOTE: set to zero for comparison with Kahl-Jackel-Lord, based on Forward price)
q    = .00;  %dividend yield (NOTE: keep this at zero for now)
T    = 0.5;    %Time (in years)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 2) CHOOSE MODEL PARAMETERS  (Levy Models)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {};

params.v0 = 0.02; % initial variance
params.theta = 0.02;   % long term variance level
params.eta = 1.6;   % rate of variance mean reversion
params.Sigmav = 0.3;   % volatility of variance
params.rho = 0;   % correlation between Brownian motions (NOTE: methods which assume rho=0 will display in output)

modelInput = getModelInput(6, T, r, q, params);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Kahl-Jackel-Lord (KJL) Approach
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/Heston/')
tic
price_KJL = Heston1993KahlJaeckelLordRev3(call, S_0,W,T,0,r,q, params.v0, params.theta, params.rho, params.eta, params.Sigmav);
time_KJL = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PROJ (Kirkby, 2015)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logN  = 12;   %Uses N = 2^logN  gridpoint
L1 = 30;

% ----------------------
N = 2^logN;    % grid roughly centered on [c1 - alph, c1 + alph]
alpha = getTruncationAlpha(T, L1, modelInput, 6);

tic
price_PROJ = PROJ_European(3, N, alpha, r, q, T, S_0, W, call, modelInput.rnCHF, modelInput.c1*T);
time_PROJ = toc;

ref = PROJ_European(3, 2^15, 2*alpha, r, q, T, S_0, W, call, modelInput.rnCHF, modelInput.c1*T);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Time-Changed Markov Chain Approximation (Cui, Kirkby, Nguyen, 2019)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../PROJ/TIME_CHANGED/European/')
addpath('../../PROJ/STOCHASTIC_VOL/Helper_Functions')
params.model = 1;

ParamsCtmc.varGridMult = .01;
ParamsCtmc.gamma = 6;  % Heston gamma = 3 is good for T ~ 1
ParamsCtmc.Nx = 100; %the number of Markov states

ProjParams.order = 3;
ProjParams.alph = 2^2;
ProjParams.N_proj = 2^8;

n = 0; % number of time steps in time disretization... set to 0 to do continuous time version
hFunc = @(u) u;   % tau = int h(X_s) ds
levyExponent = @(z) -0.5*1i*z - 0.5*z.^2;

tic
price_TCMC = PROJ_TimeChanged_Levy_European(r,q,S_0,T,W,call, levyExponent, hFunc, n, params, ParamsCtmc, ProjParams);
time_TCMC = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SV-PROJ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../PROJ/STOCHASTIC_VOL/European/')
addpath('../../PROJ/STOCHASTIC_VOL/Barrier/')
addpath('../../PROJ/STOCHASTIC_VOL/Helper_Functions/')
% This version uses the Stoch Vol pricer for Barrier options to price European (More of a multiple purpose method)
if call == 1
    down = 1; H = S_0 / 8;   % TODO: this needs to account for the variance of the underlying.
else
    down = 0; H = S_0 * 8;
end

N             = 2^10;    %number of points in density expansion... Value grid size is K:=N/2
m_0           = 40;  % number of CTMC grid points
gamma         = 5;  % CTMC grid width param
gridMethod    = 4; gridMultParam = 0.2; M = 1; psi_J = @(u)0*[u>0];
alpha         = 5;

tic
price_SVP = Barrier_StochasticVol_func(N,alpha,call,down,S_0,W,H,M,r,T,m_0,psi_J,1, params, gridMethod, gamma, gridMultParam);
time_SVP = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Carr-Madan Fourier Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/CarrMadan/')
N = 2^15;
tic
price_CM = CarrMadan_European_Price_Strikes(S_0, W, modelInput.rnCHF, N, T, r, q, call);
time_CM = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  CONV Fourier Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/CONV/')
N = 2^16;
tic
price_CONV = CONV_European_Price(S_0, W, modelInput.rnCHF, T, r, call, N, alpha);
time_CONV = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Monte Carlo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Monte_Carlo/')
addpath('../../Monte_Carlo/European/')

tic
N_sim = 10^5; M = 800; disc = exp(-r*T); scheme = 5;
Spath = Simulate_Heston_Euler_Schemes( N_sim, M, T, S_0, r, q, params, scheme);
[price_MC, stdErr] = Price_MC_European_Strikes_func(Spath, disc, call, W );
price_MC_L = price_MC - 2*stdErr; price_MC_U = price_MC + 2*stdErr;
time_MC = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPARE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n---------------------------------------------\n')
fprintf('Method       | Price       |    Err   |  CPU  \n')
fprintf('---------------------------------------------\n')
fprintf('PROJ         | %.8f  | %.2e | %.4f  \n', price_PROJ, abs(price_PROJ-ref), time_PROJ)
fprintf('KJL          | %.8f  | %.2e | %.4f  \n', price_KJL, abs(price_KJL-ref), time_KJL)
fprintf('CONV         | %.8f  | %.2e | %.4f  \n', price_CONV, abs(price_CONV-ref), time_CONV)
fprintf('Carr-Madan   | %.8f  | %.2e | %.4f  \n', price_CM, abs(price_CM-ref), time_CM)
fprintf('TC-MC (rho=0)| %.8f  | %.2e | %.4f  \n', price_TCMC, abs(price_TCMC-ref), time_TCMC)
fprintf('SV-PROJ      | %.8f  | %.2e | %.4f  \n', price_SVP, abs(price_SVP-ref), time_SVP)
fprintf('MC-Euler     |[%.3f,%.3f]| %.2e | %.4f \n', price_MC_L, price_MC_U, abs(price_MC-ref), time_MC)