text:
%% DEMO_volumetric_SED_eval
% This demo was developed as part of the paper:
% Moerman et al. "Novel Hyperelastic Models for Large Volumetric
% Deformations".
%
% The demo features:
% * Implementations of hyperelastic volumetric strain energy density
% functions (SEDs)
% * Visualizations of the SED, hydrostatic stress, and tangent as a
% function of the volume ratio.

%% Keywords
%
% * Strain energy density
% * Volumetric
% * Visualization

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
fontSizeInner=fontSize+15;
fontSizeLabel=fontSize+30;
plotColors=gjet(4);
plotColors=plotColors([1 2 4],:);

lineWidth=3;
gridAlpha=0.3;
LineWidthAxis=2;
legendHeight=0.05;
numXTicks=5;

%% Control parameters

formulationCases=1:12;  %Choose formulation 1:12
k=1; %Default bulk modulus (except for hyperfoam)

J_max=2;
numPoints=2000; %Number of points for plotting
J=linspace(0.1,J_max,numPoints)'; %The volume ratios
xtickRange=linspace(0,max(J),numXTicks); %X-axis tick range

%% Get or set formulation specific data and parameters

for formulationCase=formulationCases

    switch formulationCase
        case 1 %Hencky
            formulationName='Hencky';
            parSet(1)=k; %Bulk modulus
        case 2 %Simo
            formulationName='Simo';
            parSet(1)=k; %Bulk modulus
        case 3 %Bischoff
            formulationName='Bischoff';
            b=2; %Beta
            parSet(1)=k; %Bulk modulus
            parSet(2)=b; %Beta
        case 4 %Modified Ogden
            formulationName='Modified Ogden';
            b=2;%Beta
            parSet(1)=k; %Bulk modulus
            parSet(2)=b; %Beta
        case 5 %Hyperfoam
            formulationName='Hyperfoam';
            mu=1;
            a=2; %Alpha
            b=2; %Beta
            parSet(1)=mu; %Mu
            parSet(2)=a; %Alpha
            parSet(3)=b; %Beta
        case 6 %Doll and Schweizerhoff
            formulationName='Doll and Schweizerhoff';
            a=3; %Alpha
            b=2; %Beta
            parSet(1)=k; %Bulk modulus
            parSet(2)=a; %Alpha
            parSet(3)=b; %Beta
        case 7 %Montella et al.
            formulationName='Montella et al.';
            k1=1;
            beta1=1/8;
            k2=1;
            beta2=1/8;
            m=4;
            parSet(1)=k; %Bulk modulus
            parSet(2)=beta1;
            parSet(3)=k2;
            parSet(4)=beta2;
            parSet(5)=m;
        case 8 %Moerman 1
            formulationName='Moerman 1';
            b1=3;
            b2=2;
            parSet(1)=k; %Bulk modulus
            parSet(2)=b1; %Alpha
            parSet(3)=b2; %Beta
        case 9 %Moerman 1A
            formulationName='Moerman 1A';
            b1=3;
            b2=2;
            q=0.5;
            parSet(1)=k; %Bulk modulus
            parSet(2)=b1; %Alpha
            parSet(3)=b2; %Beta
            parSet(4)=q; %Weigthing factor
        case 10 %Moerman 1B
            formulationName='Moerman 1B';
            b1=3;
            b2=2;
            parSet(1)=k; %Bulk modulus
            parSet(2)=b1; %Alpha
            parSet(3)=b2; %Beta
        case 11 %Moerman 2
            formulationName='Moerman 2';
            J1=max(J)+0.1;
            J2=0;
            parSet(1)=k; %Bulk modulus
            parSet(2)=J1; %J1
            parSet(3)=J2; %J2
        case 12 %Moerman 2A
            formulationName='Moerman 2A';
            b1=3;
            b2=2;
            parSet(1)=k; %Bulk modulus
            parSet(2)=b1; %Alpha
            parSet(3)=b2; %Beta
        case 13 %Moerman 3
            formulationName='Moerman 3';
            J1=max(J)+0.1;
            J2=0;
            s1=0.15;
            s2=0.15;
            q1=0.9;
            q2=0.9;
            parSet(1)=k; %Bulk modulus
            parSet(2)=J1; %J1
            parSet(3)=J2; %J2
            parSet(4)=s1; %s1
            parSet(5)=s2; %s2
            parSet(6)=q2; %q1
            parSet(7)=q2; %q2
    end

    %% Calculate SED

    %Get normalized SED
    [W,S,T]=SED_eval(formulationCase,parSet,J);

    %% Visualize data

    hf=cFigure;
    ht=subtitle(formulationName);
    ht.FontSize=fontSizeLabel;
    ht.Interpreter='latex';

    subplot(1,3,1); hold on;
    set(gca,'FontSize',fontSize,'LineWidth',LineWidthAxis,'GridAlpha',gridAlpha);
    xlabel('$J$','FontSize',fontSizeLabel,'Interpreter','latex');
    ylabel('$\Psi/\kappa$','FontSize',fontSizeLabel,'Interpreter','latex');

    hp1=plot(J,W,'k-','LineWidth',lineWidth);
    hp1.Color=plotColors(1,:);

    grid on; axis tight; axis square; box on;
    xlim([0 max(J(:))]);
    set(gca,'XTick',xtickRange);

    subplot(1,3,2); hold on;
    set(gca,'FontSize',fontSize,'LineWidth',LineWidthAxis,'GridAlpha',gridAlpha);
    xlabel('$J$','FontSize',fontSizeLabel,'Interpreter','Latex');
    ylabel('$\sigma_{h}/\kappa$','FontSize',fontSizeLabel,'Interpreter','Latex');
    grid on; axis tight; axis square; box on;

    hp2=plot(J,S,'k-','LineWidth',lineWidth);
    hp2.Color=plotColors(2,:);

    grid on; axis tight; axis square; box on;
    xlim([0 max(J(:))]);
    set(gca,'XTick',xtickRange);

    subplot(1,3,3); hold on;
    set(gca,'FontSize',fontSize,'LineWidth',LineWidthAxis,'GridAlpha',gridAlpha);
    xlabel('$J$','FontSize',fontSizeLabel,'Interpreter','latex');
    ylabel('$\frac{\partial^2 \Psi}{\partial J^2} /\kappa$','FontSize',fontSizeLabel,'Interpreter','Latex');

    hp3=plot(J,T,'k-','LineWidth',lineWidth);
    hp3.Color=plotColors(3,:);

    grid on; axis tight; axis square; box on;
    xlim([0 max(J(:))]);
    ylim([0 max(T(:))]);
    set(gca,'XTick',xtickRange);

    drawnow;

end

%% Evaluate SED

function [W,S,T]=SED_eval(formulationCase,parSet,J)

switch formulationCase
    case 1 %Hencky
        k=parSet(1);
        W=k/2*log(J).^2;
        S=k*log(J)./J;
        T=(k-k*log(J))./J.^2;
    case 2 %Simo
        k=parSet(1);
        W=(k/2).*(J-1).^2;
        S=(k/2).*(2*J-2);
        T=k*ones(size(J));
    case 3 %Bischoff
        k=parSet(1);
        b=parSet(2);
        W=(k./b.^2).*(cosh(b*(J-1))-1);
        S=(k./b)   .* sinh(b*(J-1));
        T=k           .* cosh(b*(J-1));
    case 4 %Modified Ogden
        k=parSet(1);
        b=parSet(2);
        W=(k./b.^2).*(J.^-b - 1 + b.*log(J));
        S=(k./b)   .*(1./J     - J.^(-b-1));
        T=(k./b)   .*(-1./J.^2  + (b+1).*J.^(-b-2));
    case 5 %Hyperfoam
        mu=parSet(1);
        a=parSet(2);
        b=parSet(3);
        k=mu.*(b+1/3);
        W=(2*mu./(a.^2)).*( 3*(J.^(a/3)-1)...
            +(1./b.*( (J.^(-a.*b))-1 )) );
        S=(1./J).*(2*mu./a).*(J.^(a/3)...
            -J.^(-a.*b) );
        T=(1./(J.^2)).*(2*mu/a).*((a./3-1).*J.^(a./3) +...
            (a.*b+1).*J.^(-a*b) );
    case 6 %Doll and Schweizerhoff
        k=parSet(1);
        a=parSet(2);
        b=parSet(3);
        W=( (k/(a+b)).*( ((1/(a+1)).*(J.^(a+1))) + ((1/(b-1)).*(J.^(-b+1)))) )...
            -(k.*(1/(a+1)).*(1/(b-1)));
        S=(k/(a+b)).*(J.^a-J.^(-b));
        T=(k/(a+b)).*(a*J.^(a-1)+b*J.^(-b-1));
    case 7
        k=parSet(1);
        beta1=parSet(2);
        k2=parSet(3);
        beta2=parSet(4);
        m=parSet(5);

        %SED
        W1=k./(2.*beta1).*exp(beta1.*   (log(J).^2))-(k/(2*beta1));
        W2=k2./(m.*beta2).*exp(beta2.*abs(log(J).^m))-(k2/(m*beta2));
        W=W1+W2;

        %Stress
        S1=(k./J).*(exp(beta1.*    log(J).^2) .*     log(J));
        S2=(k2./J).*(exp(beta2.*abs(log(J)).^m).*(abs(log(J)).^(m-1)).*sign(log(J)));
        S=S1+S2;

        %Tangent
        T1=(k./J.^2) .* exp(beta1.*    log(J) .^2).*(beta1.*2.*    log(J) .^2 - log(J) + 1);
        T2=(k2./J.^2) .* exp(beta2.*abs(log(J)).^m).*(beta2.*m.*abs(log(J)).^m - log(J) + m -1).*abs(log(J)).^(m-2);
        T=T1+T2;

    case 8 %Moerman 1
        k=parSet(1);
        b1=parSet(2);
        b2=parSet(3);
        W=(k/4)      .*( (1/b1^2).*((J.^ b1)-1).^2 + ...
            (1/b2^2).*((J.^-b2)-1).^2 );
        S=(k/2)./J   .*( (1/b1  ).*(J.^( 2*b1) - J.^b1 ) - ...
            (1/b2  ).*(J.^(-2*b2) - J.^-b2) );
        T=(k/2)./J.^2.*( ((2-1/b1)*J.^( 2*b1)-(1-1/b1)*J.^b1) + ...
            ((2+1/b2)*J.^(-2*b2)-(1+1/b2)*J.^-b2) );
    case 9 %Moerman 1A
        k=parSet(1);
        b1=parSet(2);
        b2=parSet(3);
        q=parSet(4);

        W1=(k/(2*b1^2)).*(  q).*((J.^ b1)-1).^2;
        W2=(k/(2*b2^2)).*(1-q).*((J.^-b2)-1).^2;
        W=W1+W2;

        S1= (k/b1)./J.*(  q).*(J.^( 2*b1)-J.^b1 );
        S2=-(k/b2)./J.*(1-q).*(J.^(-2*b2)-J.^-b2);
        S=S1+S2;

        T1=(k./J.^2).*(  q).*((2-1/b1)*J.^( 2*b1)-(1-1/b1)*J.^b1 );
        T2=(k./J.^2).*(1-q).*((2+1/b2)*J.^(-2*b2)-(1+1/b2)*J.^-b2);
        T=T1+T2;
    case 10 %Moerman 1B
        k=parSet(1);
        b1=parSet(2);
        b2=parSet(3);

        L1=(J>=1);
        L2=(J<1);

        W=zeros(size(J));
        W(L1)=(k/(2*b1^2)).*((J(L1).^ b1)-1).^2;
        W(L2)=(k/(2*b2^2)).*((J(L2).^-b2)-1).^2;

        S=zeros(size(J));
        S(L1)= (k/b1).*(J(L1).^( 2*b1-1)-J(L1).^( b1-1));
        S(L2)=-(k/b2).*(J(L2).^(-2*b2-1)-J(L2).^(-b2-1));

        T=zeros(size(J));
        T(L1)=(k/b1).*((2*b1-1)*J(L1).^( 2*b1-2)-(b1-1)*J(L1).^( b1-2));
        T(L2)=(k/b2).*((2*b2+1)*J(L2).^(-2*b2-2)-(b2+1)*J(L2).^(-b2-2));
    case 11 %Moerman 2
        k=parSet(1);
        J1=parSet(2);
        J2=parSet(3);

        a1=(2/pi)*(J1-1);
        a2=(2/pi)*(J2-1);

        W1=(-k*a1.^2)*log(cos((J-1)/a1));
        W2=(-k*a2.^2)*log(cos((J-1)/a2));

        W=zeros(size(W1));
        W(J>=1)=W1(J>=1);
        W(J<1)=W2(J<1);
        W=real(W);
        W(J>=J1)=inf;
        W(J<=J2)=inf;

        S1=(k.*a1).*tan((J-1)/a1);
        S2=(k.*a2).*tan((J-1)/a2);
        S=zeros(size(S1));
        S(J>=1)=S1(J>=1);
        S(J<1)=S2(J<1);
        S(J>=J1)=inf;
        S(J<=J2)=-inf;

        T1=k*sec((J-1)/a1).^2;
        T2=k*sec((J-1)/a2).^2;
        T=zeros(size(T1));
        T(J>=1)=T1(J>=1);
        T(J<1)=T2(J<1);
        T(J>=J1)=inf;
        T(J<=J2)=inf;
    case 12 %Moerman 2A
        k=parSet(1);
        b1=parSet(2);
        b2=parSet(3);

        W1=(k./b1.^2).*(cosh(b1*(J-1))-1);
        S1=(k./b1)   .* sinh(b1*(J-1));
        T1=k            .* cosh(b1*(J-1));

        W2=(k./b2.^2).*(cosh(b2*(J-1))-1);
        S2=(k./b2)   .* sinh(b2*(J-1));
        T2=k            .* cosh(b2*(J-1));

        W3=k.*(-4/pi^2) .*log(cos(pi/2*(1-J)));
        S3=k.*(-2/pi)   .*tan(pi/2*(1-J));
        T3=k            .*sec(pi/2*(1-J)).^2;

        W=W1;
        S=S1;
        T=T1;

        W(J<1)=W2(J<1)/2+W3(J<1)/2;
        S(J<1)=S2(J<1)/2+S3(J<1)/2;
        T(J<1)=T2(J<1)/2+T3(J<1)/2;
    case 13 %Moerman 3
        k=parSet(1);
        J1=parSet(2);
        J2=parSet(3);
        s_1=parSet(4);
        s_2=parSet(5);
        q1=parSet(6);
        q2=parSet(7);

        L=J<1;

        % PART 1
        a1=(pi/2)*(1/(J1-1));
        a2=(pi/2)*(1/(J2-1));

        %SED
        W11=(-1/a1.^2)*log(cos((J-1).*a1));
        W21=(-1/a2.^2)*log(cos((J-1).*a2));
        W1=zeros(size(J));
        W1(~L)=W11(~L);
        W1(L)=W21(L);
        W1=real(W1);
        W1(J>J1)=inf;
        W1(J<J2)=inf;

        %Stress
        S11=(1./a1).*tan((J-1).*a1);
        S21=(1./a2).*tan((J-1).*a2);
        S1=zeros(size(J));
        S1(~L)=S11(~L);
        S1(L)=S21(L);
        S1(J>J1)=inf;
        S1(J<J2)=-inf;

        %Tangent
        T11=sec((J-1).*a1).^2;
        T21=sec((J-1).*a2).^2;
        T1=zeros(size(J));
        T1(~L)=T11(~L);
        T1(L)=T21(L);
        T1(J>J1)=inf;
        T1(J<J2)=inf;

        % PART 2
        b1= k/(s_1);
        b2= k/(s_2);

        W12=(1/b1.^2)*log(cosh((J-1).*b1));
        W22=(1/b2.^2)*log(cosh((J-1).*b2));
        W2=zeros(size(J));
        W2(~L)=W12(~L);
        W2(L)=W22(L);
        W2=real(W2);

        S12=1/b1*tanh((J-1)*b1);
        S22=1/b2*tanh((J-1)*b2);
        S2=zeros(size(J));
        S2(~L)=S12(~L);
        S2(L)=S22(L);

        T12=sech((J-1).*b1).^2;
        T22=sech((J-1).*b2).^2;
        T2=zeros(size(J));
        T2(~L)=T12(~L);
        T2(L)=T22(L);

        % SUM

        W(~L)=k*((1-q1)*W1(~L)+q1*W2(~L));
        S(~L)=k*((1-q1)*S1(~L)+q1*S2(~L));
        T(~L)=k*((1-q1)*T1(~L)+q1*T2(~L));

        W(L)=k*((1-q2)*W1(L)+q2*W2(L));
        S(L)=k*((1-q2)*S1(L)+q2*S2(L));
        T(L)=k*((1-q2)*T1(L)+q2*T2(L));

        %%

        W(J>=J1)=inf;
        W(J<=J2)=inf;

        S(J>=J1)=inf;
        S(J<=J2)=-inf;

        T(J>=J1)=inf;
        T(J<=J2)=inf;

end

%Normalise based on bulk-modulus
W=W/k; S=S/k; T=T/k;

end

%%
%
% <<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 triangle_ncc_rule_test05 ( )

%*****************************************************************************80
%
%% TEST05 demonstrates REFERENCE_TO_PHYSICAL_T3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  node_num = 3;

  node_xy = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    0.0, 1.0 ]';
  node_xy2 = [ ...
    1.0, 2.0; ...
    1.0, 1.0; ...
    3.0, 2.0 ]';
  point_show = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  REFERENCE_TO_PHYSICAL_T3 transforms a rule\n' );
  fprintf ( 1, '  on the unit (reference) triangle to a rule on\n' );
  fprintf ( 1, '  an arbitrary (physical) triangle.\n' );

  rule = 3;

  order_num = triangle_ncc_order_num ( rule );

  [ xy, w ] = triangle_ncc_rule ( rule, order_num );
%
%  Here is the reference triangle, and its rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The reference triangle:\n', rule );
  fprintf ( 1, '\n' );

  for node = 1 : 3
    fprintf ( 1, '  %8d  %14f  %14f\n', node, node_xy(1:2,node) );
  end

  area = triangle_area ( node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rule %d for reference triangle\n', rule );
  fprintf ( 1, '  with area = %f\n', area );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                X               Y               W\n' );
  fprintf ( 1, '\n' );

  for order = 1 : order_num
    fprintf ( 1, '  %8d  %14f  %14f  %14f\n', order, xy(1:2,order), w(order) );
  end
%
%  Transform the rule.
%
  xy2 = reference_to_physical_t3 ( node_xy2, order_num, xy );
%
%  Here is the physical triangle, and its transformed rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The physical triangle:\n' );
  fprintf ( 1, '\n' );

  for node = 1 : 3
    fprintf ( 1, '  %8d  %14f  %14f\n', node, node_xy2(1:2,node) );
  end

  area2 = triangle_area ( node_xy2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rule %d for physical triangle', rule );
  fprintf ( 1, '  with area = %f\n', area2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                X               Y               W\n' );
  fprintf ( 1, '\n' );

  for order = 1 : order_num
    fprintf ( 1, '  %8d  %14f  %14f  %14f\n', order, xy2(1:2,order), w(order) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function v2 = r8vec_any_normal ( dim_num, v1 )

%*****************************************************************************80
%
%% R8VEC_ANY_NORMAL returns some normal vector to V1.
%
%  Discussion:
%
%    If DIM_NUM < 2, then no normal vector can be returned.
%
%    If V1 is the zero vector, then any unit vector will do.
%
%    No doubt, there are better, more robust algorithms.  But I will take
%    just about ANY reasonable unit vector that is normal to V1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real V1(DIM_NUM), the vector.
%
%    Output, real V2(DIM_NUM), a vector that is
%    normal to V2, and has unit Euclidean length.
%
  if ( dim_num < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_ANY_NORMAL - Fatal error!\n' );
    fprintf ( 1, '  Called with DIM_NUM < 2.\n' );
    error ( 'R8VEC_ANY_NORMAL - Fatal error!' );
  end

  if ( r8vec_norm ( dim_num, v1 ) == 0.0 )
    v2(1) = 1.0;
    v2(2:dim_num) = 0.0;
    return
  end
%
%  Seek the largest entry in V1, VJ = V1(J), and the
%  second largest, VK = V1(K).
%
%  Since V1 does not have zero norm, we are guaranteed that
%  VJ, at least, is not zero.
%
  j = -1;
  vj = 0.0;

  k = -1;
  vk = 0.0;

  for i = 1 : dim_num

    if ( abs ( vk ) < abs ( v1(i) ) || k < 1 )

      if ( abs ( vj ) < abs ( v1(i) ) || j < 1 )
        k = j;
        vk = vj;
        j = i;
        vj = v1(i);
      else
        k = i;
        vk = v1(i);
      end

    end

  end
%
%  Setting V2 to zero, except that V2(J) = -VK, and V2(K) = VJ,
%  will just about do the trick.
%
  v2(1:dim_num) = 0.0;

  v2(j) = -vk / sqrt ( vk * vk + vj * vj );
  v2(k) =  vj / sqrt ( vk * vk + vj * vj );

  return
end


text:
function move_obj(arg)
% callback function for draggable objects
% Any object can be made draggable via
%   set(obj, 'ButtonDownFcn', 'move_obj(1)');

% using deltas allows us to drag big objects
persistent last_pos;
persistent windowbuttonmotionfcn;
persistent windowbuttonupfcn;

% handle events
switch arg
  case 2    % button motion
    pos = get(gca, 'CurrentPoint');
    pos = pos(1,:);

    switch get(gco, 'type')
      case 'text', obj_pos = get(gco, 'Pos')';
      case 'line',
	obj_pos(1,:) = get(gco,'xdata');
	obj_pos(2,:) = get(gco,'ydata');
      otherwise error(['cannot handle type ' get(gco,'type')])
    end
    % if the scale is logarithmic then the delta is a ratio
    if strcmp(get(gca,'xscale'), 'log')
      new_pos(1,:) = obj_pos(1,:) * (pos(1)/last_pos(1));
    else
      new_pos(1,:) = obj_pos(1,:) + (pos(1)-last_pos(1));
    end
    if strcmp(get(gca,'yscale'), 'log')
      new_pos(2,:) = obj_pos(2,:) * (pos(2)/last_pos(2));
    else
      new_pos(2,:) = obj_pos(2,:) + (pos(2)-last_pos(2));
    end
    switch get(gco, 'type')
      case 'text', set(gco, 'Pos', new_pos);
      case 'line',
	set(gco, 'xdata', new_pos(1,:), ...
	         'ydata', new_pos(2,:));
    end
    last_pos = pos;
  case 1    % buttondown
    % start moving
    % dragging looks better with double buffering on
    set(gcf, 'DoubleBuffer', 'on');
    last_pos = get(gca,'CurrentPoint');
    last_pos = last_pos(1,:);
    % set callbacks
    fig = gcf;
    %set(fig, 'pointer', 'fleur');
    windowbuttonmotionfcn = get(fig,'windowbuttonmotionfcn');
    windowbuttonupfcn = get(fig,'windowbuttonupfcn');
    set(fig, 'windowbuttonmotionfcn', ['move_obj(2);' windowbuttonmotionfcn]);
    set(fig, 'windowbuttonupfcn', ['move_obj(3);' windowbuttonupfcn]);
  case 3    % button up
    % finished moving
    % clear callbacks
    fig = gcf;
    %set(fig, 'pointer', 'arrow');
    set(fig, 'windowbuttonmotionfcn', windowbuttonmotionfcn);
    set(fig, 'windowbuttonupfcn', windowbuttonupfcn);
  otherwise
    error('invalid argument')
end

author: andrewssobral

repo: lrslibrary


text:
function [g1, g2] = disimXrbfKernGradient(disimKern, rbfKern, t1, t2, covGrad)

% DISIMXRBFKERNGRADIENT Compute gradient between the DISIM and RBF kernels.
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between DISIM and RBF kernels for
% the multiple output kernel.
% ARG disimKern : the kernel structure associated with the DISIM
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of DISIM kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of RBF kernel.
%
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between DISIM and RBF kernels for
% the multiple output kernel.
% ARG disimKern : the kernel structure associated with the DISIM
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of DISIM kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of RBF kernel.
%
% SEEALSO : multiKernParamInit, multiKernCompute, disimKernParamInit, rbfKernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007-2009

% KERN

arg{1} = t1;
if nargin < 5
  covGrad = t2;
  t2 = t1;
else
  arg{2} = t2;
end
if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if disimKern.inverseWidth ~= rbfKern.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end
if disimKern.rbf_variance ~= rbfKern.variance
  error('Kernels cannot be cross combined if they have different RBF variances.');
end

dim1 = size(t1, 1);
dim2 = size(t2, 1);
t1Mat = t1(:, ones(1, dim2));
t2Mat = t2(:, ones(1, dim1))';
diffT = (t1Mat - t2Mat);
l = sqrt(2/disimKern.inverseWidth);
l2 = l*l;
C_0 = sqrt(disimKern.di_variance);
C_i = sqrt(disimKern.variance);
C_j = rbfKern.variance;
D_i = disimKern.decay;
delta = disimKern.di_decay;

invLDiffT = 1/l*diffT;
halfLD_i = 0.5*l*D_i;
halfLDelta = 0.5*l*delta;

prefact = C_0 * C_i * C_j * sqrt(pi)/2 * l;

lnCommon = - log(delta - D_i);
[lnPart1, sign1] = lnDiffErfs(halfLDelta - invLDiffT, halfLDelta + t2Mat/l);
[lnPart2, sign2] = lnDiffErfs(halfLD_i + t2Mat/l, halfLD_i - invLDiffT);

lnFact1 = halfLDelta^2 - delta * diffT;
lnFact2 = halfLD_i^2 - D_i * diffT;

k = sign1 .* exp(lnCommon + lnFact1 + lnPart1) ...
    + sign2 .* exp(lnCommon + lnFact2 + lnPart2);

k = 0.5*sqrt(disimKern.variance)*sqrt(disimKern.di_variance)*rbfKern.variance*k*sqrt(pi)*l;

[dlnPart1, m1] = gradLnDiffErfs(halfLDelta - invLDiffT, halfLDelta + t2Mat/l, ...
				l/2, l/2);
[dlnPart2, m2] = gradLnDiffErfs(halfLD_i + t2Mat/l, halfLD_i - invLDiffT, ...
				l/2, l/2);

dK_dD = k .* (1./(delta-D_i)) ...
	+ prefact ...
	.* ((l*halfLD_i - diffT) .* sign2 .* exp(lnCommon + lnFact2 + lnPart2) ...
	    + dlnPart2 .* exp(lnCommon + lnFact2 - m2));
dk_dD = sum(sum(dK_dD.*covGrad));

dK_ddelta = k .* (-1./(delta-D_i)) ...
    + prefact ...
    .* ((l*halfLDelta - diffT) .* sign1 .* exp(lnCommon + lnFact1 + lnPart1) ...
	+ dlnPart1 .* exp(lnCommon + lnFact1 - m1));
dk_ddelta = sum(sum(dK_ddelta.*covGrad));

[dlnPart1, m1] = gradLnDiffErfs(halfLDelta - invLDiffT, halfLDelta + t2Mat/l, ...
				delta/2 + invLDiffT/l, delta/2 - t2Mat/l2);
[dlnPart2, m2] = gradLnDiffErfs(halfLD_i + t2Mat/l, halfLD_i - invLDiffT, ...
				D_i/2 - t2Mat/l2, D_i/2 + invLDiffT/l);

dK_dl = k/l ...
	+ prefact ...
	.* (delta*halfLDelta .* sign1 .* exp(lnCommon + lnFact1 + lnPart1) ...
	    + D_i*halfLD_i .* sign2 .* exp(lnCommon + lnFact2 + lnPart2) ...
	    + dlnPart1 .* exp(lnCommon + lnFact1 - m1) ...
	    + dlnPart2 .* exp(lnCommon + lnFact2 - m2));
dk_dl = sum(sum(dK_dl.*covGrad));

dk_dC_i = sum(sum(k.*covGrad))/C_i;
dk_dC_0 = sum(sum(k.*covGrad))/C_0;
dk_dRbfVariance = sum(sum(k.*covGrad))/rbfKern.variance;

dk_dinvWidth = -0.5*sqrt(2)/(disimKern.inverseWidth* ...
                             sqrt(disimKern.inverseWidth))*dk_dl;
dk_dDisimVariance = dk_dC_i*0.5/C_i;
dk_dDIVariance = dk_dC_0*0.5/C_0;

% only pass the gradient with respect to the inverse width to one
% of the gradient vectors ... otherwise it is counted twice.
g1 = real([dk_ddelta dk_dinvWidth dk_dDIVariance dk_dD dk_dDisimVariance 0]);
g2 = real([0 dk_dRbfVariance]);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function FIM=observedFisherInfo(z,RInv,h,JacobMat,HessMat)
%%OBSERVEDFISHERINFO Assuming that a linear or nonlinear measurement is
%           corrupted with zero-mean Gaussian noise, the observed Fisher
%           information matrix (FIM) has a standard form in terms of the
%           values of the measurement function and its first and second
%           derivatives. This function takes those values and returns the
%           observed FIM. Summing the FIMs from multiple simultaneous
%           independent measurements or measurement components returns the
%           observed FIM for the fused measurement. The inverse of the FIM
%           is the Cramér-Rao lower bound (CRLB). If only a single
%           measurement is considered, and h=z, then h, z, and HessMat can
%           all be omitted. Usualy, there is no benefit to including the
%           terms.
%
%INPUTS: z The zDimX1 measurement, or if multiple measurements of the same
%          time are to the zDimXnumMeas matrix of those measurements. This
%          can be omitted if one just wants the FIM without this.
%     RInv The zDimXzDim inverse of the covariance matrix associated with
%          the multivariate Gaussian noise corrupting z, or if multiple
%          measurements are to be fused AND RInv differs among them, then a
%          zDimXzDimXnumMeas collection of all of the inverse matrices. If
%          z is omitted and multiple measurements are fused, then RInv MUST
%          be specified as a zDimXzDimXnumMeas matrix, not as a single
%          zDimXzDim matrix.
%        h The zDimX1 value of the xDimX1 state converted into the
%          measurement domain. If z is omitted, then this is not needed.
% JacobMat The zDimXxDim Jacobian matrix of derivatives of the measurement
%          function h taken with respect to the elements of the target
%          state. This is assumed the same for all measurements fused by
%          this function.
%  HessMat The xDimXxDimXzDim matrix of second derivatives of the
%          measurement function h with respect to the elements of the state
%          x. HessMat(i,j,k) is the Hessian for the kth measurement
%          component with derivatives taken with respect to elements i and
%          j of the x vector. i and j can be equal. Note that all
%          HessMat(:,:,k) are symmetric matrices. In 3D, the order of the
%          second derivatives in each submatrix is of the form:
%                  [d^2/(dxdx), d^2/(dxdy), d^2/(dxdz);
%                   d^2/(dydx), d^2/(dydy), d^2/(dydz);
%                   d^2/(dzdx), d^2/(dzdy), d^2/(dzdz)];
%
%OUTPUTS: FIM The xDimXxDim observed Fisher information matrix.
%
%The FIM and CRLB and in many statistics texts. When considering target
%tracking, one can look at Chapter 2.7.2 of [1]. Since no expectation is
%taken for the observed Fisher information matrix, only the form in terms
%of second derivatives in [1] can be used, not the form in terms of an
%outer product of first derivatives. The use of the inverse of the observed
%FIM in characterizing the accuracy of ML estimates is discussed in [2].
%
%The observed FIM is simply the negative of the matrix of second
%derivatives of the logarithm of the likelihood function. In the problem at
%hand:
%nabla_{x}(\nabla_x)'log(p(z|x))
%where here p(z|x)=1/sqrt(det(2*pi*R))*exp(-(1/2)*(z-h(x))'*inv(R)*(z-h(x))
%The gradient of the logarithm of the likelihood function is
%nabla_{x}log(p(z|x))=-H'*inv(R)(h(x)-z)
%where H=nabla_{x} h(x)'
%The matrix of second derivatives is thus
%nabla_{x}(\nabla_x)'log(p(z|x))=-H'*inv(R)*H-C
%where the jth column of C is given by
%C(:,j)=((\partial / \partial x_j)H')*inv(R)*(h(x)-z)
%
%EXAMPLE 1:
%In this example, we consider how well the observed FIM can be used as the
%covariance matrix of a fused measurement. In this instance, with all of
%the measurement being the same accuracy, we just average the measurements
%to get a fused measurement. We then find the NEES both with and without
%using the Hessian term. It is seen that when using the Hessian term, the
%NEES is closet to 1 than when not using the Hessian term.
% numMCRuns=10000;
% numMeas=3;
% sigmaR=100;
% sigmaAz=3*(pi/180);
% sigmaEl=3*(pi/180);
% SR=diag([sigmaR;sigmaAz;sigmaEl]);
% R=SR*SR';
% RInv=inv(R);
%
% zTrue=[30e3;60*(pi/180);3*(pi/180)];
% systemType=0;
% useHalfRange=true;
% xTrue=spher2Cart(zTrue,systemType,useHalfRange);
%
% NEESWithHess=0;
% NEESWithoutHess=0;
% for k=1:numMCRuns
%     zMeas=zeros(3,numMeas);
%     for curMeas=1:numMeas
%         zMeas(:,curMeas)=zTrue+SR*randn(3,1);
%     end
%      xAvg=mean(spher2Cart(zMeas,systemType,useHalfRange),2);
%
%      h=Cart2Sphere(xAvg,systemType,useHalfRange);
%      JacobMat=calcSpherJacob(xAvg,systemType,useHalfRange);
%      HessMat=calcSpherHessian(xAvg,systemType,useHalfRange);
%
%      FIMWithHess=observedFisherInfo(zMeas,RInv,h,JacobMat,HessMat);
%      FIMWithoutHess=numMeas*observedFisherInfo([],RInv,[],JacobMat,HessMat);
%      diff=xAvg-xTrue;
%      NEESWithHess=NEESWithHess+diff'*FIMWithHess*diff;
%      NEESWithoutHess=NEESWithoutHess+diff'*FIMWithoutHess*diff;
% end
% NEESWithHess=NEESWithHess/(3*numMCRuns)
% NEESWithoutHess=NEESWithoutHess/(3*numMCRuns)
%
%EXAMPLE 2:
%In this instance, we used the observed Fisher information as a covariance
%matrix of a single spherical measurement in the absence of knowing the
%truth. Thus, we just take h(z)=z and can omit the matrix of second
%derivatives. Evaluating the NEES, one sees it is consistent (near 1, maybe
%like 0.99 or 1.001). Of course, at higher angular noise levels, a debiased
%function like spher2CartTaylor can perform better.
% numMCRuns=10000;
% sigmaR=10;
% sigmaAz=0.1*(pi/180);
% sigmaEl=0.1*(pi/180);
% SR=diag([sigmaR;sigmaAz;sigmaEl]);
% R=SR*SR';
% RInv=inv(R);
%
% zTrue=[1e3;60*(pi/180);3*(pi/180)];
% systemType=0;
% useHalfRange=true;
% xTrue=spher2Cart(zTrue,systemType,useHalfRange);
%
% NEES=0;
% for k=1:numMCRuns
%     zMeas=zTrue+SR*randn(3,1);
%     xConv=spher2Cart(zMeas,systemType,useHalfRange);
%     JacobMat=calcSpherJacob(xConv,systemType,useHalfRange);
%
%     invCRLB=observedFisherInfo([],RInv,[],JacobMat);
%
%     diff=xConv-xTrue;
%     NEES=NEES+diff'*invCRLB*diff;
% end
% NEES=NEES/(3*numMCRuns)
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation: Theory, Algorithms and
%    Software. New York: John Wiley and Sons, 2001.
%[2] B. Efron and D. Hinkley, "Assessing the accuracy of the maximum
%    likelihood estimator: Observed versus expected Fisher information,"
%    Department of Statistics, Stanford, University, Tech. Rep. 108, 8 Mar.
%    1978.
%
%December 2020 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~isempty(z))
    numZ=size(z,2);

    xDim=size(HessMat,1);
    zDim=size(HessMat,3);
    FIM=zeros(xDim,xDim);

    C=zeros(xDim,xDim);
    for curMeas=1:numZ
        if(size(RInv,3)>1)
            RInvCur=RInv(:,:,curMeas);
        else
            RInvCur=RInv;
        end
        FIM=FIM+JacobMat'*RInv*JacobMat;

        RhzVal=RInvCur*(h-z(:,curMeas));
        for k=1:xDim
            C(:,k)=reshape(HessMat(:,k,:),[xDim,zDim])*RhzVal;
        end
        FIM=FIM+C;
    end
else
    numMeas=size(RInv,3);
    xDim=size(JacobMat,2);
    FIM=zeros(xDim,xDim);
    for k=1:numMeas
        FIM=FIM+JacobMat'*RInv(:,:,k)*JacobMat;
    end
end
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 stroud_test37 ( )

%*****************************************************************************80
%
%% TEST37 tests SQUARE_UNIT_SET, SQUARE_UNIT_SUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_2D_INDEX;

  max_rule = 6;

  num = function_2d_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST37\n' );
  fprintf ( 1, '  SQUARE_UNIT_SET sets up quadrature on the unit square;\n' );
  fprintf ( 1, '  SQUARE_UNIT_SUM carries it out on the unit square.\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : max_rule

    ihi = min ( ilo + 4, max_rule );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Rule:    ' );
    for rule = ilo : ihi
      fprintf ( 1, '  %6d     ', rule );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Function\n' );
    fprintf ( 1, '\n' );

    for i = 1 : num

      FUNC_2D_INDEX = i;

      for rule = ilo : ihi

        order = square_unit_size ( rule );

        [ xtab, ytab, weight ] = square_unit_set ( rule, order );

        result(rule) = square_unit_sum ( 'function_2d', order, xtab, ytab, weight );

      end

      fname = function_2d_name ( i );

      fprintf ( 1, '  %s', fname );
      for rule = ilo : ihi
        fprintf ( 1, '  %11f', result(rule) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function extras = write_extras(fname,extras)
% Write extra bits of information
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: write_extras.m 1143 2008-02-07 19:33:33Z spm $

[pth,nam,ext] = fileparts(fname);
switch ext
case {'.hdr','.img','.nii'}
    mname = fullfile(pth,[nam '.mat']);
case {'.HDR','.IMG','.NII'}
    mname = fullfile(pth,[nam '.MAT']);
otherwise
    mname = fullfile(pth,[nam '.mat']);
end
if isstruct(extras) && ~isempty(fieldnames(extras)),
    savefields(mname,extras);
end;

function savefields(fnam,p)
if length(p)>1, error('Can''t save fields.'); end;
fn = fieldnames(p);
for i_=1:length(fn),
    eval([fn{i_} '= p.' fn{i_} ';']);
end;
if str2num(version('-release'))>=14,
    fn = {'-V6',fn{:}};
end;
if numel(fn)>0,
    save(fnam,fn{:});
end;
return;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% This code calculates the 3D distances between all possible pairs
% (combination of n epicenters taken 2 at a time) of earthquakes of
% a given dataset.
%
%
% Attributing the corresponding catalog to E
%
if ~exist('index', 'var')
    index = 1;
end

if index == 1
    E = newt2;
elseif index == 2
    E = ran;
elseif index == 3
    E = rann;
end
%
% Variables
%
N = size(E,1);				% N= # of events in the catalogue; E= Earthquake catalogue
pairdist = []; 			% pairdist= Vector of interevent distances
j = nchoosek(N,2)			% j= # of interevent distances calculated
pairdist = zeros(j,1);
depth = zeros(j,1);
k = 0;

Ho_Wb = waitbar(0,'Calculating the fractal dimension');
Hf_Cfig = gcf;
Hf_child = get(groot,'children');
set(Hf_child,'pointer','watch','papertype','A4');
%
% Calculation of the interevent distances in 2D plus the depths differences.
%
for i = 1:(N-1)

    lon1 = repmat(E(i,1), [(N-i),1]);
    lat1 = repmat(E(i,2), [(N-i),1]);
    depth1 = repmat(E(i,7), [(N-i),1]);
    lon2 = E((i+1):end, 1);
    lat2 = E((i+1):end, 2);
    depth2 = E((i+1):end, 7);
    pairdist(k+1:k + size(lon1, 1)) = distance(lat1,lon1,lat2,lon2);
    depth(k+1:k + size(lon1, 1)) = depth1-depth2;
    k = k + size(lon1,1);
    waitbar((0.5/(N-1))*i, Ho_Wb);

end
%
% Converts the interevent distances from degrees to kilometers and calculates
% the interevent distances in three dimensions.
%
pairdist = pairdist.*111;
pairdist = (pairdist.^2 + depth.^2).^0.5;
clear depth;
%
% Compute the correlation integral
%
d = 3;			%the embedding dimension
docorint;

author: CelsoReyes


text:
function [meanDistMatrix, maxDistMatrix] = mtrComputeFiberDistMatrix(fg)

% Find all pair differences between fibers in the fiber group
meanDistMatrix = zeros(length(fg.fibers));
maxDistMatrix = zeros(length(fg.fibers));

disp('Computing distances between all paths could take awhile.');
percentMark = 10;
percentInc = 10;
for ii = 1:length(fg.fibers)-1
    for jj = ii+1:length(fg.fibers)
        [indices, bestSqDistIJ] = nearpoints(fg.fibers{ii}, fg.fibers{jj});
        [indices, bestSqDistJI] = nearpoints(fg.fibers{jj}, fg.fibers{ii});
        meanDistMatrix(ii,jj) = mean([mean(sqrt(bestSqDistIJ)) mean(sqrt(bestSqDistJI))]);
        meanDistMatrix(jj,ii) = meanDistMatrix(ii,jj);
        maxDistMatrix(ii,jj) = mean(sqrt([max(bestSqDistIJ) max(bestSqDistJI)]));
        maxDistMatrix(jj,ii) = maxDistMatrix(ii,jj);
    end
    percentCompleted = round(ii/length(fg.fibers)*100);
    if( percentCompleted >= percentMark )
        disp(['Completed ' num2str(percentCompleted) ' %' ]);
        percentMark = percentMark+percentInc;
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function e2s2_storeSPGTOverlap(varargin)
% e2s2_storeSPGTOverlap(varargin)
%
% Augment the superpixel structure with the overlap with the GT blob.
%
% Copyright by Holger Caesar, 2015

% Parse input
p = inputParser;
addParameter(p, 'dataset', SiftFlowDatasetMC());
addParameter(p, 'projectName', 'WeaklySupervisedLearning');
addParameter(p, 'spName', 'Felzenszwalb2004-k100-sigma0.8-colorTypesRgb');
addParameter(p, 'gtName', 'GroundTruth');
parse(p, varargin{:});

dataset = p.Results.dataset;
projectName = p.Results.projectName;
spName = p.Results.spName;
gtName = p.Results.gtName;

% Create paths
global glFeaturesFolder;
segmentFolderSP = fullfile(glFeaturesFolder, projectName, dataset.name, 'segmentations', spName);
segmentFolderGT = fullfile(glFeaturesFolder, projectName, dataset.name, 'segmentations', gtName);

% Get image list
[imageList, imageCount] = dataset.getImageList(true);

for imageIdx = 1 : imageCount,
    printProgress(sprintf('Storing GT-SP (%s) overlap for image', spName), imageIdx, imageCount, 50);
    imageName = imageList{imageIdx};
    imageSize = dataset.getImageSize(imageName);

    % Get SP blobs
    segmentPathSP = fullfile(segmentFolderSP, [imageName, '.mat']);
    assert(exist(segmentPathSP, 'file') ~= 0);
    segmentStructSP = load(segmentPathSP, 'propBlobs');
    blobsSP = segmentStructSP.propBlobs;

    % Get GT blobs
    segmentPathGT = fullfile(segmentFolderGT, [imageName, '.mat']);
    assert(exist(segmentPathGT, 'file') ~= 0);
    segmentStructGT = load(segmentPathGT, 'propBlobs');
    blobsGT = segmentStructGT.propBlobs;

    % Compute overlap
    overlapRatiosSPGT = computeBlobOverlapSum(blobsSP, blobsGT, imageSize); %#ok<NASGU>

    % Store to disk (append)
    assert(exist(segmentPathSP, 'file') ~= 0);
    save(segmentPathSP, 'overlapRatiosSPGT', '-append');
end;

author: nightrome

repo: matconvnet-calvin


text:
function [data_rec] = gtmmapimputate(net, data)

%GTMEXPIMPUTE Impute missing data using GTM MAP estimates
%
%	Description
%	 DATA_REC = GTMEXPIMPUTE(NET, DATA) takes a GTM structure NET, and
%	imputes the missing values in DATA according to MAP estimates of hidden
% variables given observed data.
%
%	See also
%	GTMEXPIMPUTE, GTM, GTMEM, GTMLMEAN, GMLMODE, GMMPROB

% Copyright (c) Tommi Vatanen (2012)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

data_rec = data;
missing = isnan(data_rec);

%
net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
R = gmmpost(net.gmmnet, data);

R_tmp = zeros(size(R));
R_max = repmat(max(R,[],2), [1 size(R,2)]);
R_tmp(R==R_max) = 1;

% pick up map reference vectors
rec = R_tmp*net.gmmnet.centres;
data_rec(missing) = rec(missing);

author: ilarinieminen

repo: SOM-Toolbox


text:
function exercise4(varargin)
% EXERCISE4   Part 4 of the VGG CNN practical

setup ;

% -------------------------------------------------------------------------
% Part 4.1: prepare the data
% -------------------------------------------------------------------------

% Load character dataset
imdb = load('data/charsdb.mat') ;

% Visualize some of the data
figure(10) ; clf ; colormap gray ;
subplot(1,2,1) ;
vl_imarraysc(imdb.images.data(:,:,imdb.images.label==1 & imdb.images.set==1)) ;
axis image off ;
title('training chars for ''a''') ;

subplot(1,2,2) ;
vl_imarraysc(imdb.images.data(:,:,imdb.images.label==1 & imdb.images.set==2)) ;
axis image off ;
title('validation chars for ''a''') ;

% -------------------------------------------------------------------------
% Part 4.2: initialize a CNN architecture
% -------------------------------------------------------------------------

net = initializeCharacterCNN() ;

% -------------------------------------------------------------------------
% Part 4.3: train and evaluate the CNN
% -------------------------------------------------------------------------

trainOpts.batchSize = 100 ;
trainOpts.numEpochs = 15 ;
trainOpts.continue = true ;
trainOpts.gpus = [] ;
trainOpts.learningRate = 0.001 ;
trainOpts.expDir = 'data/chars-experiment' ;
trainOpts = vl_argparse(trainOpts, varargin);

% Take the average image out
imdb = load('data/charsdb.mat') ;
imageMean = mean(imdb.images.data(:)) ;
imdb.images.data = imdb.images.data - imageMean ;

% Convert to a GPU array if needed
if numel(trainOpts.gpus) > 0
  imdb.images.data = gpuArray(imdb.images.data) ;
end

% Call training function in MatConvNet
[net,info] = cnn_train(net, imdb, @getBatch, trainOpts) ;

% Move the CNN back to the CPU if it was trained on the GPU
if numel(trainOpts.gpus) > 0
  net = vl_simplenn_move(net, 'cpu') ;
end

% Save the result for later use
net.layers(end) = [] ;
net.imageMean = imageMean ;
save('data/chars-experiment/charscnn.mat', '-struct', 'net') ;

% -------------------------------------------------------------------------
% Part 4.4: visualize the learned filters
% -------------------------------------------------------------------------

figure(2) ; clf ; colormap gray ;
vl_imarraysc(squeeze(net.layers{1}.weights{1}),'spacing',2)
axis equal ; title('filters in the first layer') ;

% -------------------------------------------------------------------------
% Part 4.5: apply the model
% -------------------------------------------------------------------------

% Load the CNN learned before
net = load('data/chars-experiment/charscnn.mat') ;
%net = load('data/chars-experiment/charscnn-jit.mat') ;

% Load the sentence
[im,cmap] = imread('data/sentence-lato.png') ;
if isempty(cmap)
  im = im2single(im) ;
else
  im = im2single(ind2gray(p,cmap)) ;
end
im = 256 * (im - net.imageMean) ;

% Apply the CNN to the larger image
res = vl_simplenn(net, im) ;

% Visualize the results
figure(3) ; clf ;
decodeCharacters(net, imdb, im, res) ;

% -------------------------------------------------------------------------
% Part 4.6: train with jitter
% -------------------------------------------------------------------------

trainOpts.batchSize = 100 ;
trainOpts.numEpochs = 15 ;
trainOpts.continue = true ;
trainOpts.learningRate = 0.001 ;
trainOpts.expDir = 'data/chars-jit-experiment' ;

% Initlialize a new network
net = initializeCharacterCNN() ;

% Call training function in MatConvNet
[net,info] = cnn_train(net, imdb, @getBatchWithJitter, trainOpts) ;

% Move the CNN back to CPU if it was trained on GPU
if numel(trainOpts.gpus) > 0
  net = vl_simplenn_move(net, 'cpu') ;
end

% Save the result for later use
net.layers(end) = [] ;
net.imageMean = imageMean ;
save('data/chars-experiment/charscnn-jit.mat', '-struct', 'net') ;

% Visualize the results on the sentence
figure(4) ; clf ;
decodeCharacters(net, imdb, im, vl_simplenn(net, im)) ;

% --------------------------------------------------------------------
function [im, labels] = getBatch(imdb, batch)
% --------------------------------------------------------------------
im = imdb.images.data(:,:,batch) ;
im = 256 * reshape(im, 32, 32, 1, []) ;
labels = imdb.images.label(1,batch) ;

% --------------------------------------------------------------------
function [im, labels] = getBatchWithJitter(imdb, batch)
% --------------------------------------------------------------------
im = imdb.images.data(:,:,batch) ;
labels = imdb.images.label(1,batch) ;

n = numel(batch) ;
train = find(imdb.images.set == 1) ;

sel = randperm(numel(train), n) ;
im1 = imdb.images.data(:,:,sel) ;

sel = randperm(numel(train), n) ;
im2 = imdb.images.data(:,:,sel) ;

ctx = [im1 im2] ;
ctx(:,17:48,:) = min(ctx(:,17:48,:), im) ;

dx = randi(11) - 6 ;
im = ctx(:,(17:48)+dx,:) ;
sx = (17:48) + dx ;

dy = randi(5) - 2 ;
sy = max(1, min(32, (1:32) + dy)) ;

im = ctx(sy,sx,:) ;

% Visualize the batch:
% figure(100) ; clf ;
% vl_imarraysc(im) ;

im = 256 * reshape(im, 32, 32, 1, []) ;

author: vedaldi

repo: practical-cnn


text:
function [h,hg,htick]=terplot(number)
%FUNCTION [h,hg,htick]=TERPLOT perpares a ternary axis system that is
% needed for the ternaryc function. It returns three handels:
% - h:      to modify the patch created by the fill function;
% - hg:     to change each grid line separately (must probably be
% modified);
% - htick:  to edit the tick labels (probably very inconvinient)
%
% Uli Theune, Geophysics, University of Alberta
%
if nargin<1
    number=11;
end
h=fill([0 1 0.5 0],[0 0 0.866 0],'w','linewidth',2);
%set(h,'facecolor',[0.7 0.7 0.7],'edgecolor','w')
%set(gcf,'color',[0 0 0.3])
d1=cos(pi/3);
d2=sin(pi/3);
l=linspace(0,1,number);
hold on
for i=2:length(l)-1
   hg(i-1,3)=plot([l(i)*d1 1-l(i)*d1],[l(i)*d2 l(i)*d2],':k','linewidth',0.25);
   hg(i-1,1)=plot([l(i) l(i)+(1-l(i))*d1],[0 (1-l(i))*d2],':k','linewidth',0.25);
   hg(i-1,2)=plot([(1-l(i))*d1 1-l(i)],[(1-l(i))*d2 0],':k','linewidth',0.25);
end
hold off
axis image
axis off
% Make x-tick labels
for i=1:number
    htick(i,1)=text(l(i),-0.025,num2str(l(i)));
    htick(i,3)=text(1-l(i)*cos(pi/3)+0.025,l(i)*sin(pi/3)+0.025,num2str(l(i)));
    htick(i,2)=text(0.5-l(i)*cos(pi/3)-0.06,sin(pi/3)*(1-l(i)),num2str(l(i)));
end

set(gcf,'WindowButtonDownFcn','InitTerExpl');
global hx;

author: Sable

repo: mcbench-benchmarks


text:
% pop_expica() - export ICA weights or inverse matrix
%
% Usage:
%   >> pop_expica( EEG, whichica);             % a window pops up
%   >> pop_expica( EEG, whichica, filename );
%
% Inputs:
%   EEG         - EEGLAB dataset
%   whichica    - ['weights'|'inv'] export ica 'weights' or ica inverse
%                 matrix ('inv'). Note: for 'weights', the function
%                 export the product of the sphere and weights matrix.
%   filename    - text file name
%
% Author: Arnaud Delorme, CNL / Salk Institute, Mai 14, 2003
%
% See also: pop_export()

% Copyright (C) Mai 14, 2003, Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function com = pop_expica(EEG, whichica, filename);

com = '';
if nargin < 1
    help pop_expica;
    return;
end;

if nargin < 2
    whichica = 'weights';
end;
switch lower(whichica)
 case {'weights' 'inv'}, ;
 otherwise error('Unrecognized option for ''whichica'' parameter');
end;

if nargin < 3
	% ask user
	[filename, filepath] = uiputfile('*.*', [ 'File name for ' ...
                        fastif(strcmpi(whichica, 'inv'), 'inverse', 'weight') ' matrix -- pop_expica()']);
    drawnow;
	if filename == 0 return; end;
	filename = [filepath filename];
end;

% save datas
% ----------
if strcmpi(whichica, 'inv')
    tmpmat = double(EEG.icawinv);
else
    tmpmat = double(EEG.icaweights*EEG.icasphere);
end;
save(filename, '-ascii', 'tmpmat');

com = sprintf('pop_expica(%s, ''%s'', ''%s'');', inputname(1), whichica, filename);

return;

author: PatternRecognition

repo: OpenBMI


text:
function varargout=nmf_classic(v,winit,hinit,peval,verbose)
%
% Jean-Philippe Brunet
% Cancer Genomics
% The Broad Institute
% brunet@broad.mit.edu
%
% This software and its documentation are copyright 2004 by the
% Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
% This software is supplied without any warranty or guaranteed support whatsoever.
% Neither the Broad Institute nor MIT can not be responsible for its use, misuse,
% or functionality.
%
% NMF divergence update equations :
% Lee, D..D., and Seung, H.S., (2001), 'Algorithms for Non-negative Matrix
% Factorization', Adv. Neural Info. Proc. Syst. 13, 556-562.
%
% v (n,m) : N (genes) x M (samples) original matrix
%           Numerical data only.
%           Must be non negative.
%           Not all entries in a row can be 0. If so, add a small constant to the
%           matrix, eg.v+0.01*min(min(v)),and restart.
%
% r       : number of desired factors (rank of the factorization)
%
% verbose : prints iteration count and changes in connectivity matrix elements
%           unless verbose is 0
%
% Note : NMF iterations stop when connectivity matrix has not changed
%        for 10*stopconv interations. This is experimental and can be
%        adjusted.
%
% w    : N x r NMF factor
% h    : r x M NMF factor
%
% winit - initial value for w
% hinit - initial value for h
% peval.h_fixvec & peval.w_fixvec:
% fixvec - vector which component should be fixed: eg [2 3] will
% fix second and third component while varying the first...
fprintf('Classic NMF iterations\n')

if ~isfield(peval, 'ddterm'); peval.ddterm = 1; end %termination criterion
if ~isfield(peval, 'maxiter'); peval.maxiter = 1000; end

% test for negative values in v
if min(min(v)) < 0
    error('matrix entries can not be negative');
    return
end
if min(sum(v,2)) == 0
    error('not all entries in a row can be zero');
    return
end

[n,m]=size(v);

if ~isempty(peval.w_fixvec)
    fprintf('Fixing [ ');
    fprintf('%g ', peval.w_fixvec);
    fprintf('] component of ''w''\n');
end

if ~isempty(peval.h_fixvec)
    fprintf('Fixing [ ');
    fprintf('%g ', peval.h_fixvec);
    fprintf('] component of ''h''\n');
end

w = winit;
h = hinit;

d(1) = ddivergence(v, w*h);
htrace(1,:,:) = hinit;

w_vec = 1:size(w,2);
peval.w_dovec = find(~(w_vec==peval.w_fixvec));

h_vec = 1:size(h,1);
peval.h_dovec = find(~(h_vec==peval.h_fixvec));

for ii=2:peval.maxiter

    w_old = w;
    h_old = h;

    sumH_t = sum(h(peval.h_dovec,:),2);
    sumH_t_sq = sum(sumH_t.^2);
    sumH= sum(sumH_t);
    nh = length(peval.h_dovec);
    sparsity_h = (sqrt(nh) - sumH/sqrt(sumH_t_sq))/(sqrt(nh)-1); % [Hoyer 2004]

    gradspars=(1/(sqrt(nh)-1))*(sumH/(sumH_t_sq)^1.5*sumH_t - 1/(sqrt(sumH_t_sq)));
    gradspars_kt = repmat(gradspars, 1, m);

    x1=repmat(sum(w,1)',1,m);
    y1=w'*(v./(w*h));
    h(peval.h_dovec,:)=h(peval.h_dovec,:).*(y1(peval.h_dovec,:))./x1(peval.h_dovec,:)+peval.alpha_sparsity*gradspars_kt;
    h=max(h,eps); % adjust small values to avoid undeflow

    x2=repmat(sum(h,2)',n,1);
    y2=(v./(w*h))*h';
    w(:,peval.w_dovec)=w(:,peval.w_dovec).*(y2(:,peval.w_dovec))./x2(:,peval.w_dovec);
    w=max(w,eps); % adjust small values to avoid undeflow

    % normalization of all h:
    sumw = sum(w,1);
    w = w./repmat(sumw,n,1); %normalization of each component
%     h = h.*repmat(sumw',1,m); %to keep the multiplication equal

    %     wtrace(:,:,ii) = w;
    %     htrace(:,:,ii) = h;

    d(ii) = ddivergence(v, w*h) + peval.alpha_sparsity*sparsity_h;
    dd(ii) = abs(d(ii)-d(ii-1));
    htrace(ii,:,:)=h;
    if dd(ii) < peval.ddterm
        break
    end
    if verbose
        fprintf('Cycle %g D-divergence %g\n',ii-1,d(ii))
    end
end

peval.numiter = ii;
peval.maxiter_reached_flag = 0;
if ii == peval.maxiter
    fprintf('\nMAximum number of iteration (%g) reached! \n', peval.maxiter)
    peval.maxiter_reached_flag = 1;
end

varargout{1}=w;
varargout{2}=h;
varargout{3}=peval;
varargout{4}=d;
varargout{5}=htrace;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [pav_trans,score_bounds,llr_bounds] = pav_calibration(tar,non,small_val)
% Creates a calibration transformation function using the PAV algorithm.
% Inputs:
%   tar: A vector of target scores.
%   non: A vector of non-target scores.
%   small_val: An offset to make the transformation function
%     invertible.  small_val is subtracted from the left-hand side
%     of the bin and added to the right-hand side (and the bin
%     height is linear between its left and right ends).
% Outputs:
%   pav_trans: The transformation function.  It takes in scores and
%     outputs (calibrated) log-likelihood ratios.
%   score_bounds: The left and right ends of the line segments
%     that make up the transformation.
%   llr_bounds: The lower and upper ends of the line segments that
%     make up the transformation.

if nargin==0
    test_this()
    return
else
    assert(nargin==3)
    assert(size(tar,1)==1)
    assert(size(non,1)==1)
    assert(length(tar)>0)
    assert(length(non)>0)
end

largeval = 1e6;

scores = [-largeval tar non largeval];
Pideal = [ones(1,length(tar)+1),zeros(1,length(non)+1)];
[scores,perturb] = sort(scores);
Pideal = Pideal(perturb);
[Popt,width,height] = pavx(Pideal);
data_prior = (length(tar)+1)/length(Pideal);
llr = logit(Popt) - logit(data_prior);
bnd_ndx = make_bnd_ndx(width);
score_bounds = scores(bnd_ndx);
llr_bounds = llr(bnd_ndx);
llr_bounds(1:2:end) = llr_bounds(1:2:end) - small_val;
llr_bounds(2:2:end) = llr_bounds(2:2:end) + small_val;
pav_trans = @(s) pav_transform(s,score_bounds,llr_bounds);
end

function scr_out = pav_transform(scr_in,score_bounds,llr_bounds)
scr_out = zeros(1,length(scr_in));
for ii=1:length(scr_in)
    x = scr_in(ii);
    [x1,x2,v1,v2] = get_line_segment_vals(x,score_bounds,llr_bounds);
    scr_out(ii) = (v2 - v1) * (x - x1) / (x2 - x1) + v1;
end
end

function bnd_ndx = make_bnd_ndx(width)
len = length(width)*2;
c = cumsum(width);
bnd_ndx = zeros(1,len);
bnd_ndx(1:2:len) = [1 c(1:end-1)+1];
bnd_ndx(2:2:len) = c;
end

function [x1,x2,v1,v2] = get_line_segment_vals(x,score_bounds,llr_bounds)
p = find(x>=score_bounds,1,'last');
x1 = score_bounds(p);
x2 = score_bounds(p+1);
v1 = llr_bounds(p);
v2 = llr_bounds(p+1);
end

function test_this()
ntar = 10;
nnon = 12;
tar = 2*randn(1,ntar)+2;
non = 2*randn(1,nnon)-2;
tarnon = [tar non];

scores = [-inf tarnon inf];
Pideal = [ones(1,length(tar)+1),zeros(1,length(non)+1)];
[scores,perturb] = sort(scores);
Pideal = Pideal(perturb);
[Popt,width,height] = pavx(Pideal);
data_prior = (length(tar)+1)/length(Pideal);
llr = logit(Popt) - logit(data_prior);
[dummy,pinv] = sort(perturb);
tmp = llr(pinv);
llr = tmp(2:end-1)

pav_trans = pav_calibration(tar,non,0);
pav_trans(tarnon)

end

author: nesl


text:
function make_ref_obj(input_var)

if isa(input_var, 'ref_obj')
    return;
end

one_mutable=ref_obj;
one_mutable.ref=input_var;
assignin('caller', inputname(1), one_mutable);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function varargout = log(varargin)
%LOG   Natural logarithm of a CHEBFUN2.
%   LOG(F) is the natural logarithm of F. This function returns an error
%   if the function passes through or becomes numerically close to zero.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Copyright (C) 2016 Abhinav Tripathi
%% Copyright (C) 2016, 2017, 2022 Colin B. Macdonald
%% Copyright (C) 2020 Tianyu Chen (billchenchina)
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

function obj = check_and_convert(var_pyobj)
  persistent builtins
  persistent sp
  persistent list_or_tuple
  persistent _sym
  if isempty(builtins)
    builtins = pyeval("__builtins__");
    list_or_tuple = py.tuple({builtins.list, builtins.tuple});

    sp = py.sympy;
    _sym = py.tuple({sp.Basic, sp.MatrixBase});
  end

  if (~ py.isinstance(var_pyobj, list_or_tuple))
    var_pyobj = {var_pyobj};
  end

  obj = {};
  for i = 1:length(var_pyobj)
    x = var_pyobj{i};

    if (~ isa (x, 'pyobject'))
      obj{i} = x;
    elseif (py.isinstance(x, _sym) || isequal(x, py.None))
      obj{i} = get_sym_from_python(x);
    elseif (py.isinstance(x, list_or_tuple))
      obj{i} = check_and_convert(x);
    elseif (isa (x, 'py.str'))
      obj{i} = char (x);
    elseif (py.isinstance(x, builtins.dict))
      make_str_keys = pyeval ('lambda x: {str(k): v for k, v in x.items()}');
      x = pycall (make_str_keys, x);
      s = struct (x);
      % make sure values are converted to sym
      s = structfun (@(t) check_and_convert (t){:}, s, 'UniformOutput', false);
      obj{i} = s;
    elseif (isa (x, 'py.int'))
      if (py.isinstance(x, pyeval('bool')))
        error ('unexpected python bool')
      end
      if (abs (double (x)) > intmax ('int64'))
        error ('precision would be lost converting integer larger than %ld', ...
               intmax ('int64'))
      end
      obj{i} = int64 (x);
    else
      warning ('OctSymPy:pythonic_no_convert',
               sprintf ('something was not converted from pyobject: %s',
                       char (x)))
      obj{i} = x;
    end
  end
end

repo: octsympy


text:
function p=mysubset(small,large)
% MYSUBSET Is the small set of +ve integers a subset of the large set?
% p = mysubset(small, large)

% Surprisingly, this is not built-in.

if isempty(small)
  p = 1; % isempty(large);
else
  p = length(myintersect(small,large)) == length(small);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
 function hn=arrow3(p1,p2,s,w,h,ip,alpha,beta)
% ARROW3 (R13)
%   ARROW3(P1,P2) draws lines from P1 to P2 with directional arrowheads.
%   P1 and P2 are either nx2 or nx3 matrices.  Each row of P1 is an
%   initial point, and each row of P2 is a terminal point.
%
%   ARROW3(P1,P2,S,W,H,IP,ALPHA,BETA) can be used to specify properties
%   of the line, initial point marker, and arrowhead.  S is a character
%   string made with one element from any or all of the following 3
%   columns:
%
%     Color Switches      LineStyle            LineWidth
%     ------------------  -------------------  --------------------
%     k  blacK (default)  -  solid (default)   0.5 points (default)
%     y  Yellow           :  dotted            0   no lines
%     m  Magenta          -. dashdot           /   LineWidthOrder
%     c  Cyan             -- dashed
%     r  Red              *  LineStyleOrder            _______ __
%     g  Green                                       ^        |
%     b  Blue                                       / \       |
%     w  White                        Arrowhead    /   \   Height
%     a  Asparagus                                /     \     |
%     d  Dark gray                               /       \    |
%     e  Evergreen                              /___   ___\ __|__
%     f  Firebrick                             |    | |    |
%     h  Hot pink                              |-- Width --|
%     i  Indigo                                |    | |    |
%     j  Jade                                       | |
%     l  Light gray                                 | |
%     n  Nutbrown                                   | |
%     p  Pear                                       | |
%     q  kumQuat                      Line       -->| |<--LineWidth
%     s  Sky blue                                   | |
%     t  Tawny                                      | |
%     u  bUrgundy                                   | |
%     v  Violet                                     | |
%     z  aZure                                      | |
%     x  random                       Initial      /   \
%     o  colorOrder                   Point    -->(     )<--IP
%     |  magnitude                    Marker       \_ _/
%
%     -------------------------------------------------------------
%                          Color Equivalencies
%     -------------------------------------------------------------
%     ColorOrder     Arrow3         |     Simulink       Arrow3
%     ----------     ----------     |     ----------     ----------
%     Color1         Blue           |     LightBlue      aZure
%     Color2         Evergreen      |     DarkGreen      Asparagus
%     Color3         Red            |     Orange         kumQuat
%     Color4         Sky blue       |     Gray           Light gray
%     Color5         Violet         |
%     Color6         Pear           |
%     Color7         Dark gray      |
%     -------------------------------------------------------------
%
%   The components of S may be specified in any order.  Invalid
%   characters in S will be ignored and replaced by default settings.
%
%   Prefixing the color code with '_' produces a darker shade, e.g.
%   '_t' is dark tawny; prefixing the color code with '^' produces a
%   lighter shade, e.g. '^q' is light kumquat.  The relative brightness
%   of light and dark color shades is controlled by the scalar parameter
%   BETA.  Color code prefixes do not affect black (k), white (w), or
%   the special color switches (xo|).
%
%   ColorOrder may be achieved in two fashions:  The user may either
%   set the ColorOrder property (using RGB triples) or define the
%   global variable ColorOrder (using a string of valid color codes).
%   If the color switch is specified with 'o', and the global variable
%   ColorOrder is a string of color codes (color switches less 'xo|',
%   optionally prefixed with '_' or '^'), then the ColorOrder property
%   will be set to the sequence of colors indicated by the ColorOrder
%   variable.  The color sequence 'bersvpd' matches the default
%   ColorOrder property.  If the color switch is specified with 'o', and
%   the global variable ColorOrder is empty or invalid, then the current
%   ColorOrder property will be used.  Note that the ColorOrder variable
%   takes precedence over the ColorOrder property.
%
%   The magnitude color switch is used to visualize vector magnitudes
%   in conjunction with a colorbar.  If the color switch is specified
%   with '|', colors are linearly interpolated from the current ColorMap
%   according to the length of the associated line.  This option sets
%   CLim to [MinM,MaxM], where MinM and MaxM are the minimum and maximum
%   magnitudes, respectively.
%
%   The current LineStyleOrder property will be used if LineStyle is
%   specified with '*'.  MATLAB cycles through the line styles defined
%   by the LineStyleOrder property only after using all colors defined
%   by the ColorOrder property.  If however, the global variable
%   LineWidthOrder is defined, and LineWidth is specified with '/',
%   then each line will be drawn with sequential color, linestyle, and
%   linewidth.
%
%   W (default = 1) is a vector of arrowhead widths; use W = 0 for no
%   arrowheads.  H (default = 3W) is a vector of arrowhead heights.  If
%   vector IP is neither empty nor negative, initial point markers will
%   be plotted with diameter IP; for default diameter W, use IP = 0.
%   The units of W, H and IP are 1/72 of the PlotBox diagonal.
%
%   ALPHA (default = 1) is a vector of FaceAlpha values ranging between
%   0 (clear) and 1 (opaque).  FaceAlpha is a surface (arrowhead and
%   initial point marker) property and does not affect lines.  FaceAlpha
%   is not supported for 2D rendering.
%
%   BETA (default = 0.4) is a scalar that controls the relative
%   brightness of light and dark color shades, ranging between 0 (no
%   contrast) and 1 (maximum contrast).
%
%   ARROW3(P1,P2,S,W,H,'cone',...) will plot cones with bases centered
%   on P1 in the direction given by P2.  In this instance, P2 is
%   interpreted as a direction vector instead of a terminal point.
%   Neither initial point markers nor lines are plotted with the 'cone'
%   option.
%
%   HN = ARROW3(P1,P2,...) returns a vector of handles to line and
%   surface objects created by ARROW3.
%
%   ARROW3 COLORS will plot a table of named colors with default
%   brightness.  ARROW3('colors',BETA) will plot a table of named
%   colors with brightness BETA.
%
%   ARROW3 attempts to preserve the appearance of existing axes.  In
%   particular, ARROW3 will not change XYZLim, View, or CameraViewAngle.
%   ARROW3 does not, however, support stretch-to-fill scaling.  AXIS
%   NORMAL will restore the current axis box to full size and remove any
%   restrictions on the scaling of units, but will likely result in
%   distorted arrowheads and initial point markers.  See
%   (arrow3_messes_up_my_plots.html).
%
%   If a particular aspect ratio or variable limit is required, use
%   DASPECT, PBASPECT, AXIS, or XYZLIM commands before calling ARROW3.
%   Changing limits or aspect ratios after calling ARROW3 may alter the
%   appearance of arrowheads and initial point markers.  ARROW3 sets
%   XYZCLimMode to manual for all plots, sets DataAspectRatioMode to
%   manual for linear plots, and sets PlotBoxAspectRatioMode to manual
%   for log plots and 3D plots.  CameraViewAngleMode is also set to
%   manual for 3D plots.
%
%   ARROW3 UPDATE will restore the appearance of arrowheads and
%   initial point markers that have become corrupted by changes to
%   limits or aspect ratios.  ARROW3('update',SF) will redraw initial
%   point markers and arrowheads with scale factor SF.  If SF has one
%   element, SF scales W, H and IP.  If SF has two elements, SF(1)
%   scales W and IP, and SF(2) scales H.  If SF has three elements,
%   SF(1) scales W, SF(2) scales H, and SF(3) scales IP.  All sizes are
%   relative to the current PlotBox diagonal.
%
%   ARROW3 UPDATE COLORS will update the magnitude coloring of
%   arrowheads, initial point markers, and lines to conform to the
%   current ColorMap.
%
%   HN = ARROW3('update',...) returns a vector of handles to updated
%   objects.
%
%   EXAMPLES:
%
%     % 2D vectors
%     arrow3([0 0],[1 3])
%     arrow3([0 0],[1 2],'-.e')
%     arrow3([0 0],[10 10],'--x2',1)
%     arrow3(zeros(10,2),50*rand(10,2),'x',1,3)
%     arrow3(zeros(10,2),[10*rand(10,1),500*rand(10,1)],'u')
%     arrow3(10*rand(10,2),50*rand(10,2),'x',1,[],1)
%
%     % 3D vectors
%     arrow3([0 0 0],[1 1 1])
%     arrow3(zeros(20,3),50*rand(20,3),'--x1.5',2)
%     arrow3(zeros(100,3),50*rand(100,3),'x',1,3)
%     arrow3(zeros(10,3),[10*rand(10,1),500*rand(10,1),50*rand(10,1)],'a')
%     arrow3(10*rand(10,3),50*rand(10,3),'x',[],[],0)
%
%     % 3D animation
%     t=(0:pi/40:8*pi)'; u=cos(t); v=sin(t);
%     plot3(20*t,u,v)
%     axis([0,600,-1.5,1.5,-1.5,1.5])
%     grid on, view(35,25)
%     hold on
%     pbaspect([1.8,1.4,1])
%     arrow3(zeros(3),diag([500,1.5,1.5]),'l',0.7,[],0)
%     p=[20*t,u,v]; inc=4:1:length(t);
%     p2=p(inc,:); p1=p(inc-1,:);
%     hn=arrow3(p1(1,:),p2(1,:),'0_b',0.7);
%     for i=2:1:length(p1)
%       delete(hn)
%       hn=arrow3(p1(i,:),p2(i,:),'0_b',0.7);
%       pause(0.01)
%     end
%     hold off
%
%     % Cone plot
%     t=(pi/8:pi/8:2*pi)'; p1=[cos(t) sin(t) t]; p2=repmat([0 0 1],16,1);
%     arrow3(p1,p2,'x',2,4,'cone'), hold on
%     plot3(p1(:,1),p1(:,2),p1(:,3)), hold off
%     pause % change cone size
%     arrow3('update',[1,2])
%
%     % Just for fun
%     arrow3(zeros(100,3),50*rand(100,3),'x',8,4,[],0.95)
%     light('position',[-10 -10 -10],'style','local')
%     light('position',[60,60,60]), lighting gouraud
%
%     % ColorOrder variable, color code prefixes, and Beta
%     global ColorOrder, ColorOrder='^ui^e_hq^v';
%     theta=[0:pi/22:pi/2]';
%     arrow3(zeros(12,2),[cos(theta),sin(theta)],'1.5o',1.5,[],[],[],0.5)
%
%     % ColorOrder property, LineStyleOrder, and LineWidthOrder
%     global ColorOrder, ColorOrder=[];
%     set(gca,'ColorOrder',[1,0,0;0,0,1;0.25,0.75,0.25;0,0,0])
%     set(gca,'LineStyleOrder',{'-','--','-.',':'})
%     global LineWidthOrder, LineWidthOrder=[1,2,4,8];
%     w=[1,2,3,4]; h=[4,6,4,2];
%     arrow3(zeros(4,2),[10*rand(4,1),500*rand(4,1)],'o*/',w,h,0)
%
%     % Magnitude coloring
%     colormap spring
%     arrow3(20*randn(20,3),50*randn(20,3),'|',[],[],0)
%     set(gca,'color',0.7*[1,1,1])
%     set(gcf,'color',0.5*[1,1,1]), grid on, colorbar
%     pause % change the ColorMap and update colors
%     colormap hot
%     arrow3('update','colors')
%
%     % LogLog plot
%     set(gca,'xscale','log','yscale','log');
%     axis([1e2,1e8,1e-2,1e-1]); hold on
%     p1=repmat([1e3,2e-2],15,1);
%     q1=[1e7,1e6,1e5,1e4,1e3,1e7,1e7,1e7,1e7,1e7,1e7,1e6,1e5,1e4,1e3];
%     q2=1e-2*[9,9,9,9,9,7,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%     global ColorOrder, ColorOrder=[];
%     set(gca,'ColorOrder',rand(15,3))
%     arrow3(p1,p2,'o'), grid on, hold off
%
%     % SemiLogX plot
%     set(gca,'xscale','log','yscale','linear');
%     axis([1e2,1e8,1e-2,1e-1]); hold on
%     p1=repmat([1e3,0.05],15,1);
%     q1=[1e7,1e6,1e5,1e4,1e3,1e7,1e7,1e7,1e7,1e7,1e7,1e6,1e5,1e4,1e3];
%     q2=1e-2*[9,9,9,9,9,7,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%     arrow3(p1,p2,'x'), grid on, hold off
%
%     % SemiLogY plot
%     set(gca,'xscale','linear','yscale','log');
%     axis([2,8,1e-2,1e-1]); hold on
%     p1=repmat([3,2e-2],17,1);
%     q1=[7,6,5,4,3,7,7,7,7,7,7,7,7,6,5,4,3];
%     q2=1e-2*[9,9,9,9,9,8,7,6,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%     set(gca,'LineStyleOrder',{'-','--','-.',':'})
%     arrow3(p1,p2,'*',1,[],0), grid on, hold off
%
%     % Color tables
%     arrow3('colors')           % default color table
%     arrow3('colors',0.3)       % low contrast color table
%     arrow3('colors',0.5)       % high contrast color table
%
%     % Update initial point markers and arrowheads
%     % relative to the current PlotBox diagonal
%     arrow3('update')           % redraw same size
%     arrow3('update',2)         % redraw double size
%     arrow3('update',0.5)       % redraw half size
%     arrow3('update',[0.5,2,1]) % redraw W half size,
%                                %        H double size, and
%                                %        IP same size
%
%     See also (arrow3_examples.html), (arrow3_messes_up_my_plots.html).

%   Copyright(c)2002-2020 Version 5.18.1
%     Tom Davis (tdavis@earthsurvey.us)
%     Jeff Chang

%   Revision History:
%
%     08/08/20 - Address ColorOrder, LineStyleOrder property changes
%                in R2019b. (TD)
%     06/19/20 - Remove deprecated zbuffer renderer. (TD)
%     11/02/19 - Address ColorOrder property changes in R2019b. (TD)
%     08/28/18 - Use explicit EdgeColor. (TD)
%     01/15/13 - Use AppData instead of UserData. (TD)
%     07/27/11 - Added animation example. (TD)
%     05/13/09 - Corrected spelling errors (TD)
%     03/16/08 - Updated contact information (TD)
%     10/23/07 - Corrected zero magnitude exclusion (TD)
%     09/08/07 - Added cone plot option; removed adaptive grid
%                spacing; corrected scale factor; removed "nearly"
%                tight limits (TD)
%     07/24/07 - Ignore zero-magnitude input (TD)
%     07/08/07 - Modified named colors to match named Simulink
%                colors; added light and dark shades for basic
%                colors (ymcrgb) (TD)
%     07/01/07 - Modified named colors to match default ColorOrder
%                colors (TD)
%     06/24/07 - Error checking for empty P1, P2 (TD)
%     06/17/07 - Trim colors,W,H,IP,ALPHA to LENGTH(P1) (TD)
%     05/27/07 - Magnitude coloring and documentation revision (TD)
%     03/10/07 - Improved code metrics (TD)
%     02/21/07 - Preserve existing axis appearance;
%                use relative sizes for W, H, and IP;
%                removed version checking; minor bug fixes (TD)
%     01/09/04 - Replaced calls to LINSPACE, INTERP1, and
%                COLORMAP (TD)
%     12/17/03 - Semilog examples, CAXIS support, magnitude
%                coloring, and color updating; use CData instead
%                of FaceColor; minor bug fixes (TD)
%     07/17/03 - Changed 2D rendering from OpenGL to ZBuffer;
%                defined HN for COLORS and UPDATE options (TD)
%     02/27/03 - Replaced calls to RANDPERM, VIEW, REPMAT, SPHERE,
%                and CYLINDER; added ZBuffer for log plots, RESET
%                for CLA and CLF, and ABS for W and H (TD)
%     02/01/03 - Added UPDATE scale factor and MATLAB version
%                checking, replaced call to CROSS (TD)
%     12/26/02 - Added UserData and UPDATE option (TD)
%     11/16/02 - Added more named colors, color code prefix,
%                global ColorOrder, ALPHA , and BETA (TD)
%     10/12/02 - Added global LineWidthOrder,
%                vectorized W, H and IP (TD)
%     10/05/02 - Changed CLF to CLA for subplot support,
%                added ColorOrder and LineStyleOrder support (TD)
%     04/27/02 - Minor log plot revisions (TD)
%     03/26/02 - Added log plot support (TD)
%     03/24/02 - Adaptive grid spacing control to trade off
%                appearance vs. speed based on size of matrix (JC)
%     03/16/02 - Added "axis tight" for improved appearance (JC)
%     03/12/02 - Added initial point marker (TD)
%     03/03/02 - Added aspect ratio support (TD)
%     03/02/02 - Enhanced program's user friendliness (JC)
%                (lump Color, LineStyle, and LineWidth together)
%     03/01/02 - Replaced call to ROTATE (TD)
%     02/28/02 - Modified line plotting,
%                added linewidth and linestyle (TD)
%     02/27/02 - Minor enhancements on 3D appearance (JC)
%     02/26/02 - Minor enhancements for speed (TD&JC)
%     02/26/02 - Optimize PLOT3 and SURF for speed (TD)
%     02/25/02 - Return handler, error handling, color effect,
%                generalize for 2D/3D vectors (JC)
%     02/24/02 - Optimize PLOT3 and SURF for speed (TD)
%     02/23/02 - First release (JC&TD)

%-------------------------------------------------------------------------
% Error Checking
global LineWidthOrder ColorOrder
if nargin<8 || isempty(beta), beta=0.4; end
beta=abs(beta(1)); if nargout, hn=[]; end
if strcmpi(p1,'colors')                            % plot color table
  if nargin>1, beta=abs(p2(1)); end
  LocalColorTable(1,beta); return
end
fig=gcf; ax=gca;
if strcmpi(p1,'update'), ad=getappdata(ax,'Arrow3');    % update
  LocalLogCheck(ax);
  if size(ad,2)<13, error('Invalid AppData'), end
  setappdata(ax,'Arrow3',[]); sf=[1,1,1]; flag=0;
  if nargin>1
    if strcmpi(p2,'colors'), flag=1;               % update colors
    elseif ~isempty(p2)                            % update surfaces
      sf=p2(1)*sf; n=length(p2(:));
      if n>1, sf(2)=p2(2); if n>2, sf(3)=p2(3); end, end
    end
  end
  H=LocalUpdate(fig,ax,ad,sf,flag); if nargout, hn=H; end, return
end
InputError=['Invalid input, type HELP ',upper(mfilename),...
  ' for usage examples'];
if nargin<2, error(InputError), end
[r1,c1]=size(p1); [r2,c2]=size(p2);
if c1<2 || c1>3 || r1*r2==0, error(InputError), end
if r1~=r2, error('P1 and P2 must have same number of rows'), end
if c1~=c2, error('P1 and P2 must have same number of columns'), end
p=sum(abs(p2-p1),2)~=0; cone=0;
if nargin>5 && ~isempty(ip) && strcmpi(ip,'cone')  % cone plot
  cone=1; p=sum(p2,2)~=0;
  if ~any(p), error('P2 cannot equal 0'), end
  set(ax,'tag','Arrow3ConePlot');
elseif ~any(p), error('P1 cannot equal P2')
end
if ~all(p)
  warning('Arrow3:ZeroMagnitude','Zero magnitude ignored')
  p1=p1(p,:); p2=p2(p,:); [r1,c1]=size(p1);
end
n=r1; Zeros=zeros(n,1);
if c1==2, p1=[p1,Zeros]; p2=[p2,Zeros];
elseif ~any([p1(:,3);p2(:,3)]), c1=2; end
L=get(ax,'LineStyleOrder'); C=get(ax,'ColorOrder');
ST=get(ax,'DefaultSurfaceTag'); LT=get(ax,'DefaultLineTag');
EC=get(ax,'DefaultSurfaceEdgeColor');
if strcmp(get(ax,'nextplot'),'add') && strcmp(get(fig,'nextplot'),'add')
  Xr=get(ax,'xlim'); Yr=get(ax,'ylim'); Zr=get(ax,'zlim');
  [xs,ys,xys]=LocalLogCheck(ax); restore=1;
  if xys, mode='auto';
    if any([p1(:,3);p2(:,3)]), error('3D log plot not supported'), end
    if (xs && ~all([p1(:,1);p2(:,1)]>0)) || ...
       (ys && ~all([p1(:,2);p2(:,2)]>0))
       error('Nonpositive log data not supported')
    end
  else, mode='manual';
    if strcmp(get(ax,'WarpToFill'),'on')
      warning('Arrow3:WarpToFill',['Stretch-to-fill scaling not ',...
        'supported;\nuse DASPECT or PBASPECT before calling ARROW3.']);
    end
  end
  set(ax,'XLimMode',mode,'YLimMode',mode,'ZLimMode',mode,...
    'CLimMode','manual');
else, restore=0; cla reset; xys=0; set(fig,'nextplot','add');
  set(ax,'LineStyleOrder',L);
  if c1==2, azel=[0,90]; else, azel=[-37.5,30]; end
  setappdata(ax,'Arrow3',[]);
  set(ax,'nextplot','add','View',azel);
end

%-------------------------------------------------------------------------
% Style Control
[vc,cn]=LocalColorTable(0); prefix='';
if nargin<3, [c,ls,lw]=LocalValidateCLSW;% default color, linestyle/width
else
  [c,ls,lw]=LocalValidateCLSW(s);
  if length(c)>1, if sum('_^'==c(1)), prefix=c(1); end, c=c(2); end
  if c=='x'                              % random named color (less white)
    [~,i]=sort(rand(1,23)); c=cn(i,:);
  elseif c=='o'                                    % ColorOrder
    if ~isempty(ColorOrder)
      [c,failed]=LocalColorMap(lower(ColorOrder),vc,cn,beta);
      if failed, ColorOrderWarning=['Invalid ColorOrder ',...
        'variable, current ColorOrder property will be used'];
        warning('Arrow3:ColorOrder',ColorOrderWarning)
      else, C=c;
      end
    end, c=C;
  elseif c=='|', map=get(fig,'colormap');          % magnitude coloring
    M=(p1-p2); M=sqrt(sum(M.*M,2)); minM=min(M); maxM=max(M);
    if maxM-minM<1, minM=0; end
    set(ax,'clim',[minM,maxM]); c=LocalInterp(minM,maxM,map,M);
  elseif ~sum(vc==c), c='k'; ColorWarning=['Invalid color switch, ',...
    'default color (black) will be used'];
    warning('Arrow3:Color',ColorWarning)
  end
end
if length(c)==1                                    % single color
  c=LocalColorMap([prefix,c],vc,cn,beta);
end
set(ax,'ColorOrder',c); c=LocalRepmat(c,[ceil(n/size(c,1)),1]);
if ls~='*'                                         % LineStyleOrder
  set(ax,'LineStyleOrder',ls);
end
ls=get(ax,'LineStyleOrder');
ls=LocalRepmat(cellstr(ls),[ceil(n/size(ls,1)),1]);
if lw=='/'                                         % LineWidthOrder
  if ~isempty(LineWidthOrder)
    lw=LineWidthOrder;
  else, lw=0.5; LineWidthOrderWarning=['Undefined LineWidthOrder, ',...
    'default width (0.5) will be used'];
    warning('Arrow3:LineWidthOrder',LineWidthOrderWarning)
  end
end
lw=LocalRepmat(lw,[ceil(n/length(lw)),1]);
if nargin<4 || isempty(w), w=1; end                % width
w=LocalRepmat(abs(w(:)),[ceil(n/length(w)),1]);
if nargin<5 || isempty(h), h=3*w; end              % height
h=LocalRepmat(abs(h(:)),[ceil(n/length(h)),1]);
if nargin>5 && ~isempty(ip) && ~cone               % ip
  ip=LocalRepmat(ip(:),[ceil(n/length(ip)),1]);
  i=find(ip==0); ip(i)=w(i);
else, ip=-ones(n,1);
end
if nargin<7 || isempty(alpha), alpha=1; end
a=LocalRepmat(alpha(:),[ceil(n/length(alpha)),1]); % FaceAlpha

%-------------------------------------------------------------------------
% Log Plot
if xys
  units=get(ax,'units'); set(ax,'units','points');
  pos=get(ax,'position'); set(ax,'units',units);
  if strcmp(get(ax,'PlotBoxAspectRatioMode'),'auto')
    set(ax,'PlotBoxAspectRatio',[pos(3),pos(4),1]);
  end
  par=get(ax,'PlotBoxAspectRatio');
  set(ax,'DataAspectRatio',[par(2),par(1),par(3)]);
  % map coordinates onto unit square
  q=[p1;p2]; xr=Xr; yr=Yr;
  if xs, xr=log10(xr); q(:,1)=log10(q(:,1)); end
  if ys, yr=log10(yr); q(:,2)=log10(q(:,2)); end
  q=q-LocalRepmat([xr(1),yr(1),0],[2*n,1]);
  dx=xr(2)-xr(1); dy=yr(2)-yr(1);
  q=q*diag([1/dx,1/dy,1]);
  q1=q(1:n,:); q2=q(n+1:end,:);
else, xs=0; ys=0; dx=0; dy=0; xr=0; yr=0;
end

%-------------------------------------------------------------------------
% Line
set(ax,'DefaultLineTag','arrow3');
if ~cone
  H1=Zeros;
  for i=1:n
    if lw(i)
      H1(i)=plot3([p1(i,1),p2(i,1)],[p1(i,2),p2(i,2)],...
        [p1(i,3),p2(i,3)],ls{i},'Color',c(i,:),'LineWidth',lw(i));
    end
  end
else                                               % cone plot
  P=zeros(3*n,3); i=1:n;
  P(3*i-2,:)=p1(i,:); P(3*i-1,:)=p1(i,:); P(3*i,1)=NaN;
  H1=plot3(P(:,1),P(:,2),P(:,3));
end

%-------------------------------------------------------------------------
% Scale
if ~restore, axis tight, end
ar=get(ax,'DataAspectRatio'); ar=sqrt(3)*ar/norm(ar);
set(ax,'DataAspectRatioMode','manual');
if xys, sf=1;
else, xr=get(ax,'xlim'); yr=get(ax,'ylim'); zr=get(ax,'zlim');
  sf=norm(diff([xr;yr;zr],1,2)./ar')/72;
end

%-------------------------------------------------------------------------
% AppData
c=c(1:n,:); w=w(1:n); h=h(1:n); ip=ip(1:n); a=a(1:n);
setappdata(ax,'Arrow3',[getappdata(ax,'Arrow3');p1,p2,c,w,h,ip,a]);

%-------------------------------------------------------------------------
% Arrowhead
whip=sf*[w,h,ip];
if xys, whip=whip*sqrt(2)/72; p1=q1; p2=q2; end
w=whip(:,1); h=whip(:,2); ip=whip(:,3);
if cone                                            % cone plot
  delete(H1), H1=[];
  p2=p2./LocalRepmat(sqrt(sum(p2.*p2,2)),[1,3]);
  p2=p1+p2.*LocalRepmat(ar,[n,1]).*LocalRepmat(h,[1,3]);
end
W=(p1-p2)./LocalRepmat(ar,[n,1]);
W=W./LocalRepmat(sqrt(sum(W.*W,2)),[1,3]);         % new z direction
U=[-W(:,2),W(:,1),Zeros];
N=sqrt(sum(U.*U,2)); i=find(N<eps); j=length(i);
U(i,:)=LocalRepmat([1,0,0],[j,1]); N(i)=ones(j,1);
U=U./LocalRepmat(N,[1,3]);                         % new x direction
V=[W(:,2).*U(:,3)-W(:,3).*U(:,2),...               % new y direction
  W(:,3).*U(:,1)-W(:,1).*U(:,3),...
  W(:,1).*U(:,2)-W(:,2).*U(:,1)];
m=20;                               % surface grid spacing
set(ax,'DefaultSurfaceTag','arrow3','DefaultSurfaceEdgeColor','none');
r=[0;1]; theta=(0:m)/m*2*pi; Ones=ones(1,m+1);
x=r*cos(theta); y=r*sin(theta); z=r*Ones;
G=surface(x/2,y/2,z); dar=diag(ar);
set(G,'EdgeColor','none');
X=get(G,'XData'); Y=get(G,'YData'); Z=get(G,'ZData');
H2=Zeros; [j,k]=size(X);
for i=1:n                           % translate, rotate, and scale
  H2(i)=copyobj(G,ax);
  xyz=[w(i)*X(:),w(i)*Y(:),h(i)*Z(:)]*[U(i,:);V(i,:);W(i,:)]*dar;
  x=reshape(xyz(:,1),j,k)+p2(i,1);
  y=reshape(xyz(:,2),j,k)+p2(i,2);
  z=reshape(xyz(:,3),j,k)+p2(i,3);
  LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,...
    x,y,z,a(i),c(i,:),H2(i),2,m+1);
end
delete(G);

%-------------------------------------------------------------------------
% Initial Point Marker
if any(ip>0)
  theta=(-m:2:m)/m*pi; phi=(-m:2:m)'/m*pi/2; cosphi=cos(phi);
  x=cosphi*cos(theta); y=cosphi*sin(theta); z=sin(phi)*Ones;
  G=surface(x*ar(1)/2,y*ar(2)/2,z*ar(3)/2);
  set(G,'EdgeColor','none');
  X=get(G,'XData'); Y=get(G,'YData'); Z=get(G,'ZData');
  H3=zeros(n,1);
  for i=1:n                                        % translate
    if ip(i)>0
      H3(i)=copyobj(G,ax);
      x=p1(i,1)+X*ip(i); y=p1(i,2)+Y*ip(i); z=p1(i,3)+Z*ip(i);
      LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,...
        x,y,z,a(i),c(i,:),H3(i),m+1,m+1);
    end
  end, delete(G);
else, H3=[];
end

%-------------------------------------------------------------------------
% Finish
if restore, xr=Xr; yr=Yr; zr=Zr;
  if xys, set(ax,'DataAspectRatioMode','auto'); end
else
  axis tight
  xr=get(ax,'xlim'); yr=get(ax,'ylim'); zr=get(ax,'zlim');
  set(ax,'nextplot','replace');
end
azel=get(ax,'view');
if abs(azel(2))==90, c1=2; else, c1=3; end
set(fig,'Renderer','OpenGL');
set(ax,'LineStyleOrder',L,'ColorOrder',C,'DefaultLineTag',LT,...
  'DefaultSurfaceTag',ST,'DefaultSurfaceEdgeColor',EC,...
  'xlim',xr,'ylim',yr,'zlim',zr,'clim',get(ax,'CLim'));
if c1==3, set(ax,'CameraViewAngle',get(ax,'CameraViewAngle'),...
   'PlotBoxAspectRatio',get(ax,'PlotBoxAspectRatio'));
end
if nargout, hn=[H1(:);H2(:);H3(:)]; end

%-------------------------------------------------------------------------
% Local Functions
%-------------------------------------------------------------------------
% Update
function H=LocalUpdate(fig,ax,ad,sf,flag)
global ColorOrder
p1=ad(:,1:3); p2=ad(:,4:6); c=ad(:,7:9); a=ad(:,13);
w=sf(1)*ad(:,10); h=sf(2)*ad(:,11); ip=sf(3)*ad(:,12);
H=get(ax,'children'); tag=get(H,'tag'); type=get(H,'type');
delete(H(strcmp(tag,'arrow3') & strcmp(type,'surface')));
set(fig,'nextplot','add'); set(ax,'nextplot','add'); H1=[];
if flag, map=get(fig,'colormap');                  % update colors
  M=(p1-p2); M=sqrt(sum(M.*M,2)); minM=min(M); maxM=max(M);
  H1=H(strcmp(tag,'arrow3') & strcmp(type,'line'));
  MagnitudeWarning=['Cannot perform magnitude coloring on lines ',...
    'that\nwere drawn with a single color, linestyle, and linewidth'];
  if length(H1)>1
    for i=1:length(H1)                             % update line colors
      x=get(H1(i),'xdata'); y=get(H1(i),'ydata'); z=get(H1(i),'zdata');
      if length(x)>2                               % multiple lines
        warning('Arrow3:Magnitude',MagnitudeWarning), continue
      end
      m=sqrt((x(1)-x(2))^2+(y(1)-y(2))^2+(z(1)-z(2))^2);
      c=LocalInterp(minM,maxM,map,m); set(H1(i),'color',c);
    end
  elseif length(H1)==1
    warning('Arrow3:Magnitude',MagnitudeWarning)
  end
  c=LocalInterp(minM,maxM,map,M);
end
set(ax,'ColorOrder',c);                            % update surfaces
ColorOrder=[];
if strcmp(get(ax,'tag'),'Arrow3ConePlot')
      H=arrow3(p1,p2,'o' ,w,h,'cone',a);           % update cones
else, H=arrow3(p1,p2,'o0',w,h,    ip,a);
end,  H=[H1(:);H(:)];
set(ax,'nextplot','replace');

%-------------------------------------------------------------------------
% SetSurface
function LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,x,y,z,a,c,H,n,m)
if xys
  x=x*dx+xr(1); y=y*dy+yr(1);
  if xs, x=10.^x; end
  if ys, y=10.^y; end
end
cd=zeros(n,m,3); cd(:,:,1)=c(1); cd(:,:,2)=c(2); cd(:,:,3)=c(3);
set(H,'XData',x,'YData',y,'ZData',z,'CData',cd,'FaceAlpha',a);

%-------------------------------------------------------------------------
% ColorTable
function [vc,cn]=LocalColorTable(n,beta)
vc='kymcrgbadefhijlnpqstuvzw';                     % valid color codes
%                k               y               m               c
cn=[0.00,0.00,0.00; 1.00,1.00,0.00; 1.00,0.00,1.00; 0.00,1.00,1.00;
%                r               g               b               a
    1.00,0.00,0.00; 0.00,1.00,0.00; 0.00,0.00,1.00; 0.42,0.59,0.24;
%                d               e               f               h
    0.25,0.25,0.25; 0.00,0.50,0.00; 0.70,0.13,0.13; 1.00,0.41,0.71;
%                i               j               l               n
    0.29,0.00,0.51; 0.00,0.66,0.42; 0.50,0.50,0.50; 0.50,0.20,0.00;
%                p               q               s               t
    0.75,0.75,0.00; 1.00,0.50,0.00; 0.00,0.75,0.75; 0.80,0.34,0.00;
%                u               v               z               w
    0.50,0.00,0.13; 0.75,0.00,0.75; 0.38,0.74,0.99; 1.00,1.00,1.00];

% Named Simulink Colors (zaql)
% LightBlue = 0.38  0.74  0.99 = aZure
% DarkGreen = 0.42  0.59  0.24 = Asparagus
% Orange    = 1.00  0.50  0.00 = kumQuat
% Gray      = 0.50  0.50  0.50 = Light gray
%
% Default ColorOrder Property Colors (bersvpd)
% Color1    = 0.00  0.00  1.00 = Blue
% Color2    = 0.00  0.50  0.00 = Evergreen
% Color3    = 1.00  0.00  0.00 = Red
% Color4    = 0.00  0.75  0.75 = Sky blue
% Color5    = 0.75  0.00  0.75 = Violet
% Color6    = 0.75  0.75  0.00 = Pear
% Color7    = 0.25  0.25  0.25 = Dark gray

if n, clf reset                                    % plot color table
  name={'blacK','Yellow','Magenta','Cyan',...
    'Red','Green','Blue','Asparagus',...
    'Dark gray','Evergreen','Firebrick','Hot pink',...
    'Indigo','Jade','Light gray','Nutbrown',...
    'Pear','kumQuat','Sky blue','Tawny',...
    'bUrgundy','Violet','aZure','White'};
  c=['yptn';'gjae';'czsb';'hmvi';'qrfu';'wldk'];
  set(gcf,'DefaultAxesXTick',[],'DefaultAxesYTick',[],...
    'DefaultAxesXTickLabel',[],'DefaultAxesYTickLabel',[],...
    'DefaultAxesXLim',[0,0.75],'DefaultAxesYLim',[0,0.75],...
    'DefaultRectangleEdgeColor','none');
  for i=1:24, subplot(4,6,i); box on
    j=find(vc==c(i)); title(name{j});
    dark=LocalBrighten(cn(j,:),-beta);
    light=LocalBrighten(cn(j,:),beta);
    rectangle('Position',[0,0.00,0.75,0.25],'FaceColor',dark);
    rectangle('Position',[0,0.25,0.75,0.25],'FaceColor',cn(j,:));
    rectangle('Position',[0,0.50,0.75,0.25],'FaceColor',light);
    rectangle('Position',[0,0.00,0.75,0.75],'EdgeColor','k');
    if rem(i,6)==1
      set(gca,'YTickLabel',{'dark','normal','light'},...
        'YTick',[0.125,0.375,0.625]);
      if i==19
        text(0,-0.25,['{\bf\itARROW3}  Named Color Table  ',...
            '( \beta = ',num2str(beta),' )']);
      end
    end
  end
end

%-------------------------------------------------------------------------
% ColorMap
function [C,failed]=LocalColorMap(c,vc,cn,beta)
n=length(c); failed=0; C=zeros(n,3); i=1; j=1;
while 1
  if ~sum([vc,'_^']==c(i)), failed=1; break, end
  if sum('_^'==c(i))
    if i+1>n, failed=1; break, end
    if ~sum(vc==c(i+1)), failed=1; break, end
    cc=cn(vc==c(i+1),:); gamma=beta;
    if c(i)=='_', gamma=-beta; end
    C(j,:)=LocalBrighten(cc,gamma); i=i+2;
  else, C(j,:)=cn(vc==c(i),:); i=i+1;
  end
  if i>n, break, end, j=j+1;
end
if n>j, C(j+1:n,:)=[]; end

%-------------------------------------------------------------------------
% Brighten
function C=LocalBrighten(c,beta)
if sum([c==0,c==1])==3 && sum(c==0)<3 && sum(c==1)<3
  if beta<0
    C=(1+beta)*c;
  else
    C=c;  C(C==0)=beta;
  end
else
  C=c.^((1-min(1-sqrt(eps),abs(beta)))^sign(beta));
end

%-------------------------------------------------------------------------
% Repmat
function B=LocalRepmat(A,siz)
if length(A)==1, B(prod(siz))=A; B(:)=A; B=reshape(B,siz);
else, [m,n]=size(A); mind=(1:m)'; nind=(1:n)';
  mind=mind(:,ones(1,siz(1))); nind=nind(:,ones(1,siz(2)));
  B=A(mind,nind);
end

%-------------------------------------------------------------------------
% Interp
function v=LocalInterp(xmin,xmax,y,u)
[m,n]=size(y); h=(xmax-xmin)/(m-1); p=length(u); v=zeros(p,n);
k=min(max(1+floor((u-xmin)/h),1),m-1); s=(u-xmin)/h-k+1;
for j=1:n, v(:,j)=y(k,j)+s.*(y(k+1,j)-y(k,j)); end
v(v<0)=0; v(v>1)=1;

%-------------------------------------------------------------------------
% Check for supported log scales
function [xs,ys,xys]=LocalLogCheck(ax)
xs=strcmp(get(ax,'xscale'),'log');
ys=strcmp(get(ax,'yscale'),'log');
zs=strcmp(get(ax,'zscale'),'log');
if zs, error('Z log scale not supported'), end
xys=xs+ys;
if xys, azel=get(ax,'view');
  if abs(azel(2))~=90, error('3D log plot not supported'), end
end

%-------------------------------------------------------------------------
% Generate valid value for color, linestyle and linewidth
function [c,ls,lw]=LocalValidateCLSW(s)
if nargin<1 || isempty(s), c='k'; ls='-'; lw=0.5;
else
  % identify linestyle
  if contains(s,'--'), ls='--'; s=strrep(s,'--','');
  elseif contains(s,'-.'), ls='-.'; s=strrep(s,'-.','');
  elseif contains(s,'-'), ls='-'; s=strrep(s,'-','');
  elseif contains(s,':'), ls=':'; s=strrep(s,':','');
  elseif contains(s,'*'), ls='*'; s=strrep(s,'*','');
  else, ls='-';
  end

  % identify linewidth
  tmp=double(s);
  tmp=find(tmp>45 & tmp<58);
  if ~isempty(tmp)
    if any(s(tmp)=='/'), lw='/'; else, lw=str2double(s(tmp)); end
    s(tmp)='';
  else, lw=0.5;
  end

  % identify color
  if ~isempty(s), s=lower(s);
    if length(s)>1, c=s(1:2);
    else, c=s(1); end
  else, c='k';
  end
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function y = contract(x,i,j)
%CONTRACT Contract tensor along two dimensions (array trace).
%
%   Y = CONTRACT(X,I,J) contracts the entries of X along dimensions I
%   and J. Contraction is a generalization of matrix trace. In other
%   words, the trace is performed along the two-dimensional slices
%   defined by dimensions I and J. It is possible to implement tensor
%   multiplication as an outer product followed by a contraction.
%
%   Examples
%   X = tensor(rand(4,3,2)); Y = tensor(rand(3,2,4));
%   Z1 = ttt(X,Y,1,3); %<-- Normal tensor multiplication
%   Z2 = contract(ttt(X,Y),1,6); %<-- Outer product + contract
%   norm(Z1-Z2) %<-- Should be zero
%
%   See also TENSOR, TENSOR/TTT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% Error checking
if x.size(i) ~= x.size(j)
    error('Must contract along equally sized dimensions');
end

% Error checking
if i == j
    error('Must contract along two different dimensions');
end

% Easy case - returns a scalar
if ndims(x) == 2
    y = trace(x.data);
    return;
end

% Remaining dimensions after trace
remdims = setdiff(1:ndims(x),[i j]);

% Size for y
newsize = x.size(remdims);

% Total size of remainder
m =  prod(newsize);

% Number of items to add for trace
n = x.size(i);

% Permute trace dimensions to the end
x = permute(x, [remdims i j]);

% Reshape data to be 3D
data = reshape(x.data, m, n, n);

% Add diagonal entries for each slice
newdata = zeros(m,1);
for i = 1:n
    newdata = newdata + data(:,i,i);
end

% Reshape result
if numel(newsize) > 1
    newdata = reshape(newdata,newsize);
end
y = tensor(newdata,newsize);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%
%   Copyright (C) 2016  Starsky Wong <sununs11@gmail.com>
%
%   Note: The SIFT algorithm is patented in the United States and cannot be
%   used in commercial products without a license from the University of
%   British Columbia.  For more information, refer to the file LICENSE
%   that accompanied this distribution.

function [feat] = hist2Descr(feat,descr,descr_mag_thr)
% Function: Convert histogram to descriptor
descr = descr/norm(descr);
descr = min(descr_mag_thr,descr);
descr = descr/norm(descr);
feat.descr = descr;
end

repo: sw-sift


text:
function dist = plane_exp_point_dist_3d ( p1, p2, p3, p )

%*****************************************************************************80
%
%% PLANE_EXP_POINT_DIST_3D: distance ( explicit plane, point ) in 3D.
%
%  Discussion:
%
%    The explicit form of a plane in 3D is
%
%      the plane through P1, P2 and P3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), P3(3), three points on the plane.
%
%    Input, real P(3), the coordinates of the point.
%
%    Output, real DIST, the distance from the point to the plane.
%
  [ a, b, c, d ] = plane_exp2imp_3d ( p1, p2, p3 );

  dist = plane_imp_point_dist_3d ( a, b, c, d, p );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function geometry_test058 ( )

%*****************************************************************************80
%
%% TEST058 tests PLANE_IMP_POINT_DIST_3D, PLANE_IMP_POINT_DIST_SIGNED_3D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 4;
%
%  This is the plane Z = 10.
%
  a =    0.0;
  b =    0.0;
  c =    1.0;
  d = - 10.0;

  ptest(1:dim_num,1:ntest) = [ ...
    -12.0, 14.0,  0.0; ...
      7.0,  8.0,  9.0; ...
      1.0,  2.0, 10.0; ...
      0.0,  0.0, 12.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST058\n' );
  fprintf ( 1, '  PLANE_IMP_POINT_DIST_3D computes the distance\n' );
  fprintf ( 1, '    between an implicit plane and a point in 3D;\n' );
  fprintf ( 1, '  PLANE_IMP_POINT_DIST_SIGNED 3D computes the\n' );
  fprintf ( 1, '    signed distance between an implicit plane\n' );
  fprintf ( 1, '    and a point in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For all tests, we use the implicit plane with\n' );
  fprintf ( 1, '  (A,B,C,D) = %f  %f  %f  %f\n', a, b, c, d );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  (X,Y,Z)  DISTANCE   SIGNED_DISTANCE\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    p(1:dim_num) = ptest(1:dim_num,i);

    dist = plane_imp_point_dist_3d ( a, b, c, d, p );

    dist_signed = plane_imp_point_dist_signed_3d ( a, b, c, d, p );

    fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', p(1:dim_num), dist, dist_signed );

  end

  return
end


text:
%
%  Generate an equal ripple minimum phase filter starting with a linear
%  phase filter.
%
%  Called by dzmp.

%  written by John Pauly, 1992
%  (c) Board of Trustees, Leland Stanford Junior University

function hmp = fmp(h)

l = length(h);
if rem(l,2) == 0,
   disp('filter length must be odd');
   return;
end;
lp = 8*exp(ceil(log(l)/log(2))*log(2));
hp = [zeros(1,ceil((lp-l)/2)) h zeros(1,floor((lp-l)/2))];
hpf = fftc(hp);
hpfs = hpf-min(real(hpf))*1.000001;
hpfmp = mag2mp(sqrt(abs(hpfs)));
hpmp = ifft(fftshift(conj(hpfmp)));
hmp = hpmp(1:(l+1)/2);


text:
function Fs = fu2F7(u)

Z = lin_fm(u);

NullSp   = null(Z);
if size(NullSp,2) > 2
   Fs = [];
   return; %degenerated sample
end

F1    = reshape(NullSp(:,1),3,3);
F2    = reshape(NullSp(:,2),3,3);
p = fslcm(F1,F2);
aroots = rroots(p);

%xr = o_fslcm(F1,F2)

%aroots == xr

for i = 1:length(aroots)
   l  = aroots(i);
   Ft    = F1 * l + F2 * (1-l);
%   Ft = Ft /norm(Ft,2);
   Fs(:,:,i) = Ft;
end

author: strawlab

repo: MultiCamSelfCal