text:
function [SGP4Elements,TTEpoch1,TTEpoch2,otherInfo,checksumIsBad]=TLE2SGP4OrbEls(TLELine1,TLELine2)
%%TLE2SGP4ORBELS Convert a North American Aerospace Defense Command (NORAD)
%                two-line element (TLE) satellite ephemerides into a format
%                that can be used with the propagateOrbitSGP4 function for
%                determining the position and velocity of the satellite at
%                different times using the Simplified General
%                Perturbations 4 (SGP4) model/ the Simplified Deep Space
%                Perturbations 4 (SDP4) model for distant satellites.
%                TLEs are how the Air Force Space Command (AFSPC) has
%                publicly released information on satellite orbits for
%                years.
%
%INPUTS: TLELine1, TLELine2 The character strings for each of the two lines
%                           that make up the TLE ephemeride, including all
%                           spaces as the positions of the characters are
%                           important. With checksums, these should both be
%                           69 characters long, without checksums, 68
%                           characters long.
%
%OUTPUTS: SGP4Elements The 7X1 vector of orbital elements that can be given
%                      to the propagateOrbitSGP4 function along with the
%                      epoch time to predict the location of the satellite
%                      at various times. They are similar, but not
%                      identical to, Keplerian orbital elements of the same
%                      name. If either of the TLE strings is too short,
%                      this is an empty matrix. The elements of
%                      SGP4Elements are
%                      1) eccentricity
%                      2) inclination  (radians)
%                      3) argument of perigee (radians)
%                      4) right ascension of the ascending node (radians)
%                      5) mean anomaly (radians)
%                      6) mean motion (radians per second [TT])
%                      7) BSTAR drag term. This pseudo ballistic
%                        coefficient has units of inverse Earth radii.
%                        Normally, a ballistic coefficient BC is mass per
%                        unit area divided by a drag coefficient. The BSTAR
%                        drag term is a scaled version of the normal
%                        ballistic coefficient. That is,
%                        BC=Re*rho_0/(2*BSTAR) where Re is the radius of
%                        the Earth in meters, and rho=2.461e-5 kg/m^2.
%   TTEpoch1, TTEpoch2 The epoch time of the SGP4Elements given in
%                      terrestrial time (TT) as a two-part Julian date.
%                      This will be incorrect from 2057 onward, since the
%                      TLE data only provides the final two digits of the
%                      year starting in 1957 and the century must be
%                      guessed. If either of the TLE strings is too short,
%                      these are empty matrices.
%            otherInfo A structure holding the other information present in
%                      the TLE. If either of the TLE strings is too short,
%                      this is an empty matrix. The fields of this
%                      structure are
%                      line1Number The number starting the first TLE line.
%                                  This should be 1.
%                      line2Number The number starting the second TLE line.
%                                  This should be 2.
%                      satelliteNumber1 The NORAD catalog number of the
%                                  satellite as given by the first line of
%                                  the TLE.
%                      satelliteNumber2 The NORAD catalog number of the
%                                  satellite as given by the second line of
%                                  the TLE. This should be the same as
%                                  satelliteNumber1
%                      classification A character indicating the
%                                  classification level of the satellite.
%                                  'U'=unclassified.
%                      IDLaunchYear The last two digits of the launch year
%                                  of the satellite.
%                      IDLaunchNumber The number of the launch in in the
%                                  launch year that put up the satellite.
%                      IDPieceOfLaunch A string indicating which piece of
%                                  the launch the satellite is.
%                      firstDerivMeanMotion The first derivative of the
%                                  mean motion in radians per second
%                                  squared. This is useful with the older
%                                  orbital propagators.
%                      secondDerivMeanMotion The second derivative of the
%                                  mean motion in radians per second cubed.
%                      ephemerisType Originally, this was a number from
%                                  1-5 indicating the type of orbital
%                                  propagator used to generate the
%                                  ephemeris with 1=SGP, 2=SGP4, 3=SDP4,
%                                  4=SGP8, and 5=SDP8. However, it is
%                                  always set to zero in TLEs that are
%                                  publicly released and is thus
%                                  meaningless.
%                      elementNumber A number indicating which TLE set this
%                                  TLE is from. The number is supposed to
%                                  be incremented everytime a new set comes
%                                  out. However, in practice, it is not.
%                      revolutionNumberAtEpoch The number of orbits the
%                                  satellite has made since its launch.
%        checksumIsBad A 2X1 boolean vector indicating whether the checksum
%                      value for each of the strings is bad.
%                      checksumIsBad(1) is true if the checksum of TLELine1
%                      is invalid and checksumIsBad(2) is true if the
%                      checksum of TLELine2 is invalid. If a checksums is
%                      omitted then the value will be true. If iether of
%                      the strings is too short, then this will be an empty
%                      matrix.
%
%Note that parts of the input strings that cannot be read will be filled
%with NaN values.
%
%Information of the format of the TLE sets can be found at [1]. The
%documentation discussing how the element sets are used is in [2] and [3].
%The decimal point parsing part of this function is loosely based on David
%Vallado's public-domain code for the SGP4 propagator that was downloaded
%from http://www.centerforspace.com/downloads/
%
%A sample TLE that accompanies Vallado's code is:
%TLELine1='1 11801U          80230.29629788  .01431103  00000-0  14311-1      13'
%TLELine2='2 11801  46.7916 230.4354 7318036  47.4722  10.4117  2.28537848    13'
%
%REFERENCES:
%[1] T. C. for Space Standards and Innovation. (2012, 27 Sep.) NORAD two-
%    line element sets. [Online].
%    Available: http://www.celestrak.com/NORAD/elements/
%[2] F. R. Hoots and R. L. Roehrich, "Spacetrack report no. 3: Models for
%    propagation of NORAD element sets," Department of Defense, Tech. Rep.,
%    31 Dec. 1988. [Online].
%    Available: http://www.amsat.org/amsat/ftp/docs/spacetrk.pdf
%[3] D. A. Vallado, P. Crawford, R. Hujsak, and T. S. Kelso, "Revisiting
%    spacetrack report # 3: Rev 2," in Proceedings of the AIAA/AAS
%    Astrodynamics Specialist Conference and Exhibit, Keystone, CO, 21-24
%    Aug. 2006. [Online].
%    Available: http://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753.pdf
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.*/
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Multiplication coefficient to convert from degrees to radians.
deg2Rad=pi/180;

%If either of the strings is too short.
if(length(TLELine1)<68||length(TLELine2)<68)
    SGP4Elements=[];
    TTEpoch1=[];
    TTEpoch2=[];
    otherInfo=[];
    checksumIsBad=[];
    return;
end

%If the strings have no checksums at the end, then just mark the checksums
%as bad.
if(length(TLELine1)<69)
    checksumIsBad(1)=true;
else
    checksumIsBad(1)=~validateTLEChecksum(TLELine1);
end

if(length(TLELine2)<69)
    checksumIsBad(2)=true;
else
    checksumIsBad(2)=~validateTLEChecksum(TLELine2);
end

%Set the implied decimal points to help deal with bad values when doing a
%formatted read.
for j=11:16
    if(TLELine1(j)==' ')
        TLELine1(j)='_';
    end
end

if(TLELine1(45)~=' ')
    TLELine1(44)=TLELine1(45);
end
TLELine1(45)='.';

if(TLELine1(8)==' ')
    TLELine1(8)='U';
end

if(TLELine1(10)==' ')
    TLELine1(10)='.';
end

for j=46:50
    if (TLELine1(j)==' ')
        TLELine1(j)='0';
    end
end

if(TLELine1(52)==' ')
    TLELine1(52)='0';
end

if(TLELine1(54)~=' ')
    TLELine1(53)=TLELine1(54);
end
TLELine1(54)='.';
TLELine2(26)='.';

for j=27:33
    if(TLELine2(j)==' ')
        TLELine2(j)='0';
    end
end

if(TLELine1(63)==' ')
    TLELine1(63)='0';
end

if((length(TLELine1)<68)||(TLELine1(68)==' '))
    TLELine1(68)='0';
end

%Parse the stuff that most folks do not care about.
otherInfo.line1Number = str2double(TLELine1(1));
otherInfo.line2Number = str2double(TLELine2(1));
otherInfo.satelliteNumber1 = str2double(TLELine1(3:7));
otherInfo.satelliteNumber2 = str2double(TLELine2(3:7));

otherInfo.classification = TLELine1(8);
otherInfo.IDLaunchYear = str2double(TLELine1(10:11));
otherInfo.IDLaunchNumber= str2double(TLELine1(12:14));
otherInfo.IDPieceOfLaunch= TLELine1(15:17);

%Convert revolutions per day squared into radians per second squared,
%with 86400 seconds per TT Julian day.
otherInfo.firstDerivMeanMotion=str2double(TLELine1(34:43))*(2*pi/86400^2);
%Convert revolutions per day cubed into radians per second cubed,
%with 86400 seconds per TT Julian day.
otherInfo.secondDerivMeanMotion=str2double(TLELine1(44:50))*10.0^str2double(TLELine1(51:52))*(2*pi/86400^3);
otherInfo.ephemerisType = str2double(TLELine1(63));
otherInfo.elementNumber = str2double(TLELine1(65:68));
otherInfo.revolutionNumberAtEpoch=str2double(TLELine2(64:68));

%Parse the stuff that most folks care about
inclination = str2double(TLELine2(8:16))*deg2Rad;
RAAscendingNode = str2double(TLELine2(17:25))*deg2Rad;
eccentricity = str2double(TLELine2(26:33));
argumentOfPerigee = str2double(TLELine2(34:42))*deg2Rad;
meanAnomaly = str2double(TLELine2(43:51))*deg2Rad;
%Convert from revolutions per day to radians per second, with 86400
%seconds per TT Julian day.
meanMotion=str2double(TLELine2(52:63))*(2*pi/86400);
BSTARDrag= str2double(TLELine1(53:59))*10.0^str2double(TLELine1(60:61));

%Get the time of the orbital elements.
epochYear = str2double(TLELine1(19:20));
epochDays = str2double(TLELine1(21:32));

%Assume that the years only run from 1957->2057 and find the actual
%epoch year.
if (epochYear < 57)
    year=epochYear+2000;
else
    year=epochYear+1900;
end
%Convert into a two-part Julian date in TT.
[month,dayOfMonth]=dayOfYear2MonthDay(year,fix(epochDays));
dayFrac=epochDays-fix(epochDays);
[hour,minute,second]=fracDayOfMonth2HourMinSec(year,month,dayOfMonth,dayFrac);
[TTEpoch1,TTEpoch2]=Cal2TT(year,month,dayOfMonth,hour,minute,second);

SGP4Elements=zeros(7,1);
SGP4Elements(1)=eccentricity;
SGP4Elements(2)=inclination;
SGP4Elements(3)=argumentOfPerigee;
SGP4Elements(4)=RAAscendingNode;
SGP4Elements(5)=meanAnomaly;
SGP4Elements(6)=meanMotion;
SGP4Elements(7)=BSTARDrag;
end

function isValid=validateTLEChecksum(theTLELine)
%The validation comes from the explanation of how the checksum is formed
%from http://www.celestrak.com/NORAD/elements/ It is the last digit of the
%sum of all of the digits in the string plus an extra 1 for each '-'.
    theSum=0;
    for curCharIdx=1:68
        curChar=theTLELine(curCharIdx);
        if(curChar=='-')
            theSum=theSum+1;
        elseif(curChar>='0'&&curChar<='9')
            theSum=theSum+str2double(curChar);
        end
    end

    %Get the last digit of the sum.
    sumStr=num2str(theSum);
    isValid=(sumStr(end)==theTLELine(69));
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%SF_TET_P5 Fifth order Lagrange shape functions for tetrahedra (P5).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_TET_P5( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming fifth order P5 Lagrange shape functions on 3D tetrahedral elements
%   with values defined in the nodes, edges, faces, and cell center. XI are Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 3              Number of space dimensions
%       n_vert      scalar: 4              Number of vertices per cell
%       i_dof       scalar: 1-56           Local basis function to evaluate
%       xi          array [4,1]            Local coordinates of evaluation point
%       aInvJac     [n,12]                 Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [4,n_ldof]             Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_TET_P1

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [bOK, sFilePath] = ex_PutFile(vFilterSpec, sDialogTitle, nXPos, nYPos)
% function [bOK, sFilePath] = ex_PutFile(vFilterSpec, sDialogTitle, nXPos, nYPos)
% -------------------------------------------------------------------------------
% Wrapper for the uiputfile function of Matlab. Simplifies the output parameter analysis
%
% Input parameter:
%   vFilterSpec   Filter rules for selection, e. g. '*.mat'
%   sDialogTitle  String to appear in the titlebar
%   nXPos         X-position of dialog
%   nYPos         Y-position of dialog
%
% Output parameters:
%   bOK         Result of user interaction: 1: OK pressed, 0: Cancel pressed
%   sFilePath   Complete filepath of selected file
%
% Danijel Schorlemmer
% December 18, 2001

% Init return values
bOK = 0;
sFilePath = '';
% Invoke Matlab's original function
if ~exist('vFilterSpec')
  vFilterSpec = '';
end
if ~exist('sDialogTitle')
  sDialogTitle = '';
end
if (exist('nXPos') & exist('nYPos'))
  [sFile, sPath] = uiputfile(vFilterSpec, sDialogTitle, nXPos, nYPos);
else
  [sFile, sPath] = uiputfile(vFilterSpec, sDialogTitle);
end
% Cancel pressed?
if isequal(sFile, 0)  ||  isequal(sPath, 0)
  return;
end
% Everything ok?
sFilePath = [sPath sFile];
if length(sFilePath) > 1
  % We got everything
  bOK = 1;
end

author: CelsoReyes


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS
%
%   Authors: Arturo Gil
%   email:  date:   27/11/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 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/>.
function robot = parameters()

robot.name= 'SERIAL2_EXAMPLE';

robot.DH.theta= '[q(1) q(2)-pi/2 q(3) q(4) q(5) q(6)+pi]';
robot.DH.d='[0.630 0 0 1.8 0 0.3]';
robot.DH.a='[0.50 2.34 0.5 0 0 0]';
robot.DH.alpha= '[pi/2 0 pi/2 0 pi/2 0]';

robot.J=[];

robot.inversekinematic_fn = 'inversekinematic(robot, T)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-170) deg2rad(170); %Axis 1, minimum, maximum
                deg2rad(-65) deg2rad(85); %Axis 2, minimum, maximum
                deg2rad(-180) deg2rad(70); %Axis 3
                deg2rad(-300) deg2rad(300); %Axis 4:
                deg2rad(-120) deg2rad(120); %Axis 5
                deg2rad(-360) deg2rad(360)]; %Axis 6:

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(100); %Axis 1, rad/s
                deg2rad(90); %Axis 2, rad/s
                deg2rad(90); %Axis 3, rad/s
                deg2rad(170); %Axis 4, rad/s
                deg2rad(120); %Axis 5, rad/s
                deg2rad(190)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
            % end effectors maximum velocity
robot.linear_velmax = 1.0; %m/s, unavailable from datasheet

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-2 2 -2 2 0 2.5];
%read graphics files
robot = read_graphics(robot);

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text:
function [tout, info] = transfdiff(opt, tp, tm)
% TRANSFDIFF Transform diffusion algorithm for sequence of images.
%
% This function is part of a larger algorithm that solves the following
% registration problem:
%
% We have a sequence of images I=1,2,...,N. We want to register each image
% to its adjacent neighbours to form a volume. But instead of registering
% I=2 to I=1, and then I=3 to the result, and then I=4 to the result, etc.,
% it can be shown that all images can be aligned in parallel using an
% iterative process that we call "transform diffusion registration".
%
% TRANSFDIFF solves the "transform diffusion" step in the algorithm:
%
% Let's assume that somebody has already registered each image I to its
% neighbours I-1 and I+1. TRANSFDIFF takes that set of registration
% parameters and "diffuses" them to produce a transform for each image.
% Applying these transforms to the images produces the globally aligned
% sequence.
%
%
% TRANSFDIFF applies a diffusion process to those neighbour transforms. The
% output is the N transforms that applied to the N images best align them.
%
% [TOUT, INFO] = TRANSFDIFF(OPT, ...)
%
%   OPT is a string with the name of the transform to apply, or a struct
%   with the transform name and algorithm parameters. Options common to all
%   transformations:
%
%     'MaxIter': (def 50) Stopping criterion. The algorithm stops after
%                MaxIter diffusion iterations.
%
%     'Alpha':   (def 0.45) Diffusion coefficient. The smaller the value of
%                alpha, more iterations are needed for the same result. But
%                alpha >0.45 may not smooth high frequencies very well, and
%                values >=0.5 cause oscillations.
%
%   TOUT is an array with N output transforms for the N slices.
%
%   INFO is a struct with information about the algorithm internals:
%
%     'NumIter': Number of diffusion iterations until stop condition.
%
% -------------------------------------------------------------------------
%
% * OPT='TranslationTransform'
% * OPT.Transform='TranslationTransform'
%
% [TOUT, INFO] = TRANSFDIFF(OPT, TP)
%
%   TP is a matrix with N-1 rows. TP(I, :) is the translation from slice I
%   to slice I+1.
%
%   This case assumes that the transforms are symmetric, i.e. the
%   translation from slice I+1 to I is -TP(I, :).
%
%   Specific OPT options:
%
%     'Epsilon': (def 0.0) Stopping criterion. It stops when all
%                translation components are small:
%                for all t, abs(t)<=Epsilon.
%
% -------------------------------------------------------------------------
%
% * OPT='AffineTransform'
% * OPT.Transform='AffineTransform'
%
% * OPT='EulerTransform'
% * OPT.Transform='EulerTransform'
%
% [TOUT, INFO] = TRANSFDIFF(OPT, TP)
%
%   TP is an array of affine transforms in homogeneous coordinates. TP(I,:)
%   is the affine transform from slice I to slice I+1. Each transform has
%   the format
%
%   TP(:, :, I) = [A 0]
%                 [d 1]
%
%   where A is a matrix and d is a translation vector, and the affine
%   transform of a row vector X is defined as Y=X*TP(:,:,I).
%
%   This case assumes that the transforms are symmetric, i.e. the
%   transformation from slice I+1 to I is inv(TP(I, :)).
%
%   Specific OPT options:
%
%     'Epsilon': (def 0.0) Stopping criterion. It stops when the norm of
%                the difference between the transformation and the identity
%                matrix is small:
%                for all t, norm(t-eye(3))<=Epsilon.
%
% -------------------------------------------------------------------------
%
% * OPT.Transform='BsplineTransform'
%
% [TOUT, INFO] = TRANSFDIFF(OPT, TP, TM)
%
%   TP is a matrix with N-1 rows. Each row contains the B-spline
%   coefficients with the transformation from slice I to slice I+1. The
%   format of the B-spline coefficient vector is
%   [p0x, p1x, ..., pNx, p0y, p1y, ..., pNy].
%
%   TM is similar to TP, but TM(I, :) is the transformation from slice I+1
%   to I.
%
%   Specific OPT options:
%
%     'Epsilon': (def 0.0) Stopping criterion. It stops when all
%                coefficient displacement components are small:
%                for all t, abs(t)<=Epsilon.
%
%     'tbsp': (def []) For Choi and Lee (2000) injectivity criteria. These
%             are sufficient criteria that guarantee the B-spline doesn't
%             fold over. Transform struct in elastix format. The
%             tbsp.TransformParameters are ignored, only the grid
%             information is used. If not provided or empty, the conditions
%             are not used as a stopping criterion, and the resulting
%             B-splines may fold over.
%
%
% -------------------------------------------------------------------------
%
% See also: transdiffreg.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015-2016 University of Oxford
% Version: 0.4.4
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
switch (opt.Transform)

    case {'TranslationTransform', 'EulerTransform', 'AffineTransform'}

        narginchk(2, 2);

    case {'BSplineTransform'}

        narginchk(2, 3);

    otherwise

        error(['Transform ' opt.Transform ' not implemented'])

end
nargoutchk(0, 2);

% defaults
if (isempty(opt))
    error('OPT cannot be empty');
elseif (ischar(opt))
    aux = opt;
    clear opt
    opt.Transform = aux;
end
if (~isfield(opt, 'Epsilon'))
    opt.Epsilon = 0.0;
end
if (~isfield(opt, 'MaxIter'))
    opt.MaxIter = 50;
end
if (~isfield(opt, 'Alpha'))
    opt.Alpha = 0.45;
end
if (strcmp(opt.Transform, 'BSplineTransform'))
    if (~isfield(opt, 'tbsp'))
        opt.tbsp = [];
    end
end

if (opt.Alpha < 0 || opt.Alpha > 0.5)
    warning('alpha must be in the interval [0.0, 0.5]. alpha=0.5 can produce oscillations. We recommend alpha<=0.45')
end

% convert to array format if transforms are provided in elastix format
[tp, tp0] = elastix2mat(tp);
if (nargin > 2)
    [tm, tm0] = elastix2mat(tm);
end

% preprocessing of inputs
switch (opt.Transform)

    case {'TranslationTransform', 'BSplineTransform'}

        % compute "tm" from "tp"
        if (nargin < 3 || isempty(tm))
            tm = -tp;
        end

        % add dummy transforms at the beginning of "tm" and at the end of
        % "tp". This makes it easier to operate, because then:
        % * tp(I) and tm(I) are the neighbours of I for intermediate slices
        % * extreme slices can be dealt with in the same way as
        %   intermediate slices
        tp = tp([1:end end], :);
        tm = tm([1 1:end], :);
        tp(end, :) = tm(end, :);
        tm(1, :) = tp(1, :);

    case {'EulerTransform', 'AffineTransform'}

        % compute "tm" from "tp"
        if (nargin < 3 || isempty(tm))
            for I = 1:size(tp, 3)
                tm(:, :, I) = inv(tp(:, :, I));
            end
        end

        % add dummy transforms at the beginning of "tm" and at the end of
        % "tp". This makes it easier to operate (see note above)
        tp = tp(:, :, [1:end end]);
        tm = tm(:, :, [1 1:end]);
        tp(:, :, end) = tm(:, :, end);
        tm(:, :, 1) = tp(:, :, 1);

end

% check inputs
if any(size(tp) ~= size(tm))
    error('TP and TM must have the same size')
end

% if no transforms provided, then we don't need to run the algorithm
if (isempty(tp))
    tout = [];
    return
end

% apply registration diffusion
switch (opt.Transform)

    case {'TranslationTransform', 'BSplineTransform'}

        [tout, info] = translation_diffusion(opt, tp, tm);

    case {'EulerTransform', 'AffineTransform'}

        [tout, info] = affine_diffusion(opt, tp, tm);

    otherwise

        error('Transform not implemented')

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOCAL FUNCTIONS

%% convert transform from elastix format to array, if necessary
function [tout, t0] = elastix2mat(t)

if (isstruct(t))

    % save copy of input so that we can convert the output back to elastix
    % format
    t0 = t;

    switch (t(1).Transform)

        case {'TranslationTransform', 'BSplineTransform'}

            tout = cat(1, t(:).TransformParameters);

        % TODO: AffineTransform

        otherwise

            error('Conversion from this elastix transform to Matlab array not implemented')
    end

else % already in array format

    % no need to save a copy of the input
    t0 = [];

    % output is just the input
    tout = t;

end

end

%% translation_diffusion: Apply diffusion to translation transforms
function [tout, info] = translation_diffusion(opt, tp, tm)

% init accumulated transform to apply to each slice
tout = zeros(size(tp));

I = 0;
while (1)

    % iteration number
    I = I + 1;

    % transform to apply to each slice in this iteration (transform update)
    t = opt.Alpha * (tp + tm);

    % stopping criterion: B-spline injectivity
    if (strcmp(opt.Transform, 'BSplineTransform'))

        % find control points that don't guarantee injectivity and
        % constrain them so that they do
        [info.NConstrained(I), aux] = ...
            check_bspline_injectivity_choi2000(tout + t, opt.tbsp);

        % apply the control point modifications to the B-spline update
        t = aux - tout;

    end

    % compose with the previous transforms to obtain a total transform
    tout = tout + t;

    % update neighbour transforms ...

    % ... of internal images
    tp(1:end-1, :) = t(2:end, :) + tp(1:end-1, :) - t(1:end-1, :);
    tm(2:end, :) = t(1:end-1, :) + tm(2:end, :) - t(2:end, :);

    %... of extreme images
    tp(end, :) = tm(end, :);
    tm(1, :) = tp(1, :);

    % stopping criterion: absolute value of each translation component
    tabs = abs(t(:));
    info.MaxAbs(I) = max(tabs);
    info.MeanAbs(I) = mean(tabs);
    info.MedAbs(I) = median(tabs);
    if (info.MaxAbs(I) <= opt.Epsilon)
        break
    end

    % stopping criterion: maximum number of iterations
    if (I == opt.MaxIter)
        break
    end

end

info.NumIter = I;

% DEBUG:
% isInjective = check_bspline_injectivity_global(tout + t, opt.tbsp);
% nnz(~isInjective)

end

%% affine_diffusion: Apply diffusion to affine transforms
function [tout, info] = affine_diffusion(opt, tp, tm)

% init accumulated transform to apply to each slice
tout = eye(size(tp(:, :, 1)));
tout = repmat(tout, 1, 1, size(tp, 3));

% allocate memory for transform update
t = zeros(size(tp));

I = 0;
while (1)

    % iteration number
    I = I + 1;

    % compute the transform to apply to each slice in this iteration
    for J = 1:size(tp, 3)

        % transform to apply to each slice in this iteration (transform
        % update)
        t(:, :, J) = real(expm(opt.Alpha ...
            * (logm(tp(:, :, J)) + logm(tm(:, :, J)))));

        % compose with the previous transforms to obtain a total transform
        tout(:, :, J) = tout(:, :, J) * t(:, :, J);

    end

    % update neighbour transforms of internal images
    for J = 2:size(tm, 3)
        tm(:, :, J) = t(:, :, J) \ tm(:, :, J) * t(:, :, J-1);
    end
    for J = 1:size(tp, 3)-1
        tp(:, :, J) = t(:, :, J) \ tp(:, :, J) * t(:, :, J+1);
    end

    % update neighbour transforms of extreme images
    tp(:, :, end) = tm(:, :, end);
    tm(:, :, 1) = tp(:, :, 1);

    % stopping criteria...
    % ... norm of the transform update different from identity
    tnorm = zeros(1, size(t, 3));
    for J = 1:size(t, 3)
        tnorm(J) = norm(t(:, :, J) - eye(size(t(:, :, 1))), 2);
    end
    info.MaxTNorm(I) = max(tnorm);
    info.MeanTNorm(I) = mean(tnorm);
    info.MedTNorm(I) = median(tnorm);
    if (all(tnorm <= opt.Epsilon))
        break
    end

    % ... maximum number of iterations
    if (I == opt.MaxIter)
        break
    end

end

info.NumIter = I;

end

%% check_bspline_injectivity_global:
%
% Check B-spline global injectivity triangulating the grid and looking for
% flipped triangles.
%
% tout: (N, K)-matrix, N number of slices, K number of B-spline control
%       points
% tbsp: elastix struct with the B-spline info (we need the grid spacing)
function isInjective = check_bspline_injectivity_global(tout, tbsp)

DEBUG = false;

% number of slices
N = size(tout, 1);

% grid of control points for 0 deformation
tbsp.TransformParameters(:) = 0;
[gx0, gy0] = elastix_bspline_grid(tbsp);
tri = delaunay(gx0(:), gy0(:));

% area of the triangles
atri = trifacet_signed_area(tri, [gx0(:), gy0(:)]);

% flip triangles with negative area, so that all are positive
idx = atri < 0;
tri(idx, :) = tri(idx, end:-1:1);

% DEBUG: plot control points
if (DEBUG)
    trimesh(tri, gx0(:), gy0(:))
end

isInjective = true(1, N);
for J = 1:N

    % transfer coefficients to the B-spline
    tbsp.TransformParameters = tout(J, :);

    % apply transformation to the control points (note:
    aux = transformix_pts(tbsp, [gx0(:) gy0(:)]);
    gx1 = reshape(aux(:, 1), size(gx0, 1), size(gx0, 2));
    gy1 = reshape(aux(:, 2), size(gy0, 1), size(gy0, 2));

    % DEBUG: plot control points
    if (DEBUG)
        trimesh(tri, gx1(:), gy1(:))
    end

    % area of the triangles
    atri = trifacet_signed_area(tri, [gx1(:), gy1(:)]);

    % if any triangle is flipped, the B-spline is not injective
    isInjective(J) = all(atri > 0);

end

end

%% check_bspline_injectivity_choi2000:
%
% Check B-spline injectivity using sufficient conditions in Choi and Lee
% (2000). Constrain control points to guarantee injectivity.
%
% This function checks Choi and Lee (2000) conditions for injectivity of
% the B-spline. The conditions are sufficient, i.e. if any is fulfilled,
% then we know the B-spline is injective. If not, the B-spline may or may
% not be injective, but we'll assume it's not.
%
% tout: (N, K)-matrix, N number of slices, K number of B-spline control
%       points
% tbsp: elastix struct with the B-spline info (we need the grid spacing)
%
% nConstrained: number of control points that had to be constrained to
%               guarantee injectivity
function [nConstrained, tout] = check_bspline_injectivity_choi2000(tout, tbsp)

% constants from Choi and Lee (2000), to check violations
K2 = 2.046392675;
A2 = sqrt((3/2)^2 + (K2 - 3/2)^2);

% number of slices
N = size(tout, 1);

% Note: the grid has row->x, col->y coordinates, instead of the Matlab
% convention, which is the opposite

% number of B-spline control points
L = length(tbsp.TransformParameters)/2;
Nx = tbsp.GridSize(1);
Ny = tbsp.GridSize(2);

% loop slices
isInjective = false(Nx, Ny, N);
for J = 1:N

    %% find control points that are guaranteed to not cause injectivity
    %% problems

    % transfer coefficients to the B-spline
    tbsp.TransformParameters = tout(J, :);

    % easy nomenclature for x, y coordinates of control points
    cx = tbsp.TransformParameters(1:L);
    cy = tbsp.TransformParameters(L+1:end);

    % reshape coefficients into grid
    cx = reshape(cx, Nx, Ny);
    cy = reshape(cy, Nx, Ny);

    % normalize control point displacement to grid spacing = 1
    cx = cx / tbsp.GridSpacing(1);
    cy = cy / tbsp.GridSpacing(2);

    % Theorem 1 from Choi and Lee (2000): first sufficient condition
    % for injectivity
    cond1 = (abs(cx) < 1 / K2) & (abs(cy) < (1 / K2));

    % Theorem 2 from Choi and Lee (2000): second sufficient condition
    % for injectivity
    cond2 = cx.^2 + cy.^2 < (1 / A2)^2;

    % if either condition is fulfilled, then the B-spline is injective
    isInjective(:, :, J) = cond1 | cond2;

    %% correct control points that potential can cause non-injectivity

    % find problematic points (linear index easier than row/column
    % index)
    idx = find(~isInjective(:, :, J));

    % correction factor for the problematic coefficients
    S = 1 ./ sqrt(cx(idx).^2 + cy(idx).^2) / (A2 - 0.01);

    % correct problematic control points
    cx(idx) = cx(idx) .* S;
    cy(idx) = cy(idx) .* S;

    % transfer corrected coefficients to the output
    aux = tout(J, :);
    aux(idx) = tbsp.GridSpacing(1)*cx(idx);
    aux(L + idx) = tbsp.GridSpacing(2)*cy(idx);
    tout(J, :) = aux;

end

% count number of control points that had to be constrained to guarantee
% they are injective
nConstrained = nnz(~isInjective);

end

%% check_bspline_injectivity_chun2009:
%
% Check B-spline injectivity using sufficient conditions in Chun and
% Fessler (2009).
%
% This function checks Chun and Fessler (2009) conditions for injectivity
% of the B-spline. If all the conditions are fulfilled, then we know the
% B-spline is injective. If not, the B-spline may or may not be injective,
% but we'll assume it's not.
%
% Because the conditions check grid edge lengths, it's not trivial to
% associate the conditions to vertices.
%
% tout: (N, K)-matrix, N number of slices, K number of B-spline control
%       points
% tbsp: elastix struct with the B-spline info (we need the grid spacing)
function isInjective = check_bspline_injectivity_chun2009(tout, tbsp)

DEBUG = false;

% constants from Chun and Fessler (2009), to check for violations
kx = 0.5 - 0.01;
ky = 0.5 - 0.01;

% number of slices
N = size(tout, 1);

% number of B-spline coefficients
L = length(tbsp.TransformParameters);
Nx = tbsp.GridSize(1);
Ny = tbsp.GridSize(2);

% DEBUG: plot control points
if (DEBUG)
    % grid of control points for 0 deformation
    tbsp.TransformParameters(:) = 0;
    [gx0, gy0] = elastix_bspline_grid(tbsp);
    tri = delaunay(gx0(:), gy0(:));

    trimesh(tri, gx0(:), gy0(:))
end

% loop slices
isInjective = true;
for J = 1:N

    % transfer coefficients to the B-spline
    tbsp.TransformParameters = tout(J, :);

    % easy nomenclature for x, y coordinates of control points
    cx = tbsp.TransformParameters(1:L/2);
    cy = tbsp.TransformParameters(L/2+1:end);

    % reshape coefficients into grid
    cx = reshape(cx, Nx, Ny);
    cy = reshape(cy, Nx, Ny);

    % transpose grid so that we have x->cols, y->rows
    cx = cx';
    cy = cy';

    % grid spacing
    mx = tbsp.GridSpacing(1);
    my = tbsp.GridSpacing(2);

    % conditions (i->x, j->y):

    % -mx * kx <= c_{i+1,j}^x - c_{i,j}^x
    % we also add a column to the right to account for the lost column
    cond = (diff(cx, 1, 2) >= -mx * kx);
    isInjective = isInjective && all(cond(:));

    % -my * ky <= c_{i,j+1}^y - c_{i,j}^y
    % we also add a row to the bottom to account for the lost row
    cond = (diff(cy, 1, 1) >= -my * ky);
    isInjective = isInjective && all(cond(:));

    % |c_{i,j+1}^x - c_{i,j}^x| <= mx * kx
    % we also add a row to the bottom to account for the lost row
    cond = (abs(diff(cx, 1, 1)) <= mx * kx);
    isInjective = isInjective && all(cond(:));

    % |c_{i+1,j}^y - c_{i,j}^y| <= my * ky
    % we also add a column to the right to account for the lost column
    cond = (abs(diff(cy, 1, 2)) <= my * ky);
    isInjective = isInjective && all(cond(:));

end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% pop_snapread() - load an EEG SnapMaster file (pop out window if no arguments).
%
% Usage:
%   >> [dat] = pop_snapread( filename, gain);
%
% Graphic interface:
%   "Relative gain" - [edit box] to compute the relative gain, fisrt look at
%                   the text header of the snapmater file with a text editor.
%                   Find the recording unit, usually in volts (UNITS field).
%                   Then, find the voltage range in the "CHANNEL.RANGE" [cmin cmax]
%                   field. Finally, determine the gain of the amplifiers (directly
%                   on the machine, not in the header file).
%                   Knowing that the recording precision is 12 bits. The folowing
%                   formula
%                                    1/2^12*[cmax-cmin]*1e6/gain
%                   returns the relative gain. You have to compute it and enter
%                   it in the edit box. Enter 1, for preserving the data file units.
%                   (note that if the voltage range is not the same for all channels
%                   or if the CONVERSION.POLY field in the file header
%                   is not "0 + 1x" for all channels,  you will have to load the data
%                   using snapread() and scale manually all channels, then import
%                   the Matlab array into EEGLAB).
%
% Inputs:
%   filename       - SnapMaster file name
%   gain           - relative gain. See graphic interface help.
%
% Outputs:
%   dat            - EEGLAB data structure
%
% Author: Arnaud Delorme, CNL/Salk Institute, 13 March 2002
%
% See also: eeglab(), snapread()

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

function [EEG, command] = pop_snapread(filename, gain);
command = '';
EEG = [];

if nargin < 1
	% ask user
	[filename, filepath] = uigetfile('*.SMA', 'Choose a SnapMaster file -- pop_snapread()');
	if filename == 0 return; end;
	filename = [filepath filename];

    promptstr    = { 'Relative gain (see help)' };
    inistr       = { '400' };
    result       = inputdlg2( promptstr, 'Import SnapMaster file -- pop_snapread()', 1,  inistr, 'pop_snapread');
    if length(result) == 0 return; end;
    gain   = eval( result{1} );

end;

if exist('gain') ~= 1
    gain = 1;
end;

% load datas
% ----------
EEG = eeg_emptyset;
[EEG.data,params,events, head] = snapread(filename);

EEG.data            = EEG.data*gain;
EEG.comments        = [ 'Original file: ' filename ];
EEG.filepath        = '';
EEG.setname 		= 'SnapMaster file';
EEG.nbchan          = params(1);
EEG.pnts            = params(2);
EEG.trials          = 1;
EEG.srate           = params(3);
EEG.xmin            = 0;

A = find(events ~= 0);
if ~isempty(A)
    EEG.event = struct( 'type', mattocell(events(A), [1], ones(1,length(events(A)))), ...
                        'latency', mattocell(A(:)', [1], ones(1,length(A))) );
end;

EEG = eeg_checkset(EEG, 'eventconsistency');
EEG = eeg_checkset(EEG, 'makeur');
command = sprintf('EEG = pop_snapread(''%s'', %f);', filename, gain);

return;

author: PatternRecognition

repo: OpenBMI


text:
function CoreRunTests(testsOrDirectory)
    % CoreRunTests Runs MTest unit tests for Matlab
    %
    % Executes the specified unit tests and reports on the sucess and
    % failures.
    %
    % Syntax:
    %     CoreRunTests(testsOrDirectory);
    %
    %         testsOrDirectory - a folder or a cell array of strings, each containing the name of a test to run
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if (nargin > 0) && ~isempty(testsOrDirectory)
        if iscell(testsOrDirectory)
            testClasses = testsOrDirectory;
        elseif CoreDiskUtilities.DirectoryExists(testsOrDirectory)
            testClasses = CoreDiskUtilities.GetListOfClassFiles(testsOrDirectory, 'CoreTest');
        elseif ischar(testsOrDirectory) && exist(testsOrDirectory, 'class') == 8
            testClasses = {testsOrDirectory};
        else
            throw(MException('CoreRunTests:UnrecognisedArgument', 'Please specify a test directory or a set of test class names to run'));
        end
        RunTests(testClasses);
    else
        disp('Please specify a test directory or a set of test class names to run');
    end
end

function RunTests(testClasses)
    numTests = numel(testClasses);
    if numTests == 1
        testString = 'test';
    else
        testString = 'tests';
    end
    disp(['Running ' testString '...']);
    passed = 0;
    failed = 0;

    for testClass = testClasses
        testClassName = testClass{1};
        try
            classHandle = feval(testClassName);
            classHandle.delete();
            disp([' + Passed: ' char(testClassName)]);
            passed = passed + 1;
        catch ex
            disp(['***** FAILED: ' char(testClassName)]);
            failed = failed + 1;
        end
    end
    disp(['Complete. ' int2str(failed+passed) ' ' testString ' run. ' int2str(passed) ' passed. ' int2str(failed) ' failed. ' ]);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [tags] = sequence_get_tags(sequence, index)
% sequence_get_tags Returns all tags for a given frame
%
% Input:
% - sequence (structure): A valid sequence structure.
% - index (integer): A index of a frame.
%
% Output:
% - tags (cell): A cell array of strings. Names of tags for the given frame.

if (sequence.length < index || index < 1)
    tags = {};
else
    tags = sequence.tags.names(logical(sequence.tags.data(index, :)));
end;

author: votchallenge

repo: toolkit-legacy


text:
function classSymbols = prtPlotUtilClassSymbolsUnlabeled
% xxx Need Help xxx

classSymbols = '.';

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function epo= proc_appendEpochs(epo1, epo2)
%PROC_APPENDEPOCHS - appends two or more structs of epoched data.
%
%Synopsis:
% EPO= proc_appendEpochs(EPO1, EPO2)
% EPO= proc_appendEpochs(EPOCELL)
%
%Arguments:
% EPOCELL         - cell array of epoched data
% EPO1, EPO2      - epoched data to be appended
%
%Returns:
% EPO             - appended epoched data
%
%Description:
% Appends the epochs of 'EPO2' to 'EPO1' (or appends all epochs of ERPCELL)
%
% SEE  proc_segmentation, mrk_mergerMarkers, proc_appendCnt,
%      proc_appendChannels

%epo1 = misc_history(epo1);

misc_checkType(epo1, 'STRUCT(x clab)|CELL{STRUCT}');
misc_checkTypeIfExists('epo2', 'STRUCT(x clab)');

if iscell(epo1),
  Cepo= epo1;
  epo1= Cepo{1};
  for ii= 2:length(Cepo),
    epo1= proc_appendEpochs(epo1, Cepo{ii});
  end
  epo=epo1;
  return
end

if isempty(epo1),
  epo= epo2;
  return;
elseif isempty(epo2),
  epo= epo1;
  return;
end

% begin sthf
if isfield(epo1, 'ndims')
  ndim = max(3, epo1.ndims);
else
  ndim = max(3, ndims(epo1.x));
end

if size(epo1.x,ndim-2)~=size(epo2.x,ndim-2),
  error('interval length mismatch');
end
if size(epo1.x,ndim-1)~=size(epo2.x,ndim-1),
  error('number of channels mismatch');
end
if ( ndim == 4 ) && ( size(epo1.x,ndim-3)~=size(epo2.x,ndim-3) ),
  error('number of frequencies mismatch');
end

epo= mrkutil_appendMarkersExcludingTime(epo1, epo2);
epo.x= cat(ndim, epo1.x, epo2.x);

% We could add functionality to struct_areFieldsEqual to return as second
% argument, which fields do not match. Then we could give here a more
% precise warning.
if ~struct_areFieldsEqual(epo1, epo2, {'fs','clab', 't'})
  warning('epochs are inconsistent wrt ''fs'', ''clab'', or ''t''.');
end
epo= struct_copyFields(epo, epo1, {'fs','clab','t'});

further_fields= {'xUnit', 'yUnit'};
for Fld= further_fields,
  if ~struct_areFieldsEqual(epo1, epo2, Fld)
    warning(sprintf('epochs are inconsistent wrt ''%s''.', Fld{1}));
  end
  epo= struct_copyFields(epo, epo1, Fld);
end

if isfield(epo1, 'mrk_info') && isfield(epo2, 'mrk_info'),
  epo.mrk_info= mrkutil_appendEventInfo(epo1.mrk_info, epo2.mrk_info);
end

% The following code should be superflutious in future
if isfield(epo1, 'p') && isfield(epo2, 'p')
  epo.p= cat(ndim, epo1.p, epo2.p);
end
if isfield(epo1, 'sgnlogp') && isfield(epo2, 'sgnlogp')
  epo.sgnlogp= cat(ndim, epo1.sgnlogp, epo2.sgnlogp);
end
if isfield(epo1, 'se') && isfield(epo2, 'se')
  epo.se= cat(ndim, epo1.se, epo2.se);
end
if isfield(epo1, 'sigmask') && isfield(epo2, 'sigmask')
  epo.sigmask= cat(ndim, epo1.sigmask, epo2.sigmask);
end
if isfield(epo1, 'crit') && isfield(epo2, 'crit')
  epo.crit= cat(1, epo1.crit, epo2.crit);
end
if isfield(epo1, 'df') && isfield(epo2, 'df')
  epo.df= cat(1, epo1.df, epo2.df);
end

author: bbci


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
%
% Create random divisions of the training set for cross-validation
%
% ------------------------------------------------------------------------

% Create random divisions of the training set
id = '3';
percentage_train = 0.5;
full_set = 'train2012';

% Get full ids
im_ids = database_ids('pascal2012',full_set);

ids_a = sort(randperm(length(im_ids),floor(percentage_train*length(im_ids))));
ids_b = setdiff(1:length(im_ids),ids_a);
im_ids_a = im_ids(ids_a);
im_ids_b = im_ids(ids_b);

% Write to file
file_a = fullfile(mcg_root,'datasets', 'pascal2012','gt_sets',[full_set '_' id 'a.txt']);
file_b = fullfile(mcg_root,'datasets', 'pascal2012','gt_sets',[full_set '_' id 'b.txt']);

dlmwrite(file_a,im_ids_a,'')
dlmwrite(file_b,im_ids_b,'')

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function varargout = subsref(F, ref)
%SUBSREF   SPHEREFUNV subsref.
%
% ( )
%   F(LAM,TH) returns the values of the SPHEREFUNV F evaluated on the array
%   (LAM,TH) in spherical coordinates.
%
%   F(X,Y,Z) returns the values of the SPHEREFUNV F evaluated on the array
%   (X,Y,Z) in Cartesian coordinates.
%
%   F(k) returns the first component of F if k=1, the second if k=2, and
%   the third if k=3.
%  .
%   F.PROP returns the property PROP of F as defined by GET(F,'PROP').
%
% { }
%    Throws an error.

% Check for empty SPHEREFUNV object.
if ( isempty(F) )
   varargout = {[]};
   return
end

indx = ref(1).subs;

switch ( ref(1).type )

    case '.'
        if ( numel(ref) == 1 )
            % This is a get call to get a property.
            varargout = { get(F, indx) };
        else
            t2 = ref(2).type;
            if ( strcmp(t2,'.') )
                out = get(F, indx, ref(2).subs{:});
            else
                out = get(F, indx);
                out = out(ref(2).subs{:});
            end
            if ( numel(ref) > 2 )
                varargout = {subsref(out, ref(3:end))};
            else
                varargout = { out };
            end
        end

    case '()'
        if ( length(indx) == 2 ) % spherical coordinates
            lam = indx{1};
            th = indx{2};
            vals = feval(F, lam, th);
            varargout = { vals };
        elseif ( length(indx) == 3 ) % Cartesian coordinates
            x = indx{1};
            y = indx{2};
            z = indx{3};
            vals = feval(F, x, y, z);
            varargout = { vals };
        else
            if ( isa(indx{1},'double') )
                if all( indx{1} == 1  )
                    varargout = F.components(1);
                elseif ( all( indx{1} == 2 ) )
                    varargout = F.components(2);
                elseif ( ( all(indx{1} == 3) )  )
                    varargout = F.components(3);
                else
                    error('SPHEREFUN:SPHEREFUNV:subsref:index', ...
                        'SPHEREFUNV only contains three components.');
                end
            end
        end

    otherwise
        error('SPHEREFUN:SPHEREFUNV:subsref:unexpectedType', ...
            ['??? Unexpected index.type of ' index(1).type]);

end

% Recurse down:
if ( numel( ref ) > 1 )
   ref(1) = [];
   varargout = { subsref( varargout{ : }, ref ) };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Max filter

function y = Maxfilter(x,filter_len)

[nf,N_ch] = size(x);
hs = (filter_len-1)/2; % half side of the filter
y=x;
y(1:hs,:) = max(x(1:hs,:));
for i=hs+1:nf-hs
    y(i,:) = max(x(i-hs:i+hs,:));
end
y(end-hs+1:end,:) = max(x(end-hs+1:end,:));

author: singaxiong

repo: SignalGraph


text:
function [] = plot_results(setup_file)

    width = 24;
    height = 12;
    fontsize = 9;

    setup = ReadYaml(setup_file)

    algorithm = setup.algorithm;
    configuration_name = setup.configuration_name;
    experiment_name = setup.experiment_name;
    granularity = setup.granularity;
    dataset = setup.dataset;
    household = setup.household;
    setup_name = setup.setup_name;

    source_folder = ['results/details/', algorithm, '/', configuration_name, '/', experiment_name, '/', setup_name];
    plot_folder = ['results/plots/', algorithm, '/', configuration_name, '/', experiment_name, '/', setup_name, '/'];

    load([source_folder, '/result1.mat']);

    if ~exist(plot_folder)
        mkdir(plot_folder)
    end

    evaluation_days = result.evaluation_and_training_days{1};
    first_day = evaluation_days(1,:);
    num_days = size(evaluation_days, 1);

    total_consumption = read_smartmeter_data(dataset, household, evaluation_days, granularity, 'powerallphases');

    %%%%%%%%%%%%%%
    %% Plot events
    %%%%%%%%%%%%%%
    if isfield(result, 'events')
        i = 1;
        for appliance_name_cell = result.appliance_names
            appliance_name = cell2mat(appliance_name_cell);
            appliance_id = getApplianceID(appliance_name);
            fprintf(appliance_name);
            plug_consumption = read_plug_data(dataset, household, appliance_id, evaluation_days, granularity);

            classified_events = struct;

            % events inferred from algorithm
            % col1: event time
            % col2: appliance id
            % col3: delta
            events_inferred = result.events(result.events(:,2) == i,:);

            % get ground truth from plug-level data
            filter = struct;
            filter.p_edgeThreshold = setup.p_edgeThreshold;
            filter.p_filtering = setup.p_filtering;
            filter.p_filtLength = setup.p_filtLength;
            % col1: time start
            % col2: time stop
            % col3: consumption delta
            events_from_plug_data = getEventsFromPlugData(dataset, household, appliance_id, evaluation_days, granularity, filter);

            % class one events: correctly inferred
            % class two events: missed events
            plug_events_detected = arrayfun(@(idx) any(events_inferred(:,1) > events_from_plug_data(idx,1) - 20 & ...
                    events_inferred(:,1) < events_from_plug_data(idx,1) + 20), [1:size(events_from_plug_data, 1)]');
            avg_deviation_of_plug_data_event = mean(events_from_plug_data(:,1) - events_from_plug_data(:,2));
            classified_events.tp_events = [events_from_plug_data(plug_events_detected,1) + avg_deviation_of_plug_data_event, events_from_plug_data(plug_events_detected,3)];
            classified_events.fn_events = [events_from_plug_data(~plug_events_detected,1) + avg_deviation_of_plug_data_event, events_from_plug_data(~plug_events_detected,3)];

            % class three events: incorrectly inferred
            which_detected_events_have_ground_truth_events = arrayfun(@(t) any(events_from_plug_data(:,1) - 20 < t & t < events_from_plug_data(:,2) + 20), ...
                    events_inferred(:, 1));
            plug_data_exists = arrayfun(@(t) ~any(plug_consumption(max(1,t-5):min(t+5,length(plug_consumption))) == -1), events_inferred(:, 1));
            tmp_idx = ~which_detected_events_have_ground_truth_events & plug_data_exists;
            classified_events.fp_events = events_inferred(tmp_idx, [1,3]);

            for d = 1:num_days
                day = evaluation_days(d,:);
                idx_start = 86400 / granularity * (d-1) + 1;
                idx_stop = 86400 / granularity * d;
                date_range = datenum(day)+granularity/86400 : granularity/86400 : datenum(day)+1;

                % plot smart meter consumption
                fig = figure;
                hold on;
                plot(date_range, total_consumption(idx_start:idx_stop));
                max_cons = max(total_consumption(idx_start:idx_stop));

                % get events for current day
                tp_events = classified_events.tp_events(:,1) >= idx_start & classified_events.tp_events(:,1) <= idx_stop;
                tp_events_on = classified_events.tp_events(classified_events.tp_events(:,2) > 0 & tp_events, 1);
                tp_events_off = classified_events.tp_events(classified_events.tp_events(:,2) < 0 & tp_events, 1);
                fn_events = classified_events.fn_events(:,1) >= idx_start & classified_events.fn_events(:,1) <= idx_stop;
                fn_events_on = classified_events.fn_events(classified_events.fn_events(:,2) > 0 & fn_events, 1);
                fn_events_off = classified_events.fn_events(classified_events.fn_events(:,2) < 0 & fn_events, 1);
                fp_events = classified_events.fp_events(:,1) >= idx_start & classified_events.fp_events(:,1) <= idx_stop;
                fp_events_on = classified_events.fp_events(classified_events.fp_events(:,2) > 0 & fp_events, 1);
                fp_events_off = classified_events.fp_events(classified_events.fp_events(:,2) < 0 & fp_events, 1);

                for ev = 1:size(tp_events_on, 1)
                    time = tp_events_on(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], 'g-');
                end

                for ev = 1:size(tp_events_off, 1)
                    time = tp_events_off(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], 'g--');
                end

                for ev = 1:size(fn_events_on, 1)
                    time = fn_events_on(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], '-', 'Color', [0.5 0.5 0.5]);
                end

                for ev = 1:size(fn_events_off, 1)
                    time = fn_events_off(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], '--', 'Color', [0.5 0.5 0.5]);
                end

                for ev = 1:size(fp_events_on, 1)
                    time = fp_events_on(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], 'r-');
                end

                for ev = 1:size(fp_events_off, 1)
                    time = fp_events_off(ev, 1);
                    xval = datenum(day) + time/86400*granularity - d + 1;
                    plot([xval,xval],[0,max_cons], 'r--');
                end

                ylim([0 max_cons]);
                datetick('x','HHPM')
                % legend('Actual', 'Inferred'); % does not work
                title([appliance_name, ' - ', evaluation_days(d,:)])
                fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
                filename = ['Events_', appliance_name, '_', evaluation_days(d,:)];
                % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
                saveas(fig, [plot_folder, filename, '.png'], 'png');
                close(fig);
            end
            i = i+1;
        end
    end

    %%%%%%%%%%%%%
    %% Plot inferred consumption and pie chart
    %%%%%%%%%%%%%
    if isfield(result, 'consumption')
        %% plot inferred plug data vs ground truth plug data
        i = 1;
        for appliance_name_cell = result.appliance_names
            appliance_name = cell2mat(appliance_name_cell);
            appliance_id = getApplianceID(appliance_name);
            plug_consumption = read_plug_data(dataset, household, appliance_id, evaluation_days, granularity);
            inferred_consumption = result.consumption(i,:);
            for d = 1:num_days
                day = evaluation_days(d,:);
                idx_start = 86400 / granularity * (d-1) + 1;
                idx_stop = 86400 / granularity * d;
                date_range = datenum(day)+granularity/86400 : granularity/86400 : datenum(day)+1;

                fig = figure;
                hold on;
                plot(date_range, plug_consumption(idx_start:idx_stop));
                plot(date_range, inferred_consumption(idx_start:idx_stop), 'r');
                if strcmp(appliance_name, 'Fridge') == 1 || strcmp(appliance_name, 'Freezer') == 1
                    ylim([0 500]);
                end
                datetick('x','HHPM')
                % legend('Actual', 'Inferred'); % does not work
                title([appliance_name, ' - ', evaluation_days(d,:)])
                fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
                filename = ['Consumption_', appliance_name, '_', evaluation_days(d,:)];
                % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
                saveas(fig, [plot_folder, filename, '.png'], 'png');
                close(fig);
            end
            i = i+1;
        end

        %% plot ground truth smart meter data vs ground truth plug data and vs. ground truth inferred consumption
        i = 1;
        plug_consumption = {};
        inferred_consumption = {};
        for appliance_name_cell = result.appliance_names
            appliance_name = cell2mat(appliance_name_cell);
            appliance_id = getApplianceID(appliance_name);
            plug_consumption{i} = read_plug_data(dataset, household, appliance_id, evaluation_days, granularity);
            inferred_consumption{i} = result.consumption(i,:);
            i = i+1;
        end
        for d = 1:num_days
            day = evaluation_days(d,:);
            idx_start = 86400 / granularity * (d-1) + 1;
            idx_stop = 86400 / granularity * d;
            date_range = datenum(day)+granularity/86400 : granularity/86400 : datenum(day)+1;

            fig = figure;
            hold on;
            plot(date_range, total_consumption(idx_start:idx_stop));
            for j = 1:length(plug_consumption)
                 plot(date_range, plug_consumption{j}(idx_start:idx_stop));
            end
            datetick('x','HHPM')
            % legend('Actual', 'Inferred'); % does not work
            title(['Plug and Smart Meter - ', evaluation_days(d,:)])
            fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
            filename = ['plug-and-smartmeter_', evaluation_days(d,:)];
            % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
            saveas(fig, [plot_folder, filename, '.png'], 'png');
            close(fig);
            fig = figure;
            hold on;
            plot(date_range, total_consumption(idx_start:idx_stop));
            for j = 1:length(plug_consumption)
                plot(date_range, inferred_consumption{j}(idx_start:idx_stop));
            end
            datetick('x','HHPM')
            % legend('Actual', 'Inferred'); % does not work
            title(['Inferred and Smart Meter -  - ', evaluation_days(d,:)])
            fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
            filename = ['Smartmeter_', evaluation_days(d,:)];
            % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
            saveas(fig, [plot_folder, filename, '.png'], 'png');
            close(fig);
        end

        %% plot pie charts
        i = 1;
        plug_consumption = {};
        inferred_consumption = {};
        for appliance_name_cell = result.appliance_names
            appliance_name = cell2mat(appliance_name_cell);
            appliance_id = getApplianceID(appliance_name);
            tmp = read_plug_data(dataset, household, appliance_id, evaluation_days, granularity);
            valid_time = tmp ~= -1;
            plug_consumption{i} = mean(tmp(valid_time)) * num_days * 24;
            cons = result.consumption(i,:);
            inferred_consumption{i} = mean(cons(valid_time)) * num_days * 24;
            i = i+1;
        end
        other_plugs = mean(total_consumption) * num_days * 24 - sum(cell2mat(plug_consumption));
        other_inferred = mean(total_consumption) * num_days * 24 - sum(cell2mat(inferred_consumption));
        fig = figure;
        hold on;
        subplot(1,2,1);
        pie([cell2mat(plug_consumption), other_plugs]);
        title('Ground truth');
        subplot(1,2,2);
        pie([cell2mat(inferred_consumption), other_inferred]);
        title('Inferred');
        fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
        filename = 'Pie_Chart';
        % print('-depsc2', '-cmyk', '-r600', [plot_folder, filename, '.eps']); % if eps is needed
        saveas(fig, [plot_folder, filename, '.png'], 'png');
        close(fig);
    end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function value = g_hofstadter ( n )

%*****************************************************************************80
%
%% G_HOFSTADTER computes the Hofstadter G sequence.
%
%  Discussion:
%
%    G(N) = 0                      if N = 0
%         = N - G ( G ( N - 1 ) ), otherwise.
%
%    G(N) is defined for all nonnegative integers.
%
%    The value of G(N) turns out to be related to the Zeckendorf
%    representation of N as a sum of non-consecutive Fibonacci numbers.
%    To compute G(N), determine the Zeckendorf representation:
%
%      N = sum ( 1 <= I <= M ) F(I)
%
%    and reduce the index of each Fibonacci number by 1:
%
%      G(N) = sum ( 1 <= I <= M ) F(I-1)
%
%    However, this is NOT how the computation is done in this routine.
%    Instead, a straightforward recursive function call is defined
%    to correspond to the definition of the mathematical function.
%
%  Table:
%
%     N  G(N)  Zeckendorf   Decremented
%    --  ----  ----------   -----------
%
%     1   1    1            1
%     2   1    2            1
%     3   2    3            2
%     4   3    3 + 1        2 + 1
%     5   3    5            3
%     6   4    5 + 1        3 + 1
%     7   4    5 + 2        3 + 1
%     8   5    8            5
%     9   6    8 + 1        5 + 1
%    10   6    8 + 2        5 + 1
%    11   7    8 + 3        5 + 2
%    12   8    8 + 3 + 1    5 + 2 + 1
%    13   8    13           8
%    14   9    13 + 1       8 + 1
%    15   9    13 + 2       8 + 1
%    16  10    13 + 3       8 + 2
%    17  11    13 + 3 + 1   8 + 2 + 1
%    18  11    13 + 5       8 + 3
%    19  12    13 + 5 + 1   8 + 3 + 1
%    20  12    13 + 5 + 2   8 + 3 + 1
%    21  13    21           13
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Douglas Hofstadter,
%    Goedel, Escher, Bach,
%    Basic Books, 1979.
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%
%    Output, integer VALUE, the value of the function.
%
  if ( n <= 0 )
    value = 0;
  else
    value = n - g_hofstadter ( g_hofstadter ( n-1 ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Obs = createObservations(Sen,Opt)

% CREATEOBSERVATIONS Create Obs structure array.
%   Rob = CREATEOBSERVATIONS(Sen,Lmk) creates the Obs() structure array to
%   be used as SLAM data, from the information contained in Sen() and
%   Rob(). See the toolbox documentation for details on this structure.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

for sen = 1:numel(Sen)

    for lmk = 1:Opt.map.numLmks

        S = Sen(sen);
%         L = Lmk(lmk);

        O.sen       = sen;        % sensor index
        O.lmk       = lmk;        % landmark index.
        O.sid       = S.id;       % sensor id
        O.lid       = [];         % lmk id
        O.stype     = S.type;     % sensor type
        O.ltype     = '';         % lmk type
        O.meas.y    = [];         % observation
        O.meas.R = S.par.cov; % observation cov
        O.nom.n     = [];         % expected non-observable
        O.nom.N     = [];         % expected non-observable cov
        O.exp.e     = [];         % expectation mean
        O.exp.E     = [];         % expectation cov
        O.exp.um    = [];         % expectation uncertainty measure
        O.inn.z     = [];         % innovation
        O.inn.Z     = [];         % innovation cov
        O.inn.iZ    = [];         % inverse inn. cov
        O.inn.MD2   = 0;          % Mahalanobis distance squared
        O.app.pred  = [];         % predicted appearence
        O.app.curr  = [];         % current appearence
        O.app.sc    = 0;          % match score
        O.par       = [];         % other params
        O.vis       = false;      % lmk is visible?
        O.measured  = false;      % lmk has been measured?
        O.matched   = false;      % lmk has been matched?
        O.updated   = false;      % lmk has been updated?
        O.Jac.E_r   = [];         % Jac of expectation wrt robot state - Obs function
        O.Jac.E_s   = [];         % Jac of expectation wrt sensor
        O.Jac.E_l   = [];         % Jac of expectation wrt landmark
        O.Jac.Z_r   = [];         % Jac of innovation wrt robot state - Inn function
        O.Jac.Z_s   = [];         % Jac of innovation wrt sensor
        O.Jac.Z_l   = [];         % Jac of innovation wrt landmark

        Obs(sen,lmk) = O;

    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Launch code for static input (sample)               %
%      code by Fabrizio Conso, university of pavia, student        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                launch code example                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                       %
 format short;
 f_s=1;            % normalized sampling frequency
 f_bw=5;                    % DAC's bandwidth
 sr=2;                      % DAC's slew-rate
 nbit=12;                   % converter bits
 input=0.7;                 % sample value

 [counter,tresholds]  =Approx(input,nbit,f_bw,sr,f_s)
 [counter2,tresholds2]=Approx_2(input,nbit,f_bw,sr,f_s)
 bitstream=dec2bin(counter)
 bitstream2=dec2bin(counter2)
%                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
%% retraction of d in T_X(M_r) onto M_r
%   X=U*diag(S)*V'
%   d=U*M*V'+U_p*V'+U*V_p'
%   cf. [Van12]

function [S1,U1,V1]=rtr_fr(d,S,U,V)

n1=size(U,1);
n2=size(V,1);
n3=size(S,1);

sig=0;
if size(S,2)>1, S=diag(S); sig=1; end

eps=1e-3;

% Pu=U*U'; Pv=V*V';

M=U'*d*V;
% U_p=d*V-Pu*d*V; V_p=d'*U-Pv*d'*U;
U_p=d*V-U*M; V_p=d'*U-V*M';

[Q_u,R_u]=qr(U_p,0);
[Q_v,R_v]=qr(V_p,0);
M1=[diag(S)+M, R_v'; R_u, zeros(n3)];
[U1,S1,V1]=svd(M1);
S1=diag(S1); S1=S1(1:n3);
U1=[U,Q_u]*U1(:,1:n3);
V1=[V,Q_v]*V1(:,1:n3);

S1=max(S1,eps);

if sig==1, S1=diag(S1); end

author: andrewssobral

repo: lrslibrary


text:
function plotProjAmpVsPosition(view)
%
% plotProjAmpVsPosition(view)
%
% Plot of projected amplitude versus linear position for the current scan, for all
% pixels in the current ROI.  (The current ROI should be a line ROI,
% otherwise, this plot doesn't make much sense.)
%
%
% fwc   12/07/02    adapted to projected amplitude vs position

% Get selpts from current ROI
if view.selectedROI
  ROIcoords = getCurROIcoords(view);
else
  myErrorDlg('No current ROI');
end

% Get co and ph (vectors) for the current scan, within the
% current ROI.
%
curScan = getCurScan(view);
ph = getCurDataROI(view,'ph',curScan,ROIcoords);
amp = getCurDataROI(view,'amp',curScan,ROIcoords);
co = getCurDataROI(view,'co',curScan,ROIcoords);

% Remove NaNs from ph that may be there if ROI
% includes volume voxels where there are no data.
NaNs = find(isnan(co));
if ~isempty(NaNs)
  myWarnDlg('ROI includes voxels that have no data.  These voxels are being ignored.');
  notNaNs = find(~isnan(co));
  co = co(notNaNs);
  ph = ph(notNaNs);
  amp = amp(notNaNs);
end

% Compute the amplitude projected onto the phase of first pixel
projectedAmp = amp.*cos(ph-ph(1));

% Figure out the x-axis by the coordinates in ROIcoords
dLinePos = diff(ROIcoords');
dx = sqrt(dLinePos(:,1).^2 + dLinePos(:,2).^2);
x = [0;cumsum(dx)];

y = ph;
y = projectedAmp;

selectGraphWin

% Window header
ROIname = view.ROIs(view.selectedROI).name;
headerStr = ['Projected Amplitude vs. position, ROI ',ROIname,', scan ',num2str(curScan)];
set(gcf,'Name',headerStr);

% Plot it
fontSize = 14;
symbolSize = 4;

clf
h = plot(x, y, 'b-', x, y, 'bo','MarkerSize', symbolSize);
set(h, 'MarkerFaceColor', 'b')
ylabel('Projected Amplitude','FontSize',fontSize);
xlabel('Distance along flat line roi(pixels)');
% set(gca,'ylim',[-pi pi]);
% set(gca,'ylim',[0 2*pi]);

LEGEND(['scan ',num2str(curScan)], -1);

set(gca,'FontSize',fontSize)

% Save the data in gca('UserData')
data.position = x;
data.amp = y;
set(gca,'UserData',data);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function repeatability_generate_images()
  load('data/pascal_voc07_test_annotations.mat', 'impos');
  config = get_config();
  rep_config = repeatability_get_config();

  homographies = [];
  fields = fieldnames(rep_config);
  for field_i = 1:numel(fields)
    homographies.(fields{field_i}) = zeros(numel(impos), 3, 3);
  end

  for im_i=1:numel(impos)
    tic_toc_print('image %d/%d\n', im_i, numel(impos));
    orig_img_id = impos(im_i).im;
    im = imread(sprintf(config.pascal_images, orig_img_id));

    fields = fieldnames(rep_config);
    for field_i = 1:numel(fields)
      %tic_toc_print('running %s\n', fields{field_i});

      sub_config = rep_config.(fields{field_i});
      for i = 1:numel(sub_config.params), param = sub_config.params(i);
        % skip parameters that are not plotted
        if isfield(sub_config, 'display_points')
          if ~ismember(i, sub_config.display_points)
            continue;
          end
        end
        [transformed_im, H] = sub_config.func(im, param);
        homographies.(fields{field_i})(im_i,:,:) = H;

        img_id = sprintf(sub_config.img_id, orig_img_id, i);
        imwrite(transformed_im, sprintf(config.transformed_pascal_images, img_id));
      end
    end
  end

  save('data/pascal_voc07_test_repeatability_homographies.mat', 'homographies');
end

author: hosang

repo: detection-proposals


text:
function perm=field2d(Nx,Ny,k_avg,V_dp,clx,cly)
% All credits go to http://www.mysimlabs.com/surface_generation.html
% modified by Ali A. Eftekhari
% See the license file
% Note: a very very simple approach! Use it with utmost care :-)
Lx=1.0; % domain length
x=linspace(-Lx/2.0,Lx/2.0,Nx);
y=linspace(-Lx/2.0,Lx/2.0,Ny);
[X,Y]=ndgrid(x,y);
s=-log(1-V_dp); % standard deviation
mu=log(k_avg)-s*s/2.0; % mean for a log-random field
Z=s*randn(Nx,Ny); % normal distribution
F = exp(-(X.^2/(clx*clx/2.0)+Y.^2/(cly*cly/2.0)));
% Gaussian filter
f =2.0/sqrt(pi)*Lx/sqrt(Nx*Ny)/sqrt(clx)/sqrt(cly).*ifft2(fft2(Z).*fft2(F)); % another filter
perm=exp(mu+real(f)); % perm field

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function A = double(T)
%DOUBLE Convert a sptenmat to a sparse matrix.
%
%   A = double(T) converts T stored as a SPTENMAT to a sparse matrix.
%
%   See also SPTENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

m = prod(T.tsize(T.rdims));
n = prod(T.tsize(T.cdims));
if isempty(T.subs)
    A = sparse(m,n);
else
    A = sparse(T.subs(:,1), T.subs(:,2), T.vals, m, n);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% write_meta_tri.m
%
% write *.meta file
%
% Wan Jing
% 11/08/2008 - create
%============================================
function write_meta_tri(filename, vertices, faces)

fid = fopen(filename, 'w');
fprintf(fid, 'ObjectType = Scene\n');
fprintf(fid, 'NDims = 3\n');
fprintf(fid, 'NObjects = 1\n');
fprintf(fid, 'ObjectType = Mesh\n');
fprintf(fid, 'NDims = 3\n');
fprintf(fid, 'BinaryData = False\n');
fprintf(fid, 'TransformMatrix = 1 0 0 0 1 0 0 0 1\n');
fprintf(fid, 'Offset = 0 0 0\n');
fprintf(fid, 'CenterOfRotation = 0 0 0\n');
fprintf(fid, 'ElementSpacing = 1 1 1\n');
fprintf(fid, 'PointType = MET_FLOAT\n');
fprintf(fid, 'PointDataType = MET_FLOAT\n');
fprintf(fid, 'CellDataType = MET_FLOAT\n');
fprintf(fid, 'NCellTypes = 1\n');
fprintf(fid, 'PointDim = ID x y ...\n');
fprintf(fid, 'NPoints = %d\n', length(vertices));
fprintf(fid, 'Points = \n');
for i = 0:(length(vertices)-1)
    fprintf(fid,'%d %g %g %g \n', i, vertices(i+1,1), vertices(i+1,2), vertices(i+1, 3));
end
fprintf(fid, 'CellType = TRI\n');
nfaces = 2*length(faces);
fprintf(fid, 'NCells = %d\n', nfaces);
fprintf(fid, 'Cells = \n');
for i = 0:(length(faces)-1)
    fprintf(fid, '%d %d %d %d \n', 2*i, (faces(i+1,1)-1), (faces(i+1,2)-1), (faces(i+1,3)-1));
    fprintf(fid, '%d %d %d %d \n', 2*i+1,(faces(i+1,3)-1), (faces(i+1,4)-1), (faces(i+1,1)-1));
end
fclose(fid);

return;

author: vigente

repo: gerardus


text:
function testRansac

% real model coef
k = .5;
b = 10;
ptNum = 200;
outlrRatio = .4;
inlrStd = 5;
pts = genRansacTestPoints(ptNum,outlrRatio,inlrStd,[k b]);
figure,plot(pts(1,:),pts(2,:),'.'),hold on
X = -ptNum/2:ptNum/2;
plot(X,k*X+b,'k')
err0 = sqrError(k,b,pts(:,1:ptNum*(1-outlrRatio)))

% RANSAC
iterNum = 300;
thDist = 2;
thInlrRatio = .1;
[t,r] = ransac(pts,iterNum,thDist,thInlrRatio);
k1 = -tan(t);
b1 = r/cos(t);
plot(X,k1*X+b1,'r')
err1 = sqrError(k1,b1,pts(:,1:ptNum*(1-outlrRatio)))

% least square fitting
coef2 = polyfit(pts(1,:),pts(2,:),1);
k2 = coef2(1);
b2 = coef2(2);
plot(X,k2*X+b2,'g')
err2 = sqrError(k2,b2,pts(:,1:ptNum*(1-outlrRatio)))

end

function err = sqrError(k,b,pts)
%	Calculate the square error of the fit

theta = atan(-k);
n = [cos(theta),-sin(theta)];
pt1 = [0;b];
err = sqrt(sum((n*(pts-repmat(pt1,1,size(pts,2)))).^2));

end


text:
function view = clearData(view)
%
% view = clearData(view)
%
%AUTHOR:  Wandell
%DATE:    12.29.00
%PURPOSE:
%   Clear the data fields in a view
%
% $Author: wade $
% $Date: 2002/09/30 23:43:18 $
global mrSESSION
view.co = [];
view.ph = [];
view.amp = [];
view.map = [];
view.mapName = '';
view.tSeries = [];

% we need to clear the rm field as well, otherwise rm-definition can
% transfer (unwantedly) from datatype to datatype.
if isfield(view,'rm')
    view = rmfield(view,'rm');
end
return;


text:
function SStar=covMatEstShrinkage(x,centerData,algorithm)
%%COVMATESTSHRINKAGE When estimating a covariance matrix from a small
%           number of samples, it can occur that the resultng estimate is
%           poorly conditioned. This function estimates the covariance
%           matrix from samples using one of two shrinkage estimators that
%           are mean to avoid poor conditioning and that can produce more
%           accurate (in terms of the squared Frobenius norm) estimates
%           than one can obtain from the sample covariance matrix.
%
%INPUTS: x A pXn matrix of n samples of a p-dimensional distribution.
% centerData If the samples in x are zero mean or the mean has been
%          subtracted from them, then centerData can be false. Otherwise,
%          centerData should be true to center the data (subtract the
%          sample mean). The algorithmic choices assume that the true mean
%          has been subtracted, so subtracting the sample mean might affect
%          the performance. The default if omitted or an empty matrix is
%          passd is true.
% algorithm An optional parameter selecting the algorithm to use. Possible
%          values are:
%          0 (The default if omitted or an empty matrix is passed) Use the
%            Ledoit-Wolf algorithm of [1].
%          1 Use the empirical Bayesian estimator (the Hanff estimator),
%            which is also described in [1].
%
%OUTPUTS: SStar The pXp covariance matrix estimate.
%
%Note that the algorithms are not unbiased. However, the Ledoit-Wolf
%algorithm is shown in [1] to be asymptotically unbiased (as n gets large
%keeping p/n fixed).
%
%EXAMPLE:
%Here, we look at the performance when considering random diagonal
%covariance matrices whose elements are drawn from the log-normal
%distribution. We look at the percentage relative improvement in average
%loss (PRIAL) criterion, which is essentially the relative improvement in
%the squared Frobenius norm of the error of the sample covariance matrix
%versus the Ledoit-Wolf covariance matrix esitmate.
% numRuns=10000;
% p=20;
% n=40;
% %Parameters for the log-normal distribution from which the covariance
% %matrix is drawn.
% mu=1;
% Sigma=1;
% centerData=false;
%
% SSigmaFro2=0;
% SStarSigmaFro2=0;
% for curRun=1:numRuns
%     x=LogNormalD.rand(p,mu,Sigma);
%     X=diag(x);%The current covariance matrix.
%     SX=cholSemiDef(X,'lower');
%     x=SX*randn(p,n);%Zero mean.
%
%     S=(1/n)*(x*x');%Sample covariance matrix.
%     SStar=covMatEstShrinkage(x,centerData);
%
%     SSigmaFro2=SSigmaFro2+norm(S-Sigma,'fro')^2;
%     SStarSigmaFro2=SStarSigmaFro2+norm(SStar-Sigma,'fro')^2;
% end
% SSigmaFro2=SSigmaFro2/numRuns;
% SStarSigmaFro2=SStarSigmaFro2/numRuns;
% PRIAL=(SSigmaFro2-SStarSigmaFro2)/SSigmaFro2
%One will typically see that the PRIAL is positive, around 0.08, indicating
%that the Ledoit-Wolf covariance matrix esitmate is 8% better than the
%sample covariance estimate. It was empirically noticed that setting
%centerData=true will further improve the performance on this example even
%though the samples generated are already zero-mean.
%
%REFERENCES:
%[1] O. Ledoit and M. Wolf "A Well-conditioned estimator for large-
%    dimensional covariance matrices," Journal of Multivariate Analysis,
%    vol. 88, no. 2, pp. 365-411, 2004.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(algorithm))
    algorithm=0;
end

if(nargin<2||isempty(centerData))
    centerData=true;
end

if(centerData)
    x=bsxfun(@minus,x,mean(x,2));
end

p=size(x,1);
n=size(x,2);
%The sample covariance matrix.
S=(1/n)*(x*x');

switch(algorithm)
    case 0%The Ledoit-Wolf covariance estimator.
        %Lemma 3.2
        m=trace(S)/p;

        %Lemma 3.3
        d2=norm(S-diag(m),'fro')^2;

        %Lemma 3.4
        b2Bar=0;
        for k=1:n
            b2Bar=b2Bar+norm(x(:,k)*x(:,k)'-S,'fro')^2;
        end
        b2Bar=(1/n^2)*b2Bar;
        b2=min(b2Bar,d2);

        %Lemma 3.5
        a2=d2-b2;

        %Equation 14
        SStar=(b2/d2)*diag(m)+(a2/d2)*S;
    case 1%The empirical Bayesian estimator (Haff estimator).
        mEB=det(S)^(1/p);
        if(mEB<=0||~isreal(mEB))
            mEB=trace(S)/p;
        end

        SStar=((p*n-2*n-2)/(p*n^2))*diag(mEB)+(n/(n+1))*S;
    otherwise
        error('Unknown algorithm specified.')
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary