text:
%% DEMO_febio_0037_lattice_test_octet_truss_01
% Below is a demonstration for:
%
% * Building the geometry for the octet truss lattive with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement and stress results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * compression, tension, compressive, tensile
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * Lattice
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=1;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
markerSize=25;
markerSize2=25;
cMap=gjet(4);

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress
febioLogFileName_stiffness=[febioFebFileNamePart,'_stiffness_out.txt']; %Log file name for exporting stiffness

%Specifying dimensions and number of elements
sampleSize=10;
latticeType=1;
elementType='hex8'; %'hex8'
strutThickness=0.5;

%Define applied displacement
appliedStrain=0.3; %Linear strain (Only used to compute applied stretch)
loadingOption='compression'; % or 'tension'
switch loadingOption
    case 'compression'
        stretchLoad=1-appliedStrain; %The applied stretch for uniaxial loading
    case 'tension'
        stretchLoad=1+appliedStrain; %The applied stretch for uniaxial loading
end
displacementMagnitude=(stretchLoad*sampleSize)-sampleSize; %The displacement magnitude

%Material parameter set
E_youngs1=0.1; %Material Young's modulus
nu1=0.4; %Material Poisson's ratio

% FEA control settings
numTimeSteps=20; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=(1/numTimeSteps); %Maximum time step size
min_residual=1e-20;
symmetric_stiffness=0;
runMode='external'; %'internal' or 'external'

%%

%Specifying dimensions and number of elements
r=0.5; %Radii, results in a width of 1
n=3;
nCopies=n*ones(1,3); %Number of offset copies
d=2*r; %Diameter
w=(n-1)*d; %sampleSize

shrinkFactor=strutThickness./((sampleSize./n).*(sqrt(2)./2));

%% Create lattice

switch latticeType
    case 1 %Octet truss
        [E,V,C,F,CF]=rhombicDodecahedronMesh(r,nCopies);
        V=V./(n-1);
        V=V*sampleSize;

        [indBoundary]=tesBoundary(F);
        cPar.shrinkFactor=shrinkFactor; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=1;
        cPar.latticeSide=2; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [Es,Vs,Cs]=element2lattice(E,V,cPar); %Get lattice structure

        logicKeep1=~(Vs(:,1)<=-1e-3);
        logicKeep2=~(Vs(:,2)<=-1e-3);
        logicKeep3=~(Vs(:,3)<=-1e-3);
        logicKeep4=~(Vs(:,1)>=sampleSize+1e-3);
        logicKeep5=~(Vs(:,2)>=sampleSize+1e-3);
        logicKeep6=~(Vs(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(Es),2)>=4 &...
            sum(logicKeep2(Es),2)>=4 &...
            sum(logicKeep3(Es),2)>=4 &...
            sum(logicKeep4(Es),2)>=4 &...
            sum(logicKeep5(Es),2)>=4 &...
            sum(logicKeep6(Es),2)>=4;

        Es=Es(logicKeepEs,:);
        Cs=Cs(logicKeepEs,:);
        [Es,Vs,indFix]=patchCleanUnused(Es,Vs);

        % [Es,Vs,~,~]=subHex(Es,Vs,1,1);
        % Cs=repmat(Cs,8,1);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs);
        Fb=Fs(indB,:);
    case 2 %Rhombic dodecahedron mesh ("dual" of octet truss lattice)
        [E,V,C,F,CF]=rhombicDodecahedronMesh(r,nCopies);
        V=V./(n-1);
        V=V*sampleSize;

        [indBoundary]=tesBoundary(F);
        cPar.shrinkFactor=0.15; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=3;
        cPar.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [Es,Vs,Cs]=element2lattice(E,V,cPar); %Get lattice structure

        logicKeep1=~(Vs(:,1)<=-1e-3);
        logicKeep2=~(Vs(:,2)<=-1e-3);
        logicKeep3=~(Vs(:,3)<=-1e-3);
        logicKeep4=~(Vs(:,1)>=sampleSize+1e-3);
        logicKeep5=~(Vs(:,2)>=sampleSize+1e-3);
        logicKeep6=~(Vs(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(Es),2)>=4 &...
            sum(logicKeep2(Es),2)>=4 &...
            sum(logicKeep3(Es),2)>=4 &...
            sum(logicKeep4(Es),2)>=4 &...
            sum(logicKeep5(Es),2)>=4 &...
            sum(logicKeep6(Es),2)>=4;

        Es=Es(logicKeepEs,:);
        Cs=Cs(logicKeepEs,:);
        [Es,Vs,indFix]=patchCleanUnused(Es,Vs);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs);
        Fb=Fs(indB,:);
end

if strcmp(elementType,'hex20')
    [Es,Vs,~,Fb]=hex8_hex20(Es,Vs,{},Fb);
end
%%
% Visualizing input mesh and lattic structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

% Fst=[Fs(:,[1 2 3]); Fs(:,[3 4 1]);];
% indB=tesBoundary(Fst,Vs);
% Fbt=Fst(indB,:);

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,Vs,'bw','k',1);
% plotV(Vs(Fb(:),:),'r.');
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%% DEFINE BC's

% Define node set logics
indAll=(1:1:size(Vs,1))';
logicBoundary=ismember(indAll,Fb);

Z=Vs(:,3);
logicTop=Z>=(sampleSize-eps(sampleSize))& logicBoundary;
logicBottom=Z<=eps(sampleSize) & logicBoundary;

X=Vs(:,1);
logicSide1=X>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide2=X<=eps(sampleSize)& logicBoundary;

Y=Vs(:,2);
logicSide3=Y>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide4=Y<=eps(sampleSize)& logicBoundary;

bcPrescribeListCell{1}=find(logicSide1)';
bcPrescribeListCell{2}=find(logicSide2)';
bcPrescribeListCell{3}=find(logicSide3)';
bcPrescribeListCell{4}=find(logicSide4)';
bcPrescribeListCell{5}=find(logicTop)';
bcPrescribeListCell{6}=find(logicBottom)';

%% Smoothing lattice

% indKeep=unique([bcPrescribeListCell{:}]);
% [Fb_clean,Vb_clean,indFix]=patchCleanUnused(Fb,Vs);
%
% cPar.Method='HC';
% cPar.n=6;
%
% cPar.RigidConstraints=indFix(indKeep);
% % cPar.RigidConstraints=cPar.RigidConstraints(cPar.RigidConstraints>0);
%
% [Vb_clean]=tesSmooth(Fb_clean,Vb_clean,[],cPar);
% ind=Fb(:);
% ind=unique(ind(:));
% Vs(ind,:)=Vb_clean;

% cFigure; hold on;
% gpatch(Fb,Vs,'bw','k',1);
% % patchNormPlot(Fs,Vs);
% % plotV(Vs(indKeep,:),'k.','MarkerSize',25)
% axisGeom(gca,fontSize);
% camlight headlight; lighting flat;
% drawnow;

%%

%Prescribed displacement nodes
bcPrescribeList=find(logicTop);
bcSupportList=find(logicBottom);

%%
% Visualizing input mesh and lattice structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,Vs,'bw');
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%%
% Visualize BC's

hf=cFigure; hold on;
title('Boundary conditions model','FontSize',fontSize);
gpatch(Fb,Vs,'w','none',0.5);
hl2(1)=plotV(Vs(bcPrescribeList,:),'r.','MarkerSize',markerSize2);
hl2(2)=plotV(Vs(bcSupportList,:),'b.','MarkerSize',markerSize2);
legend(hl2,{'BC prescribe','BC support'});
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Check porosity
% Note it may be better to use the convex hull volume rather than the cube
% for the octet truss
vol_lattice=sum(hexVol(Es,Vs)); %Volume of hexahedra
porosity_lattice=vol_lattice./sampleSize.^3; %Porosity

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='neo-Hookean';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.E=E_youngs1;
febio_spec.Material.material{1}.v=nu1;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(Vs,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=Vs; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type=elementType; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(Es,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=Es; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

nodeSetName2='bcPrescribeList';
febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{2}.dofs='x,y';

febio_spec.Boundary.bc{3}.ATTR.type='prescribe';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{3}.dof='z';
febio_spec.Boundary.bc{3}.scale.ATTR.lc=1;
febio_spec.Boundary.bc{3}.scale.VAL=displacementMagnitude;
febio_spec.Boundary.bc{3}.relative=0;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

% febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
% febio_spec.Output.logfile.element_data{1}.ATTR.data='sz';
% febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
% febio_spec.Output.logfile.element_data{1}.VAL=1:size(Es,1);

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode=runMode; %Run in external or in matlab terminal

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    % Importing nodal displacements from a log file
    [timeVec, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements
    timeVec=[0; timeVec(:)]; %Time

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN(:,3).^2,2));
    Vs_def=Vs+DN;

    %     % Importing element stress from a log file
    %     [time_mat, E_stress_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress)); %Nodal forces
    %     time_mat=[0; time_mat(:)]; %Time
    %     stress_cauchy_sim=[0; mean(squeeze(E_stress_mat(:,end,:)),1)'];
    %%
    % Importing nodal forces from a log file

    [dataStruct]=importFEBio_logfile(fullfile(savePath,febioLogFileName_force),1,1); %Nodal forces

    %Access data
    timeVec=dataStruct.time;
    f_sum_x=squeeze(sum(dataStruct.data(bcPrescribeList,1,:),1));
    f_sum_y=squeeze(sum(dataStruct.data(bcPrescribeList,2,:),1));
    f_sum_z=squeeze(sum(dataStruct.data(bcPrescribeList,3,:),1));

    %%
    % Visualize force data

    displacementApplied=timeVec.*displacementMagnitude;

    cFigure; hold on;
    xlabel('$u$ [mm]','Interpreter','Latex');
    ylabel('$F_z$ [N]','Interpreter','Latex');
    hp=plot(displacementApplied(:),f_sum_z(:),'b-','LineWidth',3);
    grid on; box on; axis square; axis tight;
    set(gca,'FontSize',fontSize);
    drawnow;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(Fb,Vs_def,DN_magnitude,'k',1); %Add graphics object to animate
    %     gpatch(Fb,Vs,'kw','none',0.25); %A static graphics object
    hp.FaceColor='interp';

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    clim([0 max(DN_magnitude)]);
    axis([min([Vs_def(:,1);Vs(:,1)]) max([Vs_def(:,1);Vs(:,1)])...
        min([Vs_def(:,2);Vs(:,2)]) max([Vs_def(:,2);Vs(:,2)])...
        min([Vs_def(:,3);Vs(:,3)]) max([Vs_def(:,3);Vs(:,3)]) ]); %Set axis limits statically
    %     view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude
        Vs_def=Vs+DN; %Current nodal coordinates

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={Vs_def,DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

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 value = zeta ( p )

%*****************************************************************************80
%
%% ZETA estimates the Riemann Zeta function.
%
%  Discussion:
%
%    For 1 < P, the Riemann Zeta function is defined as:
%
%      ZETA ( P ) = Sum ( 1 <= N < Infinity ) 1 / N**P
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, real P, the power to which the integers are raised.
%    P must be greater than 1.  For integral P up to 20, a
%    precomputed value of ZETA is returned; otherwise the infinite
%    sum is approximated.
%
%    Output, real VALUE, an approximation to the Riemann
%    Zeta function.
%
  if ( p <= 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ZETA - Fatal error!\n' );
    fprintf ( 1, '  Exponent P <= 1.0.\n' );
    error ( 'ZETA - Fatal error!' );
  elseif ( p == 2.0 )
    value = pi^2 / 6.0;
  elseif ( p == 3.0 )
    value = 1.2020569032;
  elseif ( p == 4.0 )
    value = pi^4 / 90.0;
  elseif ( p == 5.0 )
    value = 1.0369277551;
  elseif ( p == 6.0 )
    value = pi^6 / 945.0;
  elseif ( p == 7.0 )
    value = 1.0083492774;
  elseif ( p == 8.0 )
    value = pi^8 / 9450.0;
  elseif ( p == 9.0 )
    value = 1.0020083928;
  elseif ( p == 10.0 )
    value = pi^10 / 93555.0;
  elseif ( p == 11.0 )
    value = 1.0004941886;
  elseif ( p == 12.0 )
    value = 1.0002460866;
  elseif ( p == 13.0 )
    value = 1.0001227133;
  elseif ( p == 14.0 )
    value = 1.0000612482;
  elseif ( p == 15.0 )
    value = 1.0000305882;
  elseif ( p == 16.0 )
    value = 1.0000152823;
  elseif ( p == 17.0 )
    value = 1.0000076372;
  elseif ( p == 18.0 )
    value = 1.0000038173;
  elseif ( p == 19.0 )
    value = 1.0000019082;
  elseif ( p == 20.0 )
    value = 1.0000009540;
  else

    zsum = 0.0;
    n = 0;

    while ( 1 )

      n = n + 1;
      zsum_old = zsum;
      zsum = zsum + 1.0 / n^p;

      if ( zsum <= zsum_old )
        break
      end

    end

    value = zsum;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% TOOL
%
% Files
%   findedge         - highlights edges in certain range
%   findedge3        - highlights edges in certain range in 3-D
%   findedgedof      - highlights some dofs associated to edges
%   findedgedof3     - highlights some dofs associated to edges
%   findelem         - highlights some elements
%   findelem3        - highlights the elements in the range
%   findnode         - highlights nodes in certain range.
%   findnode3        - highlights some nodes for 3-D meshes.
%   findnodedof      - highlights dof associated to nodes.
%   findquadelem     - highlights some elements
%   hg               - the interface of command hg in matlab
%   includepath      - inlcudes pathes
%   latexmatrix      - output the matrix in latex formate.
%   latexmesh        - draw a mesh using pstrick
%   latextable       - output the matrix in table formate in latex.
%   showboundary     - boundary mesh plot
%   showboundary3    - boundary surface mesh plot
%   showcoloring     - plot coloring vertices set
%   showgraph        - displays a planar graph
%   showgraphmatrix  - displays a planar graph
%   showmatrix       - displays a matrix.
%   showmesh         - displays a triangular mesh in 2-D.
%   showmesh3        - displays a tetrahedron mesh in 3-D.
%   showmeshdensity  - show the density of the mesh
%   showmeshquality  - show the qualify of the mesh
%   showquadmesh     - displays a QUAD mesh in 2-D
%   showquadsolution - plots the solution u on a quad mesh in 2-D.
%   showrate         - rate of an err sequence
%   showrate2        - rate of two error sequences
%   showrate3        - rate of two error sequences
%   showresult       - display the mesh and the solution
%   showresult3      - display the mesh and the solution in 3-D
%   showsolution     - plots the solution u on a triangular mesh in 2-D.
%   showsolution3    - plots the solution u on a tetrahedron mesh in 3-D.
%   spsub2ind        - SPARSESUB2IND find values of a sparse matrix

repo: ifem


text:
function r = vertcat(c1,c2)

%tstoolbox/@core/vertcat
%   Syntax:
%     * r = vertcat(c1,c2)
%
%   Input Arguments:
%     * c1,c2 - core objects
%
%   catenate two timeseries verticaly
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

r = core([data(c1) ; data(c2)]);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [hd_estimates]=hdols123(const,exo,betahat,k,n,p,D,m,T,X,Y,data_exo)

contributors = n + const + length(exo) + 1; %variables + constant + exogenous + initial conditions
hd_estimates=cell(contributors+2,n); %shocks+constant+initial values+exogenous+unexplained+to be explained by shocks only
hd_estimates2=cell(contributors+2,n); %contributors + to be explained by shocks only (unexplained part, other part)

identified=n;

%% compute the historical decomposition, if its not yet computed
%===============================================
%%%%% could take this part from previous computations if possible
Bfull=reshape(betahat,k,n);                    %get the Bfull matrix
B=Bfull(1:n*p,:);                           %drop the coefficients for all exogenous variables from the matrix
B_comp = [B'; eye(n*(p-1)) zeros(n*(p-1),n)];%put into companion form
EPS=Y-X*Bfull;                               %get reduced form residuals
ETA=(D\EPS');                              %get structural shocks

%% Compute historical decompositions
%===============================================

% Contribution of each shock
    aux_D = zeros(n*p,n); %auxilary D matrix, that is consistent with companion matrix
    aux_D(1:n,:) = D; %set the first entrys equal to Dinv
    Selec = [eye(n) zeros(n,(p-1)*n)]; %selection matrix
    HDestimates_store = zeros(p*n,T+1,n); %cell aray to store results
    HDestimates = zeros(n,T+1,n);
    for j=1:n % for each variable %%%%%
        ETA_comp = zeros(n,T+1); % structural shock matrix that is consistent with companion form
        if j <= identified % if j is an identified shock
        ETA_comp(j,2:end) = ETA(j,:); % fill in the entry for shock n, leave the first entry blank
        end
        for i = 2:T+1
            HDestimates_store(:,i,j)=aux_D*ETA_comp(:,i) + B_comp*HDestimates_store(:,i-1,j); %recursively sum over shock impulse at period i on variable j and the previous period
            HDestimates(:,i,j)=Selec*HDestimates_store(:,i,j); %select the entry corresponding to the current period
        end
    end

%% contribution of the initial values
    HDinitial_storage   = zeros(p*n,T+1);
    HDinitial_estimates = zeros(n, T+1);
    Xnoexo = X(:,1:n*p);
    HDinitial_storage(:,1) = Xnoexo(1,:)'; %set the initial values to the first row of X (n*P)+exo
    HDinitial_estimates(:,1) = Selec*HDinitial_storage(:,1); %select the initial values for the first n variables (i.e. the values at Y_{t-1}
    for i = 2:T+1 %loop over periods and compute the impact of the initial conditions recursively
        HDinitial_storage(:,i) = B_comp*HDinitial_storage(:,i-1); %compute the impact of those values (which in principle consist of past shocks)
        HDinitial_estimates(:,i) = Selec*HDinitial_storage(:,i);
    end

if const==1
%%  Contribution of the Constant
    HDconstant_storage = zeros(p*n,T+1);
    HDconstant_estimates = zeros(n, T+1);
    Coefficients = zeros(p*n,1);
        Coefficients(1:n,:) = Bfull(n*p+1,:);
        for i = 2:T+1 %loop over periods
            HDconstant_storage(:,i)=Coefficients+B_comp*HDconstant_storage(:,i-1);
            HDconstant_estimates(:,i)=Selec*HDconstant_storage(:,i);
        end
end
%% Contribution of exogenous variables
 if m > 1
    HDexo_storage = zeros(p*n,T+1);
    HDexo_estimates = zeros(n,T+1);
    Coefficients_exo = zeros(p*n,(m-1)*(1));
    data_exocut=data_exo(p+1:end,:); %cut initial conditions from exogenous
    Coefficients_exo(1:n,:) = Bfull(n*p+const+1:end,:)'; %get the corresponding coefficients
        for i = 2:T+1
            HDexo_storage(:,i)=Coefficients_exo*data_exocut(i-1,:)'+B_comp*HDexo_storage(:,i-1);
            HDexo_estimates(:,i)=Selec*HDexo_storage(:,i);
        end
 end

%% put these values into the corresponding cell for hd_estimates such that
% for variable x (hd_estimates(x,n+1)) = HDinitial_estimates(x,:)
%reorganize storage
        for jj=1:n %for variables
            for kk=1:T+1 %for periods
                for ii=1:n %for shock contributions
                    hd_estimates2{ii,jj}(1,kk)=HDestimates(jj,kk,ii);
                end
                hd_estimates2{n+1,jj}(1,kk)=HDinitial_estimates(jj,kk);
                if const==1
                hd_estimates2{n+2,jj}(1,kk)=HDconstant_estimates(jj,kk);
                end
                if m>1
                hd_estimates2{n+3,jj}(1,kk)=HDexo_estimates(jj,kk);
                end
            end
        end

 HDsum = zeros(T+1,n); %if we sum over all variables this should give Y
 for jj=1:n %loop over variables (columns)
 sumvariable = zeros(1,T+1);
 for kk=1:T+1 %loop over periods
     sumperiod=0;
  for ii=1:contributors %loop over contributors (rows)
      sumperiod=sumperiod+hd_estimates2{ii,jj}(1,kk);
  end
  sumvariable(1,kk)=sumperiod;
 end
 HDsum(:,jj)=sumvariable(1,:);  %%%%% this should be Y?
 end

 %% determine the unexplained part (if model is not fully identified) %%%%% does that mean in IRFt5???
 aux=zeros(1,n);
 unexplained = Y-HDsum(2:end,:);
 unexplained = [aux; unexplained];
for jj=1:n
    hd_estimates2{contributors+1,jj}=unexplained(:,jj)';
end

%% finally substract the sum of the contribution of the (residual?)
%exogenous, constant,initial conditions from Y to get the
%part that was left to be explained by the shocks (for plotting reasons)
 Exosum = zeros(T+1,n); %if we sum over all variables this should give Y
 for jj=1:n %loop over variables (columns)
 sumvariable = zeros(1,T+1);
 for kk=1:T+1 %loop over periods
     sumperiod=0;
  for ii=n+1:contributors %loop over contributors (rows)
      sumperiod = sumperiod+hd_estimates2{ii,jj}(1,kk);
  end
  sumvariable(1,kk)=sumperiod;
 end
 Exosum(:,jj)=sumvariable(1,:);
 end

 %% determine the part that was left to be explained by the shocks
 aux = zeros(1,n);
 tobeexplained = [aux; Y-Exosum(2:end,:)];
for jj=1:n
    hd_estimates2{contributors+2,jj}=tobeexplained(:,jj)';
end

%drop the initial entry for each cell in hd_estimates
for jj=1:n % over variables
    for ii=1:contributors+2 %all contributions
        hd_estimates{ii,jj}=hd_estimates2{ii,jj}(2:end);
    end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function r = acoth(a)
%ACOTH        Hessian (elementwise) inverse hyperbolic cotangent
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = acoth(a.x);
    f = 1 / ( 1 - sqr(a.x) );
    r.dx = a.dx * f;
    r.hx = ( a.hx + reshape( (a.x*r.dx) * a.dx.' , size(a.hx) ) ) * f;

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = acoth(a.x);
    if issparse(a.hx)               % input sparse

      ax = 1 ./ ( 1 - sqr(full(a.x(:))) );
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = a.x(:).*ax;
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,r.dx);
      end
      [ia,ja,sa] = find(a.hx);      % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = acoth(a.x);
      ax = a.x(:).';
      f = 1 ./ ( 1 - sqr(ax) );
      f = f(ones(N*N,1),:);
      r.dx = a.dx .* f(1:N,:);
      adx = repmat(ax,N,1) .* r.dx;
      r.hx = ( a.hx + adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:) ) .* f;

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%DEMO_PASSGP Demonstration of PASS-GP method for GP classification
%
%  Description
%    Here we demonstrate PASS-GP method for Gaussian Processes
%    classification. Data used is 2-dimensional toy data with Gaussian
%    bumbs defining classes. We demonstrate with both fixed and not fixed
%    sizes of active set for PASS-GP.
%
%    PASS-GP uses a predictive active set selection method by Henao &
%    Winther (2012) to select a subset of training data to be used for
%    inference in classification problems.
%
%  Reference:
%    Ricardo Henao & Ole Winther (2012). Predictive active set selection
%    methods for Gaussian processes. Neurocomputing 80 (2012), 10-18.
%
%  See also PASSGP
%
% Copyright (c) 2013 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Generate toy data
prevstream=setrandstream(0);
[x1,x2]=meshgrid(-5:0.1:5,-5:0.1:5);
x=[x1(:) x2(:)]; x=x(randperm(size(x,1),3000),:);
y=2.*mnorm_pdf(x, [0 0], [0.5 0;0 0.5]) + mnorm_pdf(x, [3 3], [0.5 0;0 0.5]) + mnorm_pdf(x, [-3 -3], [0.5 0;0 0.5]);
y=y+mnorm_pdf(x, [3 -3], [0.5 0;0 0.5])+mnorm_pdf(x, [-3 3], [0.5 0;0 0.5]);
y=y+0.03.*randn(size(y,1),1);
y(y>0.15)=1; y(y<=0.15)=-1;

[xt1, xt2]=meshgrid(-5:0.23:5,-5:0.23:5);
xt=[xt1(:) xt2(:)];
yt=ones(size(xt,1),1);

[n, nin] = size(x);

% Define covariance and likelihood functions and create the model
gpcf = gpcf_sexp();
lik=lik_probit();
gp=gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6);

opt=optimset('TolX',1e-3,'TolFun',1e-3,'display','on');
w0=gp_pak(gp);

% fPASS-GP with fixed size of 800 points in active set and data divided to
% 10 subsets with 4 sweeps over data
start=tic;[gp, indA]=passgp(gp, x, y, 'opt', opt, 'npass', 4, 'ninit', 800, 'nsub', 10, 'display', 'on', 'fixed', 'on', 'pexc', 0.1, 'optimn', 2);time=toc(start);
tt=time;
[Eft, Varft, lpyt, Eyt, Varyt]=gp_pred(gp, x(indA,:), y(indA,:), xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (fpass-gp)')
param=gp_pak(gp);

% PASS-GP with inclusion threshold 0.65, deletion threshold 0.99, intial
% size of 400 points in active set and 3 sweeps over data.
gp=gp_unpak(gp,w0);
start=tic;[gp, indA2]=passgp(gp, x, y, 'opt', opt, 'pinc', 0.65, 'pdel', 0.99, 'npass', 3, 'ninit', 400, 'nsub', 10, 'display', 'on', 'optimn', 2);time=toc(start);
tt2=time;
[Eft2, Varft2, lpyt2, Eyt2, Varyt2]=gp_pred(gp, x(indA2,:), y(indA2,:), xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt2),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (pass-gp)')
param2=gp_pak(gp);

% Full Gaussian process for comparison
gp=gp_unpak(gp,w0);
opt.Display='iter';
start=tic;gp=gp_optim(gp,x,y,'opt',opt);tt3=toc;
[Eft3, Varft3, lpyt3, Eyt3, Varyt3]=gp_pred(gp, x, y, xt, 'yt', yt);
figure, [cc,hh]=contour(reshape(xt(:,1),size(xt1,1), size(xt1,1)), reshape(xt(:,2),size(xt1,1), size(xt1,1)), reshape(exp(lpyt3),size(xt1,1), size(xt1,1)), [0.1 0.9]);
clabel(cc,hh); title('Pr(y==1) (full gp)')

% Display some statistics
mlpd_fpassgp=mean(mean(lpyt,2))
time_fpassgp=mean(tt)
mlpd_passgp=mean(mean(lpyt2,2))
time_passgp=mean(tt2)
mlpd_full=mean(lpyt3)
time_full=tt3

% Plot data and active sets for both methods
figure(4), subplot(1,2,1),  plot(x(y==1,1),x(y==1,2),'or',x(y==-1,1),x(y==-1,2),'ob');
hold all; plot(x(indA,1), x(indA,2), '.k'); title('Data and active set (fpass-gp)')
legend('y=1', 'y=-1', 'Active set for fpass-gp');
subplot(1,2,2),  plot(x(y==1,1),x(y==1,2),'or',x(y==-1,1),x(y==-1,2),'ob');
hold all; plot(x(indA2,1), x(indA2,2), '.k'); title('Data and active set (pass-gp)')
legend('y=1', 'y=-1', 'Active set for pass-gp');
setrandstream(prevstream);

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function acc = LRSDL_pred_LC(Y, D, D0, CoefMM0, opts, label_test)
    nClasses = size(CoefMM0, 2);
    N = size(Y,2);
    acc = [];
    optsX.max_iter = 500;
    for lambda = [0.001, 0.005, 0.01, 0.05]
        E = zeros(nClasses, size(Y,2));
        %% ========= LC ==============================
        for c = 1: nClasses
            Di = [get_block_col(D, c, opts.D_range) D0];
            [Xi, ~] = lasso_fista(Y, Di, [], lambda, optsX);
            R = Y - Di*Xi;
            E(c,:) = 0.5*sum(R.^2,1) + lambda*sum(abs(Xi),1);
        end
        [~, pred] = min(E);
        aaaa = double(sum(pred == label_test))/N;
        acc = [acc aaaa];
        fprintf('lambda = %.4f, acc: %f\n', lambda, aaaa);
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function cost = getCost(x)
    % cost = exp(x);
    cost = x^2;
    cost = x;
end

author: UMich-BipedLab


text:
function ChanLoc = channel_project_scalp(Vertices, ChanLoc)
% CHANNEL_PROJECT_SCALP: Project EEG electrodes on a scalp surface,
% radially from the center of mass of the scalp surface points.
%
% INPUT:
%     - Vertices : [Mx3] positions of the scalp vertices
%     - ChanLoc  : [Nx3] positions of the EEG electrodes

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2014

% Center the surface on its center of mass
center = mean(Vertices, 1);
Vertices = bst_bsxfun(@minus, Vertices, center);
% Parametrize the surface
p   = .2;
th  = -pi-p   : 0.01 : pi+p;
phi = -pi/2-p : 0.01 : pi/2+p;
rVertices = tess_parametrize(Vertices, th, phi);

% Process each sensor
for iChan = 1:size(ChanLoc,1)
    % Get the closest surface from the point
    c = ChanLoc(iChan,:);
    % Center electrode
    c = c - center;
    % Convert in spherical coordinates
    [c_th,c_phi,c_r] = cart2sph(c(1), c(2), c(3));
    % Interpolate
    c_r = interp2(th, phi, rVertices, c_th, c_phi);
    % Project back in cartesian coordinates
    [c(1),c(2),c(3)] = sph2cart(c_th, c_phi, c_r);
    % Restore initial origin
    c = c + center;
    ChanLoc(iChan,:) = c;
end

author: brainstorm-tools


text:
function value = p27_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P27_F evaluates the integrand for problem 27.
%
%  Dimension:
%
%    N arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    The integral depends on a parameter R and vector C(1:N).
%
%    R defaults to 0.3.
%
%    The reference suggests choosing C at random in [0,1]
%    and then multiplying by the normalizing factor (25/N).
%    C(1:N) defaults to 1/N.
%
%    To get or set R, call P27_R8.
%    To get or set C, call P27_R8VEC.
%
%  Integrand:
%
%    cos ( 2 * pi * R + sum ( c(1:n) * x(1:n) ) )
%
%  Exact Integral:
%
%    2^N * cos ( 2 * pi * R + 0.5 * sum ( c(1:n) ) )
%      * product ( sin ( 0.5 * c(1:n) ) / c(1:n) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    [Integral #1]
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    D Reidel, 1987, pages 337-340,
%    LC: QA299.3.N38.
%
%    Thomas Patterson,
%    [Integral #5],
%    On the Construction of a Practical Ermakov-Zolotukhin
%    Multiple Integrator,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    D. Reidel, 1987, pages 269-290,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  r = 0.0;
  r = p27_r8 ( 'G', 'R', r );

  c = [];
  c = p27_r8vec ( 'G', 'C', dim_num, c );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    arg = 2.0 * pi * r + c(1:dim_num,1)' * x(1:dim_num,point);
    value(point) = cos ( arg );

  end

  p27_i4 ( 'I', '#', point_num );

  return
end


text:
function ht = my_xticklabels(varargin)

%MY_XTICKLABELS replaces XTickLabels with "normal" texts
%   accepting multiline texts and TEX interpreting
%   and shrinks the axis to fit the texts in the window
%
%    ht = my_xticklabels(Ha, xtickpos, xtickstring)
% or
%    ht = my_xticklabels(xtickpos, xtickstring)
%
%  in:    xtickpos     XTick positions [N*1]
%        xtickstring   Strings to use as labels {N*1} cell of cells
%
% Examples:
% plot(randn(20,1))
% xtl = {{'one';'two';'three'} '\alpha' {'\beta';'\gamma'}};
% h = my_xticklabels(gca,[1 10 18],xtl);
% % vertical
% h = my_xticklabels([1 10 18],xtl, ...
%     'Rotation',-90, ...
%     'VerticalAlignment','middle', ...
%     'HorizontalAlignment','left');

% Pekka Kumpulainen 12.2.2008

textopts = {};
if length(varargin{1})==1 && ...
        ishandle(varargin{1}) && ...
        strcmpi(get(varargin{1},'Type'),'axes');
    Ha = varargin{1};
    xtickpos = varargin{2};
    xtickstring = varargin{3};
    if nargin > 3
        textopts = varargin(4:end);
    end
else
    Ha = gca;
    Hfig = get(Ha,'Parent');
    xtickpos = varargin{1};
    xtickstring = varargin{2};
    if nargin > 2
        textopts = varargin(3:end);
    end
end

set(Ha,'XTick',xtickpos, 'XTickLabel','')
h_olds = findobj(Ha, 'Tag', 'MUXTL');
if ~isempty(h_olds)
    delete(h_olds)
end

%% Make XTickLabels
NTick = length(xtickpos);
Ybot = min(get(gca,'YLim'));
ht = zeros(NTick,1);
for ii = 1:NTick
    ht(ii) = text('String',xtickstring{ii}, ...
        'Units','data', ...
        'VerticalAlignment', 'top', ...
        'HorizontalAlignment', 'center ', ...
        'Position',[xtickpos(ii) Ybot], ...
        'Tag','MUXTL');
end
if ~isempty(textopts)
    set(ht,textopts{:})
end

%% squeeze axis if needed

set(Ha,'Units','pixels')
Axpos = get(Ha,'Position');
% set(Hfig,'Units','pixels')
% Figpos = get(Hfig,'Position');

set(ht,'Units','pixels')
TickExt = zeros(NTick,4);
for ii = 1:NTick
    TickExt(ii,:) = get(ht(ii),'Extent');
end

needmove = -(Axpos(2) + min(TickExt(:,2)));

if needmove>0;
    Axpos(2) = Axpos(2)+needmove+2;
    Axpos(4) = Axpos(4)-needmove+2;
    set(Ha,'Position',Axpos);
end

set(Ha,'Units','normalized')
set(ht,'Units','normalized')

author: Sable

repo: mcbench-benchmarks


text:
function dist = segment_point_dist_2d ( p1, p2, p )

%*****************************************************************************80
%
%% SEGMENT_POINT_DIST_2D: distance ( line segment, point ) in 2D.
%
%  Discussion:
%
%    A line segment is the finite portion of a line that lies between
%    two points.
%
%    The nearest point will satisfy the condition
%
%      PN = (1-T) * P1 + T * P2.
%
%    T will always be between 0 and 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), the endpoints of the line segment.
%
%    Input, real P(2,1), the point whose nearest neighbor on the line
%    segment is to be determined.
%
%    Output, real DIST, the distance from the point to the line segment.
%

%
%  If the line segment is actually a point, then the answer is easy.
%
  if ( p1(1:2,1) == p2(1:2,1) )

    t = 0.0;

  else

    bot = sum ( ( p2(1:2,1) - p1(1:2,1) ).^2 );

    t = ( p(1:2,1) - p1(1:2,1) )' * ( p2(1:2,1) - p1(1:2,1) ) / bot;

    t = max ( t, 0.0 );
    t = min ( t, 1.0 );

  end

  pn(1:2,1) = p1(1:2,1) + t * ( p2(1:2,1) - p1(1:2,1) );

  dist = sqrt ( sum ( ( pn(1:2,1) - p(1:2,1) ).^2 ) );

  return
end


text:
function inplanes = rtGetInplanes(examNum,ipSeries,doCrop);
%
% inplanes = rtGetInplanes(examNum,inplaneSeriesNum,[doCrop]);
%
% For pseudo real-time, copy over the DICOM files
% from the scanner for the selected inplane
% anatomical series (to the session Raw/Inplane dir),
% then read them in as in InitAnatomy and save/return
% the anat struct.
%
% ras 04/05
if ieNotDefined('examNum')
    examNum = input('Enter exam #: ');
end
if ieNotDefined('ipSeries')
    ipSeries = input('Enter series for Inplanes: ');
end
if ieNotDefined('doCrop')
    doCrop = input('Crop Inplanes? [1 for yes, 0 for no]: ');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run a unix command to ssh over to    %
% the control machine (lcmr3) and get  %
% the DICOM images:                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first, the command to be executed at lcmr3:
xferCmd = sprintf(['/export/home/sdc/bin/lx_ximg ',...
                   'E%iS%iIall -d /usr/g/mrraw'],...
                   examNum,ipSeries);
% now, the ssh command to call this remotely
cmd = sprintf('ssh sdc@lcmr3 %s',xferCmd);

% execute
unix(cmd)

% now, make sure there's a local Raw/Inplane dir:
if ~exist(fullfile(pwd,'Raw'),'dir')
    mkdir Raw
    disp('Made Raw')
end

if ~exist(fullfile(pwd,'Raw','Anatomy'),'dir')
    mkdir Raw Anatomy
    disp('Made Raw/Anatomy')
end

if ~exist(fullfile(pwd,'Raw','Anatomy','Inplane'),'dir')
    mkdir Raw/Anatomy Inplane
    disp('Made Raw/Anatomy/Inplane')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Xfer the files to the local Raw dir %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the file names
pattern = sprintf('/lcmr3/mrraw/E%iS%i*.MR.dcm',examNum,ipSeries);
w = dir(pattern);
nFiles = length(w);
if nFiles==0
    error('No matching files found!')
end
for i = 1:nFiles
    srcFile = sprintf('E%iS%iI%i.MR.dcm',examNum,ipSeries,i);
    tgtFile = sprintf('Raw/Anatomy/Inplane/I%03d.dcm',i);
    cmd = sprintf('cp /lcmr3/mrraw/%s %s',srcFile,tgtFile);
    unix(cmd);
    fprintf('Got %s \n',tgtFile);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now get the inplanes for mrVista %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[anat, inplanes] = InitAnatomy(pwd,'Raw',doCrop);

if doCrop==0
    %%% GUM
    inplanes.crop = [1 1; inplanes.fullSize];
    inplanes.cropSize = [inplanes.fullSize];
else
    [anat, inplanes] = CropInplanes('Raw',anat,inplanes,[1 1]);
end

% save the anat image locally
if ~exist(fullfile(pwd,'Inplane'),'dir')
    disp('Making Inplane directory');
    mkdir Inplane
end
save('Inplane/anat','anat');
disp('Saved Inplane/anat.mat file.')

% update mrSESSION
if exist('mrSESSION.mat','file')
    load mrSESSION
    mrSESSION.inplanes = inplanes;
    save mrSESSION mrSESSION -append
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ samplingMask ] = readSamplingMaskFromFile(sFilename)
%This function reads out a sampling pattern from a .txt-file

if(nargin < 1)
    file = fopen('samplingPattern.txt','r');
else
    file = fopen(sFilename,'r');
end

endOfDoc = false;
i = 0;

while(~endOfDoc)
    actualLine = fgetl(file);

    if(actualLine == -1)
        endOfDoc = true;
        break;
    end
    i = i + 1;
    lines{i} = actualLine;

    for j=1:numel(lines{i})
        actualChar = lines{i}(j);

        if(actualChar == '-')
            samplingMask(i, j) = 0;
        elseif (actualChar == 'x' || actualChar == '*' || actualChar == 'O' || actualChar == '0')
            samplingMask(i, j) = 1;
        else
            error('not working, because there is an unknown charakter');
        end
    end
end

if(exist('drawMask','file'))
    drawMask(samplingMask);
end

fclose(file);

end

author: thomaskuestner


text:
% SCRIPT TEST FOR THE UR10 ROBOT KINEMATICS

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

close all

robot = load_robot('UR', 'UR10');
% adjust 3D view as desired
adjust_view(robot)

% test both solutions: based on the transpose an on the Moore-Penrose
robot.inversekinematic_fn = 'inverse_kinematics_jacobian(robot, T, q)';
q0 = [0.1 -pi/2 -pi/2 0.1 0.1 0.1]';
%q0 = [0.2  0.2 0.2 0.2 0.2 0.2]';
%q0 = [0.2  0.2 0.2 0.2 0.2 0.2]';

T0 = directkinematic(robot, q0)

fprintf('\nSimple test: try to reach T')
T_target = [-0.3390   -0.3390    0.8776    0.2447;
            0.8469    0.2962    0.4416   -0.4373;
            -0.4096    0.8929    0.1867    1.2568;
              0         0         0    1.0000];

qinv = inversekinematic(robot, T_target, q0)

T_reached = directkinematic(robot, qinv)
'diff T-T_target'
T_reached-T_target
% Plot manipulatiliby ellipse
drawrobot3d(robot, qinv)

repo: ARTE


text:
function [f,nodes] = evalSectionsEquispacedFFT(SO3F,varargin)
% Evaluate an @SO3FunHarmonic on @ODFSections by using an equispaced grid
% along the other 2 Euler angles
%     $$(\alpha_a,\beta_b,\gamma_c) = (\frac{2\pi a}{H_1},\frac{\pi b}{H_2-1},\gamma_c)$$
% where $a=0,...,H_1-1$, $b=0,...,H_2-1$ and $c=0,...,S_{num}$.
%
% Therefore we transform the SO(3) Fourier series to an usual Fourier series
% equivalent as in the function <SO3FunHarmonic.evalV2.html |evalV2|>.
% But we use an 2-variate equispaced FFT instead of the NFFT analogously to
% <SO3FunHarmonic.evalEquispacedFFT.html |evalEquispacedFFT|>
%
% Syntax
%   f = evalSectionsEquispacedFFT(SO3F)
%   f = evalSectionsEquispacedFFT(SO3F,'resolution',2.5*degree)
%   [f,nodes] = evalSectionsEquispacedFFT(SO3F,'resolution',[2*degree,2.5*degree])
%   f = evalSectionsEquispacedFFT(SO3F,oS)
%
% Input
%  SO3F - @SO3FunHarmonic
%  oS - @ODFSection (phi2,gamma,phi1)
%
% Options
%  'resolution' - shape constant along Euler angles. (default = 2.5°)
%
% Output
%  nodes - @orientation
%  f - values at this grid points
%
% See also
% SO3FunHarmonic/evalV2 SO3FunHarmonic/evalEquispacedFFT SO3FunHarmonic/eval

N = SO3F.bandwidth;

if isa(varargin{1},'ODFSections'), oS = varargin{1}; else, oS =[]; end
[a_max,b_max,g_max] = fundamentalRegionEuler(SO3F.CS,SO3F.SS,varargin{:});
if isempty(oS)
  gamma = mod((0:5)*g_max/6+pi/2,2*pi);
  shift = 1;
elseif isa(oS,'gammaSections')
  gamma = oS.gamma;
  shift = 0;
elseif isa(oS,'phi2Sections')
  gamma = mod(oS.phi2+pi/2,2*pi);
  shift = 1;
elseif isa(oS,'phi1Sections')
  % phi1-sections are analogous to alpha sections. We transform SO3F(R(a,b,g))
  % to SO3F(R(g,b,a)).
  fhat = SO3F.fhat;
  for n=0:N
    ind = deg2dim(n)+1 : deg2dim(n+1);
    fhat(ind) = (-1).^((1:2*n+1)+(1:2*n+1).') .* reshape(fhat(ind),2*n+1,2*n+1).';
  end
  a_max = g_max;
  SO3F = SO3FunHarmonic(fhat,SO3F.SS,SO3F.CS);
  gamma = mod(oS.phi1-pi/2,2*pi);
  shift = -1;
end

% get grid size along 1st and 2nd Euler angles
res = get_option(varargin,'resolution',2.5*degree);
if length(res)==1
  res = [1,1]*res;
end
% [a_max,b_max,~] = fundamentalRegionEuler(SO3F.CS,SO3F.SS,varargin{:});
s = ceil([a_max,b_max]./res);
res = [a_max,b_max]./s;
H = 2*pi./res;
if any(s<=1)
  error('The resolution is to big.')
end

shiftGrid = get_option(varargin,'shiftGrid',[0,0]);
if check_option(varargin,'shiftGrid')
  s = s-1;
end

% compute ghat
if SO3F.isReal

  % compute ghat -> k x j x l
  % with  k = -N:N
  %       j = -N:N    -> use ghat(k,-j,l) = (-1)^(k+l) * ghat(k,-j,l)   (*)
  %       l =  0:N    -> use ghat(-k,-j,-l) = conj(ghat(k,j,l))        (**)
  % flags: 2^0 -> use L_2-normalized Wigner-D functions
  %        2^2 -> fhat are the fourier coefficients of a real valued function
  %        2^4 -> use right and left symmetry
  flags = 2^0+2^2+2^4;
  sym = [min(SO3F.SRight.multiplicityPerpZ,2),SO3F.SRight.multiplicityZ,...
         min(SO3F.SLeft.multiplicityPerpZ,2),SO3F.SLeft.multiplicityZ];
  ghat = representationbased_coefficient_transform(N,SO3F.fhat,flags,sym);
  ghat = permute(ghat,[3,2,1]);
  % correct ghat by i^(-k+l)
  z = - (0:N)'+ reshape(-N:N,1,1,[]) + shift * (0:N)';
  ghat = ghat.*(1i).^z;

  if check_option(varargin,'shiftGrid')
    ghat = ghat.*exp(-1i * (shiftGrid(1)*(0:N)' + shiftGrid(2)*(-N:N)) );
  end

else

  % compute ghat
  % flags: 2^0 -> use L_2-normalized Wigner-D functions
  %        2^4 -> use right and left symmetry
  flags = 2^0+2^4;
  sym = [min(SO3F.SRight.multiplicityPerpZ,2),SO3F.SRight.multiplicityZ,...
         min(SO3F.SLeft.multiplicityPerpZ,2),SO3F.SLeft.multiplicityZ];
  ghat = representationbased_coefficient_transform(N,SO3F.fhat,flags,sym);
  ghat = permute(ghat,[3,2,1]);
  % correct ghat by i^(-k+l)
  z = - (-N:N)' + reshape(-N:N,1,1,[]) + shift * (-N:N)';
  ghat = ghat.*(1i).^z;

  if check_option(varargin,'shiftGrid')
    ghat = ghat.*exp(-1i * (shiftGrid(1)*(-N:N)' + shiftGrid(2)*(-N:N)) );
  end

end

% use rotational symmetries around Z-axis to speed up (cut zeros in ghat)
SRightZ = SO3F.SRight.multiplicityZ;
SLeftZ = SO3F.SLeft.multiplicityZ;
H(1) = H(1) / SLeftZ;
if SLeftZ>1 || SRightZ>1
  if SO3F.isReal
    ind1 =  (0:SLeftZ:N)+1;
  else
    ind1 =  [-flip(SLeftZ:SLeftZ:N),(0:SLeftZ:N)] + (N+1);
  end
  ind3 =  [-flip(SRightZ:SRightZ:N),(0:SRightZ:N)] + (N+1);
  ghat = ghat(ind1,:,ind3);
end

% For small H we go through a smaller FFT(H) several times. Hence we reduce
% the size of the fourier coefficient matrix ghat by adding the coefficients
% with same complex exponentials.
sz = size(ghat,1,2);
if any(H<sz)
  dim = ceil(sz./H);
  A = zeros([dim.*H,size(ghat,3)]);
  A(1:size(ghat,1),1:2*N+1,1:size(ghat,3)) = ghat;
  % Note that H(1) and H(2) should be bigger than 1 to avoid errors by squeezing
  ghat = squeeze(sum(reshape(A,H(1),dim(1),H(2),dim(2),[]),[2,4]));
  clear A;
end

% Do 2-variate fft for all Euler angles gamma
f = zeros([s+1,length(gamma)]);
for g = 1:length(gamma)
  exponentials = exp( -1i*gamma(g) * reshape([-flip(SRightZ:SRightZ:N),(0:SRightZ:N)],1,1,[]) );
  A = sum(ghat.*exponentials,3);
  % fft
  f2 = fft2(A,H(1),H(2));
  % cut beta to [0,pi]
  if s(1) == H(1)
    f2 = f2([1:s(1),1],1:s(2)+1);
  else
    f2 = f2(1:s(1)+1,1:s(2)+1);
  end
  % shift the summation of fft from [-|_N/r_|:|_N/r_|]x[-N:N]x[-|_N/s_|:|_N/s_|] to
  % [0:2*|_N/r_|]x[0:2N]x[0:2*|_N/s_|]. With r- & s-fold rotational symmetry
  % around Z-axis and |_ ... _| denotes the round off operator.
  if SO3F.isReal
    z = (0:s(2)) * (N/H(2));    % oder fftshift
    f(:,:,g) = 2*real( exp(2i*pi*z) .* f2 );
  else
    z = (0:s(1)).' * (floor(N/SLeftZ)/H(1)) + (0:s(2)) * (N/H(2));
    f(:,:,g) = exp(2i*pi*z) .* f2;
  end
end

% compute the corresponding nodes
if nargout>1
  if check_option(varargin,'shiftGrid')
    grid = combvec((0:s(1))*res(1)-shift*pi/2,(0:s(2))*res(2),gamma);
    grid = grid + [shiftGrid,0];
  else
    grid = combvec((0:s(1))*res(1)-shift*pi/2,(0:s(2))*res(2),gamma);
  end
  if isa(oS,'phi1Sections')
    grid = grid(:,[3,2,1]);
    grid = orientation.byEuler(grid,'nfft',SO3F.SS,SO3F.CS);
  else
    grid = orientation.byEuler(grid,'nfft',SO3F.CS,SO3F.SS);
  end
  nodes = reshape(grid,size(f));
end

end

function v = combvec(v1,v2,v3)
% This function does the same as the combvec() function from Digital
% Processing Toolbox and additionaly transpose the output matrix.
% It is a simple version for the special case of three input vectors.
  l1 = length(v1);
  l2 = length(v2);
  l3 = length(v3);

  V1 = repmat(v1.',l2*l3,1);
  V2 = repmat(v2,l1,l3);
  V3 = repmat(v3,l1*l2,1);

  v = [V1,V2(:),V3(:)];
end

author: mtex-toolbox

repo: mtex

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


text:
function [Cbe_new, Ve_new, ecef_new]=strapdown_ecef_dcm(Cbe, Ve, ecef, a, w, dt)
%Gravity (most time consuming part of ecef implementations)
Llh=ecef2geo_v000(ecef,0);
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(Llh);
Cne=pos2Cne_v000(Llh(1), Llh(2));
ge=-Cne(:,3)*g;

%Update attitude
rot=w*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_a=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);

rot=-([0;0;WIE_E])*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_b=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);

Cbe_new=mx_b*Cbe*mx_a;

%%Update Velocity
%Update Vel
vel_inc1=Cbe*a*dt;
vel_inc2=(-ge+2*cross(Ve,[0;0;WIE_E]))*dt;
Ve_new=Ve+vel_inc1+vel_inc2;

%Update_pos
ecef_new=ecef+Ve*dt;

repo: ekfmonoslam


text:
function [W, p, t] = build8(im, angles)
t_s = toc;
%Pad image
[im_pad,D] = impad(im,'diag',0);
sz = size(im);
m = sz(1); %rows
n = sz(2); %cols
n_proj = length(angles);
N = m*n; %number of unknown variables
M = D*n_proj; %number of equations
%Weighting Factor Matrix
try
    W = zeros(M,N);
catch expr
    fprintf([expr.message '\nGenerates a sparse.\r'])
    W = sparse(M,N);
end
%proj_info = zeros(M,2); %[theta No.]
p = zeros(M,1); %projection vector
ix = false(M,1);
rc = [m;n]/2+0.5;
fprintf('Building Weight Matrix...\r')
[x,y] = ind2sub(sz,1:N);
%xy = [x;y]; %x,y coordinate
%x,y coordinate with respect to rotation center
xy_c = [x;y]-repmat(rc,1,N);
d = 1; %side length of square
A0 = d^2;
H = @heaviside;
try
    A = zeros(D,N);
catch expr
    fprintf(['Building A...\n' ...
        expr.message '\nGenerating a sparse...\r'])
    A = sparse(D,N);
end
for kp = 1:n_proj
    A(:) = 0; %#
    im_rot = imrotate(im_pad,-angles(kp),'bilinear','crop');
    pvec = sum(im_rot,1);
    t = -angles(kp)/180*pi;
    R = [cos(t) -sin(t);sin(t) cos(t)];
    fprintf('\nAngle No.%d(%d Degree)\r',kp,angles(kp))
    xy_rot = R*xy_c+D/2+0.5;
    idx = round(xy_rot(2,:));
    ixM = D*(kp-1)+idx; %corresponding indice in W and p matrix
    %Calculate Area
    x1 = d*min(abs([sin(t) cos(t)]))+eps;
    x2 = d*max(abs([sin(t) cos(t)]))-eps;
    l = x1+x2;
    h = d^2/x2;
    %A0 = d^2;
    A1 = x1*h/2;
    G1 = @(x) H(x)-H(x-x1);
    G2 = @(x) H(x-x1)-H(x-x2);
    G3 = @(x) H(x-x2)-H(x-l);
    f_A = @(x) (x/x1).^2*A1.*G1(x)+...
        (A1+h*(x-x1)).*G2(x)+...
        (A0-((l-x)/x1).^2.*A1).*G3(x)+...
        A0*H(x-l);
    x0 = xy_rot(2,:)-l/2;
    S1 = f_A((idx-0.5)-x0)/A0;
    S2 = f_A((idx+0.5)-x0)/A0;
    for kn = 1:N
        if idx(kn)-1 >= 1
        A(idx(kn)-1,kn) = S1(kn);
        end
        A(idx(kn)  ,kn) = S2(kn)-S1(kn);
        if idx(kn)+1 <= D
        A(idx(kn)+1,kn) = 1-S2(kn);
        end
        p(ixM(kn)) = pvec(idx(kn));
        ix(ixM(kn)) = true;
    end
    W((D*(kp-1)+(1:D)),:) = A; %#
end
% %
%Delete all-zero rows in W
fprintf('\nDelete all-zero rows in W...\r')
%ix = sum(W,2)==0;
try
    W(~ix,:) = [];
catch expr
    fprintf([expr.message '\nConverting to sparse...\r'])
    W = sparse(W);
    W(~ix,:) = [];
end
p(~ix) = [];
n_eq = size(W,1);
fprintf('\nFinish building A. \n%d equations in total.\r',n_eq);
figure('name','Sparsity pattern of matrix W');
spy(W)
t_e = toc;
t = t_e-t_s;
fprintf('\nTotal: %f seconds\n\r',t);
end


text:
function varExplained = rmVarExplained(view, roi, modelNum);
% Compute the proportion variance explained by a retinotopy model fit.
%
% varExplained = rmVarExplained([view], [coords or roi], [modelNum]);
%
%
% ras, 12/2006
if notDefined('view'),      view = getCurView;                          end
if notDefined('modelNum'),  modelNum = viewGet(view, 'rmModelNum');     end
if notDefined('coords'),    roi = rmGetCoords(view, modelNum);       end

roi = tc_roiStruct(view, roi);

verbose = prefsVerboseCheck;
if verbose, hwait = mrvWaitbar(0, 'Loading % Variance Explained...'); end

params = viewGet(view, 'rmParams');
[tSeries coords params] = rmLoadTSeries(view, params, roi);
if verbose, mrvWaitbar(.5, hwait); end

fit = rmPredictedTSeries(view, coords, modelNum);
if verbose,  mrvWaitbar(.7, hwait); end

residual = tSeries - fit;
if verbose,  mrvWaitbar(1, hwait); end

nVoxels = size(tSeries, 2);
for v = 1:nVoxels
    varExplained(v) = 1 - (sum(residual(:,v).^2) ./ sum(tSeries(:,v).^2));
end

return


text:
function F_BBHEwVEDg_gray = BBHEwVEDg_gray( F )

% clear all;
% close all;
% clc;
%x0=imread('eswar.jpg');
x0=F;
%x0=imread('11.bmp');
% x0=imread('12.bmp');
%y0=rgb2ycbcr(x0);
luma=x0;%(:,:,1);
[m,n]=size(luma);
% cb=y0(:,:,2);
% cr=y0(:,:,3);
k=mean(mean(luma));
r=round(k);
l=length(luma(luma<=r));
         listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
  sum=0;
 l1=length(luma(luma<=k3));%find all  values in luma falls below fist mean
 listindex1=find(luma<=k3);%all  values in luma falls below fist mean given in their locations
 k4=reshape(luma(listindex1),l1,1);%arrange all fist mean values in a vector
     xpdf=hist(k4,[0:k3]);%pdf from 0:r
      xpdf=xpdf/l1;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form 0 to r.
%       plot(xpdf);
%       xlabel('gray levels up to mean');
%       ylabel('pdf up to mean');
%       title('histogram for half an image up to mean');
      sk=xpdf*triu(ones(k3+1));
%       figure(2);
%       plot(sk);
%       xlabel('gray levels upto 1st mean');
%       ylabel('cdf upto mean');
%       title('cdf for half  of an image up to ist mean');
      alpha=0.6;
      for l2=0:k3
              list1=find(k4==l2);%find value in an vector i.e converted from matrix
              list(list1)=alpha*sk(l2+1)*(k3+1)+(1-alpha)*(k3+1);
              %list(list1)=sk(l2+1)*(k3+1);%map dont disturb to get bhe as
              %it is 13/3/2011
              ert(l2+1)=alpha*sk(l2+1)*(k3+1)+(1-alpha)*(k3+1);
      end
          p=zeros(m,n);
      p(listindex1)= list;
%      figure(3);
%      imshow(p);
%      xlabel('gray levels up to first mean');
%       ylabel('luma component equilized image up to first mean');
%      title('processed luma image up to first mean');
     k=mean(mean(luma));
r=round(k);
l=length(luma(luma<=r));
        listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
%   sum=0;
   b=k3;
%  l2=length(luma(luma>k3));%find all  values in luma falls below fist mean
  listindex2=find((luma>k3)&(luma<=r));%all  values in luma falls below fist mean given in their locations
  l2=length(listindex2);
  k5=reshape(luma(listindex2),l2,1);%arrange all fist mean values in a vector
      x2pdf=hist(k5,[k3+1:r]);%pdf from 0:r
       x2pdf=x2pdf/l2;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form 0 to r.
%        figure(4);
%        plot(x2pdf);
%        xlabel('gray levels 2nd mean');
%        ylabel('pdf of 2nd mean');
%        title('histogram for 2nd mean');
       sk2=x2pdf*triu(ones(r-k3));
%        figure(5);
%        plot(sk2);
%        xlabel('gray levels of 2nd mean');
%        ylabel('cdf upto mean');
%        title('cdf for half  of an image of 2nd mean');
       k2u=1;
       for l3=k3+1:r
               list2=find(k5==l3);%find value in an vector i.e converted from matrix
               list1(list2)=alpha*sk2(k2u)*(r-k3)+(1-alpha)*(k3+1);
               %list1(list2)=(k3+1)+(sk2(k2u))*(r-k3);%map dont disturb to
               %get BHE 13/3/2011
               ert(l3)=alpha*sk2(k2u)*(r-k3)+(1-alpha)*(k3+1);
               k2u=k2u+1;
       end
           p1=zeros(m,n);
       p1(listindex2)= list1;
%       figure(6);
%       imshow(p1);
%      xlabel('gray levels up to first 2nd mean');
%       ylabel('luma component equilized image 2nd mean');
%      title('processed luma image up to 2nd mean');
 %lupper30=length(luma(luma>r);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 lupper30=length(luma(luma>r));
%for i=0:r
    %if(luma(luma<=r))
        listindexupper3=find(luma>r);
   % end
%end
 k1upper=reshape(luma(listindexupper3),lupper30,1);
 mean3=mean(k1upper);
  r3=round(mean3);
 %length30=length((luma<=r3)&(luma>r));
  listindexupper30=find((luma>r)&(luma<=r3));
  length30=length(listindexupper30);
  k30upper=reshape(luma(listindexupper30),length30,1);

  %length30=length((luma<=r3)&(luma>r));
%  sum=0;
      xpdfupper30=hist(k30upper,[r+1:r3]);%pdf from r+1:r3
      xpdfupper30=xpdfupper30/length30;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form r+1 to 255.
%       figure(7);
%      plot(xpdfupper30);
%      xlabel('gray levels   3rd mean');
%      ylabel('pdf of 3rd mean');
%      title('histogram for upper half an image  3rd mean');
      skupper30=xpdfupper30*triu(ones(r3-r));
%       figure(8);
%       plot(skupper30);
%       xlabel('gray levels after mean');
%       ylabel('cdf after mean');
%       title('cdf for upper half  of an image after mean');
      k3u=1;
      for k3upper=(r+1):r3
              list1upper30=find(k30upper==k3upper);%find value in an vector i.e converted from matrix
              listnew(list1upper30)=(r+1)+skupper30(k3u)*(r3-r);%map
              ert(k3upper)=(r+1)+skupper30(k3u)*(r3-r);
         k3u=k3u+1;
     end

          p2=zeros(m,n);

     p2(listindexupper30)= listnew;
%     figure(9);
%      imshow(p2);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     listindexupper40=find(luma>r3);
     lupper40=length(listindexupper40);
   % end
%end
 k1upper40=reshape(luma(listindexupper40),lupper40,1);
 %sum=0;
 %for i=0:r
     xpdfupper40=hist(k1upper40,[r3+1:255]);%pdf from r+1:255
     xpdfupper40=xpdfupper40/lupper40;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form r+1 to 255.
%      figure(10);
%      plot(xpdfupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('pdf after 4mean');
%      title('histogram for upper half an image after 4mean');
     skupper40=xpdfupper40*triu(ones(255-r3));
%      figure(11);
%      plot(skupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('cdf after 4mean');
%      title('cdf for upper half  of an image after 4mean');
     k4u=1;
     for k4upper=(r3+1):255
         %if(xpdfupper(k2upper)>r)
             list1upper40=find(k1upper40==k4upper);%find value in an vector i.e converted from matrix
             %for k2u=1:58
             listnew4(list1upper40)=(r3+1+skupper40(k4u)*(255-r3));%map
         %end
         ert(k4upper)=(r3+1+skupper40(k4u)*(255-r3));
         k4u=k4u+1;
     end

%      for i=0:l-1
          p3=zeros(m,n);
%          if (p(listindex))
%              p(:)=list;
%          end
%      end

     p3(listindexupper40)= listnew4;
%     figure(12);
%     imshow(p3);
%     xlabel('gray levels after 4mean');
%      ylabel('luma component equilized image after 4mean');
%      title('processed luma image after 4mean');
     om=p+p1+p2+p3;
     F_BBHEwVEDg_gray=om;
     %ommmmm=p1+p;
%     figure(13);
%  imshow(uint8(om));
    % colormap('gray');
%       xlabel('gray level');
%      ylabel('combined lower and upper half luma component equilized image');
%      title('combined luma image');
%     figure(8);
    %image(om);
%     for j1=0:255
%         count=0;
%         for i1=0:m*n-1
%             if om(i1+1)==j1
%                 count=count+1;
%             end
%         end
%             prob(j1+1)=count/m*n;
%     end
%     figure(16);
%     plot(prob);
%
%    for j2=0:255
%         count1=0;
%         for i2=0:m*n-1
%             if luma(i2+1)==j2
%                 count1=count1+1;
%             end
%         end
%             prob2(j2+1)=count1/m*n;
%     end
%     figure(17);
%     plot(prob2);
%     xlabel('gray levels after mean');
%      ylabel('luma component equilized image after mean');
%      title('processed luma image after mean');
%     ommmmm=p1+p;
%     figure(7);
%      colormap('gray');
%       xlabel('gray level');
%      ylabel('combined lower and upper half luma component equilized image');
%      title('combined luma image');
%     image(ommmmm);
%      cat1=cat(3,om,cb,cr);
%     figure(14);
%     imshow(cat1);
% %     xlabel('gray level(ycbcr)');
% %      ylabel('combined lower and upper half luma,cromablue,croma red component equilized image');
% %      title('luma croma b and r color processed image');
%     catconversion=ycbcr2rgb(cat1);
%     figure(15);
%     imshow(catconversion);
%     xlabel('gray level(rgb)');
%      ylabel('combined lower and upper half  RGB component equilized image');
%      title('converted from ycbcr2rgb color(RGB) processed image');

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function s7_resave_final_check_for_errors()

stp = load('annot_full_Philippe.mat'); % train+test
phill_annot = stp.videos;
clear stp;

stp = load('annot_apt.mat'); % train+test
apt_annot = stp.videos;
clear stp;

philList = {};
for i=1:length(phill_annot)
    philList{i,1} = phill_annot(i).name; % store all the cor vid in a cellarray
end

aptList = {};

for i=1:length(apt_annot)
    aptList{i,1} = apt_annot(i).name; % store all the cor vid in a cellarray
end

images_dir = '/mnt/sun-alpha/datasets/UCF101/images';
save_path = 'saved/';

load('testnames.mat')
% load('merged_annot_remianing.mat')
load('checkpoint_improved_after_remaining.mat')
load('newAnnotats.mat')
vc = 0; tdata = cell(2,1); %v=1;
while v <= length(merged_annot)

    num_imgs = merged_annot(v).num_imgs;
    videoname = merged_annot(v).name;

    if isgood(v)
        [base_tubes,cc,ca] = verifyTubes(merged_annot(v).tubes,num_imgs,v,videoname);
        vc = vc +cc;
        if ca

            index = find(strcmp(philList,videoname));
            phil_tubes = phill_annot(index).tubes;

            index = find(strcmp(aptList,videoname));
            apt_tubes = apt_annot(index).tubes;
            base_tubes = merged_annot(v).tubes;
            if v == 1855
                base_tubes.sf = 1;
                base_tubes.ef = num_imgs;
                base_tubes.boxes = repmat(base_tubes.boxes,num_imgs,1);
                base_tubes.class = base_tubes.class;
            end
            checkboxes(base_tubes,phil_tubes,apt_tubes,num_imgs,images_dir,save_path,videoname)

            fid = fopen('temp.txt','w');

            fid = fopen('temp.txt','w');

            for t = 1 : length(base_tubes)
                fprintf(fid,'b %d %d %d\n', t, base_tubes(t).sf, base_tubes(t).ef);
            end

            for t = 1 : length(phil_tubes)
                fprintf(fid,'p %d %d %d\n', t, phil_tubes(t).sf, phil_tubes(t).ef);
            end

            for t = 1 : length(apt_tubes)
                fprintf(fid,'a %d %d %d\n', t, apt_tubes(t).sf, apt_tubes(t).ef);
            end

            fclose(fid);
            fid = fopen('temp.txt','r');

            tline = fgets(fid);
            newcell = cell(1);
            count = 1;
            newcell{count} = tline;
            while ischar(tline)
                tline = fgets(fid);
                if ischar(tline)
                    count = count + 1;
                    newcell{count} = tline;
                    disp(tline)
                end
            end
            fclose(fid);
            tdata{v} = newcell;
            class = phil_tubes(1).class;
            base_tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class);
            merged_annot(v).tubes = base_tubes;
            save('newAnnotats.mat','merged_annot','v');
        end

        merged_annot(v).tubes = base_tubes;

    end
    v=v+1;
end
save('newAnnotats.mat','merged_annot','v');

% fprintf('corrected %d tubes\n',vc);
% save('../testlist','testlist');
% annot = trainannot;
% save('../trainAnnot.mat','annot');
% annot = testannot;
% save('../testAnnot.mat','annot')
% annot = combinedannot;
% save('../finalAnnots.mat','annot')

function [newtubes,cc,ca] = verifyTubes(tubes, num_imgs,v,videoname)
cc = 0;tc = 0;ca=0;
newtubes = struct();
if isfield(tubes,'ef')>0
    for t = 1 : length(tubes)
        ef = tubes(t).ef;
        sf = tubes(t).sf;
        if ef>sf
            fdiff = abs(num_imgs-ef);
            if fdiff==1 || fdiff==2
                cc = cc+1;
                tubes(t).ef = num_imgs;
                for kk = 1:fdiff
                    tubes(t).boxes = [tubes(t).boxes;tubes(t).boxes(end,:)];
                end
            end
            tc = tc +1;
            newtubes(tc).ef = tubes(t).ef;
            newtubes(tc).sf = tubes(t).sf;
            newtubes(tc).boxes = tubes(t).boxes;
            newtubes(tc).class = tubes(t).class;
        else
            fprintf('%d %s\n',v,videoname);
            ca=1;
        end
    end
else
    fprintf('tube is empty %d %s\n',v,videoname);
end

function checkboxes(base_tubes,phil_tubes,apt_tubes,num_imgs,images_dir,save_path,videoName)

for i = 1 :1: num_imgs
    % fprintf('i=%d\n', i);
    philboxes = [];
    baseboxes = [];
    aptboxes = [];
    for t = 1 : length(phil_tubes)
        if phil_tubes(t).sf<= i && phil_tubes(t).ef >=i
            offset = uint16(phil_tubes(t).sf-1);
            philboxes = [philboxes;phil_tubes(t).boxes(i-offset,:),t];
        end
    end
    if isfield(base_tubes,'sf')
        for t = 1 : length(base_tubes)
            if i >= base_tubes(t).sf && i <= base_tubes(t).ef
                offset = uint16(base_tubes(t).sf-1);
                baseboxes = [baseboxes;base_tubes(t).boxes(uint8(i-offset),:),t];
            end
        end
    end

    for t = 1 : length(apt_tubes)
        if i >= apt_tubes(t).sf && i <= apt_tubes(t).ef
            offset = uint16(apt_tubes(t).sf-1);
            aptboxes = [aptboxes;apt_tubes(t).boxes(uint8(i-offset),:),t];
        end
    end

    image_name = sprintf('%s/%s/%05d.jpg',images_dir,videoName,i);
    %     if ~isempty(philboxes) || ~isempty(baseboxes)
    plotBoxes(philboxes,baseboxes,aptboxes,image_name, videoName, save_path, i,num_imgs);
    %     pause(0.000000001);
    %     end
end

function plotBoxes(philboxes,baseboxes,aptboxes,img_file, videoName, savepath, framenr,num_imgs)

if ~exist(savepath,'dir')
    mkdir(savepath);
end

savefile = sprintf('%s%05d.jpg',savepath,framenr);

%img_w=500;
%img_h=500;

hold off
im = imread(img_file);
imshow(im);
%set(gca,'Units','normalized','Position',[0 0 1 1]);  %# Modify axes size
%set(gcf,'Units','pixels','Position',[200 200 img_w img_h]);  %# Modify figure size
hold on;
legends = cell(1);
count =1 ;
steps = 1:15:440;

for  b = 1 : size(philboxes,1)
    rectangle('Position', philboxes(b,1:4), 'EdgeColor','b','LineWidth',5.0 );
    text(double(philboxes(b,1)+15), double(philboxes(b,2)-15), sprintf('%d',uint8(philboxes(b,5))),'Color','blue','FontSize',14);
    count = count + 1;
end

for  b = 1 : size(baseboxes,1)
    rectangle('Position', baseboxes(b,1:4), 'EdgeColor','r','LineWidth',3.0 );
    text(double(baseboxes(b,1)), double(baseboxes(b,2)-15), sprintf('%d',uint8(baseboxes(b,5))),'Color','red','FontSize',14);
    %     text(-50,steps(count),sprintf('base %d',b),'Color','red');
    count = count + 1;
end

for  b = 1 : size(aptboxes,1)
    rectangle('Position', aptboxes(b,1:4), 'EdgeColor','g','LineWidth',1.0 );
    text(double(aptboxes(b,1))+30, double(aptboxes(b,2)-15), sprintf('%d',uint8(aptboxes(b,5))),'Color','green','FontSize',14);
    %     text(-50,steps(count),sprintf('apt %d',b),'Color','green');
    count = count + 1;
end

% title([num2str(num_imgs),' ',num2str(framenr)])
% pause(0.000001)

frm_save = getframe(gcf); %# Capture the current window
imwrite(frm_save.cdata, savefile);
% imwrite(im,savefile);
%disp('ok');

function gtVidInd = getVidInd(video,videoName)
for i=1:length(video)
    vidid = video(i).name;
    if strcmp(vidid,videoName)
        gtVidInd = i;
        break;
    end
end

function tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class)
tubes = struct();
alltubes = {base_tubes,phil_tubes,apt_tubes};
for t = 1 : length(newcell)
    tline = newcell{t};
    stl = strsplit(tline,' ');

    if length(stl)<6
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,0);
    else
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,1);
    end
    tubes(t).sf = sf;
    tubes(t).ef = ef;
    tubes(t).class = class;
    tubes(t).boxes = boxes;
end

function [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,isone)
tube  = [];
annotType = {'b','p','a'};
ind = find(strcmp(annotType,stl{1}));

annot = alltubes{ind};

tid = str2num(stl{2});

sf = max(str2num(stl{3}),annot(tid).sf);

ef = min(str2num(stl{4}),annot(tid).ef);
st = sf-annot(tid).sf+1;
et = ef-sf+st; ct = 1;
boxes = annot(tid).boxes(st:et,:);
if isone
    ind = find(strcmp(annotType,stl{5}));
    annot = alltubes{ind};
    tid = str2num(stl{6});
    sfa = str2num(stl{7});
    if (sfa-ef)>1
        for kk = 1:(sfa-ef-1)
            boxes = [boxes;boxes(end,:)];
        end
    end
    ef = str2num(stl{8});
    st = sfa-annot(tid).sf+1;
    et = ef-sfa+st; ct = 1;
    boxes = [boxes;annot(tid).boxes(st:et,:)];
end

author: gurkirt


text:
function [ triangle_index, alpha, beta, gamma, edge, step_num ] = ...
  triangulation_search_delaunay ( node_num, node_xy, triangle_order, ...
  triangle_num, triangle_node, triangle_neighbor, p )

%*****************************************************************************80
%
%% TRIANGULATION_SEARCH_DELAUNAY searches a Delaunay triangulation for a point.
%
%  Purpose:
%
%    The algorithm "walks" from one triangle to its neighboring triangle,
%    and so on, until a triangle is found containing point P, or P is found
%    to be outside the convex hull.
%
%    The algorithm computes the barycentric coordinates of the point with
%    respect to the current triangle.  If all three quantities are positive,
%    the point is contained in the triangle.  If the I-th coordinate is
%    negative, then P lies on the far side of edge I, which is opposite
%    from vertex I.  This gives a hint as to where to search next.
%
%    For a Delaunay triangulation, the search is guaranteed to terminate.
%    For other triangulations, a cycle may occur.
%
%    Note the surprising fact that, even for a Delaunay triangulation of
%    a set of points, the nearest point to P need not be one of the
%    vertices of the triangle containing P.
%
%    The code can be called for triangulations of any order, but only
%    the first three nodes in each triangle are considered.  Thus, if
%    higher order triangles are used, and the extra nodes are intended
%    to give the triangle a polygonal shape, these will have no effect,
%    and the results obtained here might be misleading.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2012
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the vertices.
%
%    Input, integer TRIANGLE_ORDER, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles in the triangulation.
%
%    Input, integer TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM),
%    the nodes that make up each triangle.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle
%    neighbor list.
%
%    Input, real P(2), the coordinates of a point.
%
%    Output, integer TRIANGLE_INDEX, the index of the triangle where the
%    search ended.  If a cycle occurred, then TRIANGLE_INDEX = -1.
%
%    Output, real ALPHA, BETA, GAMMA, the barycentric
%    coordinates of the point relative to triangle TRIANGLE_INDEX.
%
%    Output, integer EDGE, indicates the position of the point P in
%    triangle TRIANGLE_INDEX:
%    0, the interior or boundary of the triangle;
%    -1, outside the convex hull of the triangulation, past edge 1;
%    -2, outside the convex hull of the triangulation, past edge 2;
%    -3, outside the convex hull of the triangulation, past edge 3.
%
%    Output, integer STEP_NUM, the number of steps.
%
  persistent triangle_index_save;

  dim_num = 2;

  step_num = - 1;
  edge = 0;

  if ( length ( triangle_index_save ) == 0 )
    triangle_index_save = -1;
  end

  if ( triangle_index_save < 1 | triangle_num < triangle_index_save )
    triangle_index = floor ( ( triangle_num + 1 ) / 2 );
  else
    triangle_index = triangle_index_save;
  end

  while ( 1 )

    step_num = step_num + 1;

    if ( triangle_num < step_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGULATION_SEARCH_DELAUNAY - Fatal error!\n' );
      fprintf ( 1, '  The algorithm seems to be cycling.\n' );
      triangle_index = -1;
      alpha = -1.0;
      beta = -1.0;
      gamma = -1.0;
      edge = -1;
      return
    end
%
%  Get the vertices of triangle TRIANGLE_INDEX.
%
    a = triangle_node(1,triangle_index);
    b = triangle_node(2,triangle_index);
    c = triangle_node(3,triangle_index);
%
%  Using vertex C as a base, compute the distances to vertices A and B,
%  and the point P.
%
    dxa = node_xy(1,a) - node_xy(1,c);
    dya = node_xy(2,a) - node_xy(2,c);

    dxb = node_xy(1,b) - node_xy(1,c);
    dyb = node_xy(2,b) - node_xy(2,c);

    dxp = p(1)         - node_xy(1,c);
    dyp = p(2)         - node_xy(2,c);

    det = dxa * dyb - dya * dxb;
%
%  Compute the barycentric coordinates of the point P with respect
%  to this triangle.
%
    alpha = ( dxp * dyb - dyp * dxb ) / det;
    beta =  ( dxa * dyp - dya * dxp ) / det;
    gamma = 1.0 - alpha - beta;
%
%  If the barycentric coordinates are all positive, then the point
%  is inside the triangle and we're done.
%
    if ( 0.0 <= alpha && 0.0 <= beta && 0.0 <= gamma )
      break
    end
%
%  At least one barycentric coordinate is negative.
%
%  If there is a negative barycentric coordinate for which there exists
%  an opposing triangle neighbor closer to the point, move to that triangle.
%
%  (Two coordinates could be negative, in which case we could go for the
%  most negative one, or the most negative one normalized by the actual
%  distance it represents).
%
    if ( alpha < 0.0 && 0 < triangle_neighbor(2,triangle_index) )
      triangle_index = triangle_neighbor(2,triangle_index);
      continue;
    elseif ( beta < 0.0 && 0 < triangle_neighbor(3,triangle_index) )
      triangle_index = triangle_neighbor(3,triangle_index);
      continue;
    elseif ( gamma < 0.0 && 0 < triangle_neighbor(1,triangle_index) )
      triangle_index = triangle_neighbor(1,triangle_index);
      continue;
    end
%
%  All negative barycentric coordinates correspond to vertices opposite
%  sides on the convex hull.
%
%  Note the edge and exit.
%
    if ( alpha < 0.0 )
      edge = -2;
      break
    elseif ( beta < 0.0 )
      edge = -3;
      break
    elseif ( gamma < 0.0 )
      edge = -1;
      break
    end

  end

  triangle_index_save = triangle_index;

  return
end


text:
function s = char(tpl)
%TEMPLATE Convert a template object in a one line description string
%  S = CHAR(TPL) is a class convertor from Template to a string, used
%  in online display.
%
%  See also DISPLAY

%  Copyright (C) 2003 Guillaume Flandin <Guillaume@artefact.tk>
%  $Revision: 1.0 $Date: 2003/05/05 22:19:51 $

s = ['Template Object: root ''',...
		tpl.root,''', ',...
		num2str(length(tpl.file)), ' files, ',...
		num2str(length(tpl.varkeys)), ' keys, ',...
		tpl.unknowns, ' unknowns.'];

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function Y = Rotation_Velocity_ToTireFixed(X,A_K_F,A_Rfl_F,A_Rfr_F,A_down_up)
%% Input parameters:
% A_K_F [---]    Rotation Matrix from RearAxle Fixed Coordinate System to VehicleFixed Coordinate System
% X     [---]    Original Input Vectors [x1 x2 x3 x4] [3x4]
% Output parameters:
% Y     [---]    Transformed Output Vectors [y1 y2 y3 y4] [3x4]

%% Output of transformed vectors
% Quarter Vehicle Model along Vehicle Fixed z-Axis: Rotation from VehicleFixed Coordinate System into TireFixed Coordinate System
% Y=[(A_down_up*A_Rfl_F*transpose(A_K_F)*X(:,1)) (A_down_up*A_Rfr_F*transpose(A_K_F)*X(:,2)) (A_down_up*transpose(A_K_F)*X(:,3)) (A_down_up*transpose(A_K_F)*X(:,4))];

% Quarter Vehicle Model along Rear Axle Fixed z-Axis: Rotation from RearAxleFixed Coordinate System into TireFixed Coordinate System
Y=[(A_down_up*A_Rfl_F*X(:,1)) (A_down_up*A_Rfr_F*X(:,2)) (A_down_up*X(:,3)) (A_down_up*X(:,4))];

% Quarter Vehicle Model according MATLAB-Documentation
% Y=[(A_down_up*A_Rfl_F*X(:,1)) (A_down_up*A_Rfr_F*X(:,2)) (A_down_up*X(:,3)) (A_down_up*X(:,4))];

end

author: TUMFTM


text:
function [w,b,Yt] = ridgeregress(X,Y,gam, Xt)
% Linear ridge regression
%
% >> [w, b]     = ridgeregress(X, Y, gam)
% >> [w, b, Yt] = ridgeregress(X, Y, gam, Xt)
%
% Ordinary Least squares with a regularization parameter (gam).
%
% Full syntax
%
% >> [w, b, Yt] = ridgeregress(X, Y, gam, Xt)
%
% Outputs
%   w     : d x 1 vector with the regression coefficients
%   b     : bias term
%   Yt(*) : Nt x 1 vector with predicted outputs of test data
% Inputs
%   X     : N x d matrix with the inputs of the training data
%   Y     : N x 1 vector with the outputs of the training data
%   gam   : Regularization parameter
%   Xt(*) : Nt x d matrix with the inputs of the test data
%
% See also:
% bay_rr,bay_lssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

if size(X,1)~=size(Y,1),
  error('X and Y need to have the same number of data points');
end
if size(Y,2)~=1,
  error('Only handling one-dimensional output');
end
if nargin==4 & size(Xt,2)~=size(X,2),
  error('Training input and test inputs need to have the same dimension');
end

[nD,nx] = size(X);
if nx>nD, warning('dim datapoints larger than number of datapoints...');end

Xe = [X ones(nD,1)];
%H = [ Xe'*Xe + gam^-1*[eye(nx) zeros(nx,1); zeros(1,nx+1)]];
H = Xe'*Xe + inv(gam).*eye(nx+1);

sol = pinv(H)*Xe'*Y;
w = sol(1:end-1);
b = sol(end);

if nargin<4, return; end
Yt = Xt*w+b;

author: peterhcharlton

repo: RRest


text:
%INITWIDEANGLEPROJMAP  Initializes maps for cv.remap for wide-angle
%
%     [map1, map2] = cv.initWideAngleProjMap(cameraMatrix, distCoeffs, imageSize, destImageWidth)
%     [map1, map2, scale] = cv.initWideAngleProjMap(cameraMatrix, distCoeffs, imageSize, destImageWidth)
%     [...] = cv.initWideAngleProjMap(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __cameraMatrix__ Input camera matrix `A = [f_x 0 c_x; 0 f_y c_y; 0 0 1]`
% * __distCoeffs__ Input vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12 or 14
%   elements. If the vector is empty, the zero distortion coefficients are
%   assumed.
% * __imageSize__ image size `[w,h]`.
% * __destImageWidth__
%
% ## Output
% * __map1__ The first output map. See `M1Type`.
% * __map2__ The second output map. See `M1Type`.
% * __scale__
%
% ## Options
% * __M1Type__ Type of the first output map, default `single2`. See
%   cv.convertMaps. Accepted types are:
%   * __int16__ first output map is a MxNx2 `int16` array, second output map
%     is MxNx1 `uint16` (fixed-point representation).
%   * __single1__ first output map is a MxNx1 `single` matrix, second output
%     map is MxNx1 `single` (separate floating-point representation).
%   * __single2__ first output map is a MxNx2 `single` matrix, second output
%     map is empty (combined floating-point representation).
% * __ProjType__ projection type, default 'EqRect'. One of:
%   * __Ortho__
%   * __EqRect__
% * __Alpha__ default 0
%
% See also: cv.initUndistortRectifyMap, cv.remap, cv.convertMaps
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function M = ply_read(filename)
% Read PLY-formatted data from disk
% FORMAT M = ply_read(filename)
%
% filename - PLY-formatted file name
% M        - data structure
%__________________________________________________________________________
%
% Stanford Triangle Format Specification:
% https://en.wikipedia.org/wiki/PLY_%28file_format%29
%__________________________________________________________________________
% Copyright (C) 2017-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: ply_read.m 7676 2019-10-22 10:19:29Z guillaume $

fid = fopen(filename,'rt');
if fid == -1
    error('Cannot open %s.',filename);
end

M = struct('vertices',[],'faces',[]);

%-Read header
isvertex = false;
nval = 0;
while true
    l = fgetl(fid);
    if strcmp(l,'end_header'), break; end
    [l,r] = strtok(l);
    if strcmp(l,'element')
        [l,r] = strtok(r);
        switch l
            case 'vertex'
                nv = str2double(r);
                isvertex = true;
            case 'face'
                nf = str2double(r);
                isvertex = false;
        end
    elseif strcmp(l,'property')
        if isvertex
            nval = nval + 1;
        end
    end
end

%-Read data
M.vertices = fscanf(fid,'%f',[nval nv])';
M.vertices = M.vertices(:,1:3);
M.faces    = fscanf(fid,'%d',[4 nf])';
M.faces    = M.faces(:,2:4) + 1;

fclose(fid);

author: spm


text:
function [X_norm, mu, sigma] = featureNormalize(X)
% FEATURENORMALIZE(X) returns a normalized version of X where
% the mean value of each feature is 0 and the standard deviation
% is 1. This is often a good preprocessing step to do when
% working with learning algorithms.

mu = mean(X);
sigma = std(X);

t = ones(length(X), 1);
X_norm = (X - (t * mu)) ./ (t * sigma); % Vectorized

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function p = presolve(p)

model = p.model;

p.model.Q = p.model.Q*0;
p.model.F_struc(1:p.dimin(1),:)=[];
p.model.K.f = p.model.K.f-p.dimin(1);
redundant = sparse(zeros(p.model.K.l,1));
for i = 1:p.model.K.l
    b = p.model.F_struc(p.model.K.f+i,1);
    p.model.F_struc(p.model.K.f+i,1) = b+1;
    p.model.c = p.model.F_struc(p.model.K.f+i,2:end)';

     eval(['output = ' p.model.solver.call '(p.model);']);
     if output.problem == 0
         if p.model.c'*output.Primal > -b
             redundant(i) = 1
         end
     end

     p.model.F_struc(p.model.K.f+i,1) = b;
end
p.model.c = c;
p.model.Q = c;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [sm,pro]=oper(A,B)

% This function computes the
% sum and the product of 2 matrices

sm=A+B;

pro=A*B;

% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni