text:
classdef ControlMode < Simulink.IntEnumType

    enumeration
        None(0)
        Manual(1)
        Acro(2)
        Stabilize(3)
        ALTCTL(4)
        POSCTL(5)
        Offboard(6)
    end

    methods (Static)

        function defaultValue = getDefaultValue()
            % GETDEFAULTVALUE  Returns the default enumerated value.
            %   If this method is not defined, the first enumeration is used.
            defaultValue = ControlMode.None;
        end

        function dScope = getDataScope()
            % GETDATASCOPE  Specifies whether the data type definition should be imported from,
            %               or exported to, a header file during code generation.
            dScope = 'Auto';
        end

        function desc = getDescription()
            % GETDESCRIPTION  Returns a description of the enumeration.
            desc = 'enumeration of control mode';
        end

        function fileName = getHeaderFile()
            % GETHEADERFILE  Returns path to header file if non-empty.
            fileName = '';
        end

        function flag = addClassNameToEnumNames()
            % ADDCLASSNAMETOENUMNAMES  Indicate whether code generator applies the class name as a prefix
            %                          to the enumeration.
            flag = true;
        end

    end
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function HeadModel = in_headmodel_bst(HeadModelFile, ApplyOrient, varargin)
% IN_HEADMODEL_BST: Wrapper for in_bst_headmodel (added for backward compatibility)

% @=============================================================================
% 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, 2016

% Display warning
warning('Use in_bst_headmodel instead.');
% Just call in_bst_headmodel()
if (nargin > 2)
    HeadModel = in_bst_headmodel(HeadModelFile, ApplyOrient, varargin{:});
elseif (nargin == 2)
    HeadModel = in_bst_headmodel(HeadModelFile, ApplyOrient);
elseif (nargin == 1)
    HeadModel = in_bst_headmodel(HeadModelFile);
end

author: brainstorm-tools


text:
function handles = dtiShowInplaneImages(handles,xSliceAxes,xSliceRgb,ySliceAxes,ySliceRgb,zSliceAxes,zSliceRgb)
%   Display the images in the main Matlab mrDiffusion (dtiFiberUI) window.
%
%   dtiShowInplaneImages(handles,xSliceAxes,xSliceRgb, ...
%      ySliceAxes,ySliceRgb,zSliceAxes,zSliceRgb)
%
% See also: dtiGetCurSlices
%
% HISTORY:
%  Authors: Dougherty, Wandell
%  2006.11.09 RFD: we now save and reset the axis limits, in case the use
%  has zoomed in on an axis.
%
% Stanford VISTA Team

% Axial slice.  The z-dimension is the inferior/superior
axes(handles.z_cut);
% We don't want to reset the axis limit when the figure is first
% initialized, so we do a simple test for that initial condition.
if(isempty(get(handles.z_cut,'Children')))
    handles.z_cut_img = image(zSliceAxes(1,:), zSliceAxes(2,:), zSliceRgb);
    axis(handles.z_cut, 'equal', 'tight', 'xy');
else
    a = axis(handles.z_cut);
    handles.z_cut_img = image(zSliceAxes(1,:), zSliceAxes(2,:), zSliceRgb);
    axis(handles.z_cut, 'equal', 'tight', 'xy', a);
end
ylabel('Y (mm)'); xlabel('X (mm)');

% Coronal slice.
axes(handles.y_cut);
ySliceRgb = permute(ySliceRgb,[2,1,3]);
if(isempty(get(handles.y_cut,'Children')))
    handles.y_cut_img = image(ySliceAxes(1,:), ySliceAxes(2,:), ySliceRgb);
    axis(handles.y_cut, 'equal', 'tight', 'xy');
else
    a = axis(handles.y_cut);
    handles.y_cut_img = image(ySliceAxes(1,:), ySliceAxes(2,:), ySliceRgb);
    axis(handles.y_cut, 'equal', 'tight', 'xy', a);
end
ylabel('Z (mm)'); xlabel('X (mm)');

% Sagittal slice.
% We want the x-cut (sagittal view) to have the nose to the left
axes(handles.x_cut);
xSliceRgb = permute(xSliceRgb,[2,1,3]);
if(isempty(get(handles.x_cut,'Children')))
    handles.x_cut_img = image(xSliceAxes(1,:), xSliceAxes(2,:), xSliceRgb);
    axis(handles.x_cut, 'equal', 'tight', 'xy');
else
    a = axis(handles.x_cut);
    handles.x_cut_img = image(xSliceAxes(1,:), xSliceAxes(2,:), xSliceRgb);
    axis(handles.x_cut, 'equal', 'tight', 'xy', a);
end
ylabel('Z (mm)'); xlabel('Y (mm)');

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function X = SynthMMV(Y, HforVec, err, normfac)

% Solves problem of the form:
% min ||X||_2,1 s.t. ||Y-HX|| < err

% Inputs
% Y - n x r observation vectors
% H - n x N projection operator
% err - error tolerance
% normfac - maximum eigenvalue of projection operator

% Outputs
% X - N x r input vectors to be recovered

% Converting input matrix to operator
explicitA = ~(ischar(HforVec) || isa(HforVec, 'function_handle'));
if (explicitA)
    HOp = opMatrix(HforVec); clear HforVec
else
    HOp = HforVec;
end

r = size(Y,2); % number of observations/inputs

H = opMatWrapper(HOp, r); % wrapper for handling matrix inputs

alpha = 1.05*(normfac^2);

X = H(Y,2); % Initialize X
N = size(X,1); % length of each input vector

maxIter = 100; % Define the maximum number of iterations
tol = 1e-4; % tolerance level
decfac = 0.5; % decrease factor for lambda

lambdaInit = decfac*max(max((abs(X)))); lambda = lambdaInit;

while lambda > lambdaInit*tol
    iter = 0;
    while iter < maxIter
        iter = iter + 1;
        for i = 1:N
            D(i,:) = (1/norm(X(i,:))).*abs(X(i,:));
        end

        B = X + (1/alpha)*H(Y-H(X,1),2);

        Xvec = SoftTh(B(:),(lambda/alpha).*D(:));
        X = reshape(Xvec,[N r]);

    end
    if norm(Y-H(X,1),'fro') < err
        break;
    end
    lambda = decfac*lambda;
end

    function  z = SoftTh(s,thld)
        z = sign(s).*max(0,abs(s)-thld);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function P=maxpooling_s2(A)

[h,w]=size(A);

if mod(h,2)==1
    A=[A;zeros(1,w)];
    h=h+1;
end
if mod(w,2)==1
    A=[A,zeros(h,1)];
    w=w+1;
end

A1=A(1:2:h-1,1:2:w-1);
A2=A(1:2:h-1,2:2:w);
A3=A(2:2:h,1:2:w-1);
A4=A(2:2:h,2:2:w);

P1=max(A1,A2);
P2=max(A3,A4);
P=max(P1,P2);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [fid, tree, dir] = fiff_open(fname)
%
% [fid, tree, dir] = fiff_open(fname)
%
% Open a fif file and provide the directory of tags
%
% fid     the opened file id
% tree    tag directory organized into a tree
% dir     the sequential tag directory
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.7  2009/03/30 11:37:37  msh
%   Added copying of measurement info blocks from the original like in mne_browse_raw
%
%   Revision 1.6  2008/11/16 21:31:23  msh
%   Added mne_transform_coordinates and new coordinate frame definitions
%
%   Revision 1.5  2006/05/03 19:03:19  msh
%   Eliminated the use of cast function for Matlab 6.5 compatibility
%
%   Revision 1.4  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/17 15:01:34  msh
%   More small improvements.
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_open';
verbose=false;
fid = fopen(fname,'rb','ieee-be');

if (fid < 0)
    error(me,'Cannot open file %s', fname);
end;
%
%   Check that this looks like a fif file
%
tag = fiff_read_tag_info(fid);
if tag.kind ~= FIFF.FIFF_FILE_ID
    error(me,'file does not start with a file id tag');
end
if tag.type ~= FIFF.FIFFT_ID_STRUCT
    error(me,'file does not start with a file id tag');
end
if tag.size ~= 20
    error(me,'file does not start with a file id tag');
end
tag = fiff_read_tag(fid);
if tag.kind ~= FIFF.FIFF_DIR_POINTER
    error(me,'file does have a directory pointer');
end
if nargout == 1
    fseek(fid,0,'bof');
    return;
end
%
%   Read or create the directory tree
%
if verbose
    fprintf(1,'\tCreating tag directory for %s...',fname);
end
dirpos = double(tag.data);
if dirpos > 0
    tag = fiff_read_tag(fid,dirpos);
    dir = tag.data;
else
    k = 0;
    fseek(fid,0,'bof');
    dir = struct('kind',{},'type',{},'size',{},'pos',{});
    while tag.next >= 0
        pos = ftell(fid);
        tag = fiff_read_tag_info(fid);
        k = k + 1;
        dir(k).kind = tag.kind;
        dir(k).type = tag.type;
        dir(k).size = tag.size;
        dir(k).pos  = pos;
    end
end
%
%   Create the directory tree structure
%
tree = fiff_make_dir_tree(fid,dir);
if verbose
    fprintf(1,'[done]\n');
end
%
%   Back to the beginning
%
fseek(fid,0,'bof');
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [P,Pi] = computePandPi(Dir_alpha,Dir2d_alpha)

K = size(Dir2d_alpha,2);
if length(size(Dir2d_alpha))==3
    Q = length(size(Dir2d_alpha));
else
    Q = 1;
end
P = zeros(K,K,Q);
if Q==1, Pi = zeros(1,K);
else, Pi = zeros(K,Q);
end
for i = 1:Q
    for j = 1:K
        PsiSum = psi(sum(Dir2d_alpha(j,:,i)));
        for k = 1:K
            P(j,k,i) = exp(psi(Dir2d_alpha(j,k,i))-PsiSum);
        end
        P(j,:,i) = P(j,:,i) ./ sum(P(j,:,i));
    end
    if Q==1
        PsiSum = psi(sum(Dir_alpha));
        for k = 1:K
            Pi(k) = exp(psi(Dir_alpha(k))-PsiSum);
        end
        Pi = Pi ./ sum(Pi);
    else
        PsiSum = psi(sum(Dir_alpha(:,i)));
        for k = 1:K
            Pi(k,i) = exp(psi(Dir_alpha(k,i))-PsiSum);
        end
        Pi(:,i) = Pi(:,i) ./ sum(Pi(:,i));
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
% segmentiris - peforms automatic segmentation of the iris region
% from an eye image. Also isolates noise areas such as occluding
% eyelids and eyelashes.
%
% Usage:
% [circleiris, circlepupil, imagewithnoise] = segmentiris(image)
%
% Arguments:
%	eyeimage		- the input eye image
%
% Output:
%	circleiris	    - centre coordinates and radius
%			          of the detected iris boundary
%	circlepupil	    - centre coordinates and radius
%			          of the detected pupil boundary
%	imagewithnoise	- original eye image, but with
%			          location of noise marked with
%			          NaN values
%
% Author:
% Libor Masek
% masekl01@csse.uwa.edu.au
% School of Computer Science & Software Engineering
% The University of Western Australia
% November 2003

function [circleiris, circlepupil, imagewithnoise] = segmentiris(eyeimage)

% define range of pupil & iris radii

%CASIA
lpupilradius = 28;
upupilradius = 75;
lirisradius = 80;
uirisradius = 150;

%    %LIONS
%    lpupilradius = 32;
%    upupilradius = 85;
%    lirisradius = 145;
%    uirisradius = 169;

% define scaling factor to speed up Hough transform
scaling = 0.4;

reflecthres = 240;

% find the iris boundary by Daugman's intefrodifferential operaror

        [rowp, colp, rp] = SearchInnerBoundary(double(eyeimage));

        [row, col, r] = SearchOuterBoundary(double(eyeimage), rowp, colp, rp);
        circleiris = [row, col, r];
        rowd = double(row);
        cold = double(col);
        rd = double(r);

        irl = round(rowd-rd);
        iru = round(rowd+rd);
        icl = round(cold-rd);
        icu = round(cold+rd);

        imgsize = size(eyeimage);

        if irl < 1
            irl = 1;
        end

        if icl < 1
            icl = 1;
        end

        if iru > imgsize(1)
            iru = imgsize(1);
        end

        if icu > imgsize(2)
            icu = imgsize(2);
        end

        % to find the inner pupil, use just the region within the previously
        % detected iris boundary
        imagepupil = double(eyeimage( irl:iru,icl:icu));

        rowp = double(rowp)-irl;
        colp = double(colp)-icl;
        r = double(rp);

        row = double(irl) + rowp;
        col = double(icl) + colp;

        row = round(row);
        col = round(col);

        circlepupil = [row col r];

% set up array for recording noise regions
% noise pixels will have NaN values
imagewithnoise = double(eyeimage);

%find top eyelid
topeyelid = uint8(imagepupil(1:(rowp-r),:));
lines = findline(topeyelid);

if size(lines,1) > 0
    [xl yl] = linecoords(lines, size(topeyelid));
    yl = double(yl) + irl-1;
    xl = double(xl) + icl-1;

    yla = max(yl);

    y2 = 1:yla;

    ind3 = sub2ind(size(eyeimage),yl,xl);
    imagewithnoise(uint32(ind3)) = NaN;
    imagewithnoise(round(y2), round(xl)) = NaN;
end

%find bottom eyelid
bottomeyelid = uint8(imagepupil((rowp+r):size(imagepupil,1),:));
lines = findline(bottomeyelid);

if size(lines,1) > 0

    [xl yl] = linecoords(lines, size(bottomeyelid));
    yl = double(yl)+ irl+rowp+r-2;
    xl = double(xl) + icl-1;

    yla = min(yl);

    y2 = yla:size(eyeimage,1);

    ind4 = sub2ind(size(eyeimage),yl,xl);
    imagewithnoise(uint32(ind4)) = NaN;
    imagewithnoise(round(y2), round(xl)) = NaN;

end

% %For CASIA, eliminate eyelashes by thresholding
ref = eyeimage < 80;
coords = find(ref==1);
imagewithnoise(coords) = NaN;

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris


text:
function dgrav=model_grav(sys,rr,rs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%obtains the gravitational delay correction for the effect of general
%relativity (red shift) to the GPS signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
dgrav=0;

if norm(rr)<=0||norm(rs)<=0,return;end

rec_module=sqrt(rr(1)^2+rr(2)^2+rr(3)^2);
sat_module=sqrt(rs(1)^2+rs(2)^2+rs(3)^2);
distance  =norm(rr-rs);

switch sys
    case glc.SYS_GLO,MU=glc.MU_GLO;
    case glc.SYS_GAL,MU=glc.MU_GAL;
    case glc.SYS_BDS,MU=glc.MU_BDS;
    otherwise,       MU=glc.MU_GPS;
end

dgrav=2*MU/glc.CLIGHT^2*log((rec_module+sat_module+distance)/(rec_module+sat_module-distance));

return

repo: GINav


text:
function [inflMap, colXCoord, rowYCoord, mi] = getLSInfluenceMapFactorMovie(LS)
%"getLSInfluenceMap"
%   Gets an image of the influence generated by the beam described in LS.
%   Use getDICOMLeafPositions to generate LS.
%
%JRA&KZ 02/8/05
%
%Usage:
%   function inflMap = getLSInfluenceMap(LS);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Maximum precision of leaf position, in mm. Varian End and Side Accuracy 1.0 mm at
%isocenter. End and Side Repeatability 0.5 mm
precision = .5;

%Get x max, min and round to precision value.

%If X jaws dosn't exist in DICOM
if ~isfield(LS,'xLimits')
    xMax = ceil(max(vertcat(LS.xLeafPositions{:}),[],1) / precision) * precision;
    xMin = floor(min(vertcat(LS.xLeafPositions{:}),[],1) / precision) * precision;
    LS.xLimits{1}(1) = xMin;
    LS.xLimits{1}(2) = xMax;
end

xMax = ceil(max(vertcat(LS.xLimits{:}),[],1) / precision) * precision;
xMin = floor(min(vertcat(LS.xLimits{:}),[],1) / precision) * precision;
fieldSize.x = max(xMax) - min(xMin);
fieldLim.x  = [max(xMax) min(xMin)];

yMax = ceil(max(vertcat(LS.yLimits{:}),[],1) / precision) * precision;
yMin = floor(min(vertcat(LS.yLimits{:}),[],1) / precision) * precision;
fieldSize.y = max(yMax) - min(yMin);
fieldLim.y  = [max(yMax) min(yMin)];

yRes = precision;
nyElements = ceil(fieldSize.y/yRes);
xRes = precision;
nxElements = ceil(fieldSize.x/xRes);

inflMap=zeros(nyElements, nxElements);
colDividerXCoord = linspace(fieldLim.x(2), fieldLim.x(1), nxElements+1);
rowDividerYCoord = linspace(fieldLim.y(2), fieldLim.y(1), nyElements+1);

if isfield(LS, 'yLeafPositions')
    rowLeafPositions = round(interp1(rowDividerYCoord, 1:nyElements+1, LS.yLeafPositions,'linear', 'extrap'));
    rowLeafPositions = clip(rowLeafPositions, 1, nyElements+1, 'limits');
    leafBoundariesToKeep = [diff(rowLeafPositions)>0;true];
    rowLeafPositions = rowLeafPositions(leafBoundariesToKeep);
    leavesToKeep = leafBoundariesToKeep(1:end-1);
else
    LS.xLeafPositions{1} = [xMin xMax];
    LS.meterSetWeight = {1};
    rowLeafPositions = [1 nyElements+1];
    leavesToKeep = 1;
end

if length(LS.meterSetWeight) == 1
    doses = LS.meterSetWeight{:};
else
    doses = [0 diff([LS.meterSetWeight{:}])];
end

backupMap = inflMap;
%h = waitbar(0,['Generating Fluence Map From MLC Positions For Beam ',num2str(beamIndex)],'Name','Please wait...');

for i=1:length(LS.xLeafPositions)
    mapMovie = backupMap;
    inflMap = backupMap;
    nLeaves = length(LS.xLeafPositions{i})/2;

    if length(LS.xLimits) > 1
        jpL = LS.xLimits{i}(1);
        jpR = LS.xLimits{i}(2);
    else
        jpL = LS.xLimits{1}(1);
        jpR = LS.xLimits{1}(2);
    end

    lpL = LS.xLeafPositions{i}(1:nLeaves);
    lpR = LS.xLeafPositions{i}(nLeaves+1:end);
    lpLK = lpL(leavesToKeep);
    lpRK = lpR(leavesToKeep);

    MLCopeningSize(:,i) = lpRK - lpLK;

    lpLCols = interp1(colDividerXCoord, 1:nxElements+1, lpLK, 'linear', 'extrap');
    lpRCols = interp1(colDividerXCoord, 1:nxElements+1, lpRK, 'linear', 'extrap');

    %Column divider positions of jaws.
    jpLCol = interp1(colDividerXCoord, 1:nxElements+1, jpL, 'linear', 'extrap');
    jpRCol = interp1(colDividerXCoord, 1:nxElements+1, jpR, 'linear', 'extrap');

    jpLCol = round(jpLCol);
    jpRCol = round(jpRCol);

    lpLCols = clip(lpLCols, jpLCol, jpRCol, 'limits');
    lpRCols = clip(lpRCols, jpLCol, jpRCol, 'limits');

    lpLCols = round(lpLCols);
    lpRCols = round(lpRCols);

    A1 = 0.0013;
    A2 = 0.078;
    K = 1.5;
    LAMDA = 7.69;

    for j=1:length(lpLCols)
        %HCF from output ratio for MLC fields Zhu, MedPhys
        F_X = abs(lpLCols(j) - (lpRCols(j)-1))*precision/10;
        F_Y = abs(rowLeafPositions(j+1) - rowLeafPositions(j))*precision/10;
        F_X = abs(lpLCols(j) - (lpRCols(j)-1))*precision/10;
        F_Y = abs(rowLeafPositions(j+1) - rowLeafPositions(j))*precision/10;
        F = (1+K)*F_X * F_Y/(K*F_X + F_Y);
        HCF = (1+A1*F)*(1+A2*(erf(F/LAMDA))^2)/((1+A1*10)*(1+A2*(erf(10/LAMDA))^2));

        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpLCols(j):lpRCols(j)-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpLCols(j):lpRCols(j)-1) + HCF*doses(i);
        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, jpLCol:lpLCols(j)-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, jpLCol:lpLCols(j)-1);
        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpRCols(j):jpRCol-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpRCols(j):jpRCol-1);
    end

    %waitbar(i/length(LS.xLeafPositions));
      % frame = inflMap;
       %imagesc(inflMap);
       %mi(:,:,i) = inflMap;
       mapMovie = inflMap;
       mapMovie(mapMovie == 0) = 1;
       mapMovie(mapMovie ~= 1) = 2;
       %colormap([0 0 0; 1 1 1]);
       %mi(:,:,i) = inflMap;
       mi(i) = im2frame(mapMovie, [0 0 1; 1 0 0]);
%        drawnow;
%        pause(.006);
end
%close(h);
colXCoord = colDividerXCoord(1:end-1) + precision/2;
rowYCoord = rowDividerYCoord(1:end-1) + precision/2;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function showmesh3(node,elem,expr,varargin)
%% SHOWMESH3 displays a tetrahedron mesh in 3-D.
%
%    showmesh3(node,elem) displays a 3-dimensional tetrahderon mesh given
%    by node and elem matrices; see meshdoc for the data structure:
%    node and elem.
%
%    showmesh3(node,elem,expr) displays parts of the mesh specificed by the
%    expression. For example, showmesh3(node,elem,'~(x>=0 & y>=0)') only
%    shows the tetrahedron not in the first quadrant.
%
%    showmesh3(node,elem,expr,'param','value','param','value'...) allows
%    additional patch param/value pairs to be used when displaying the
%    mesh. For example, the default transparency parameter is set to 0.5.
%    You can overwrite this value by using the param pair ('FaceAlpha',
%    value). The value has to be a number between 0 and 1. Other parameters
%    include: 'Facecolor', 'Edgecolor' etc.
%
%    For meshes with large data, the 3-D graphics is very slow. You may use
%    showboundary3 to display the boundary surface mesh only.
%
%   Example:
%     % A mesh for a cube
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     [node,elem] = uniformbisect3(node,elem);
%     subplot(1,2,1);
%     showmesh3(node,elem); pause(1)
%     subplot(1,2,2);
%     showmesh3(node,elem,'~(x>=0 & y>=0 & z>=0)','FaceAlpha',0.25);
%     axis on; view([59,20])
%
%   See also showboundary3, showsolution3, showmesh.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin >= 3) && (any(expr))
    x = node(:,1);  y = node(:,2);  z = node(:,3); %#ok<NASGU>
    incl = find(eval(expr));
    elem = elem(any(ismember(elem,incl),2),:);
end
h = tetramesh(elem(:,1:4),node,ones(size(elem,1),1));
if nargin > 3
    set(h,varargin{1:end})
else % default display properties
    set(h,'FaceAlpha',0.4);
end
view(3);
axis off; axis equal; axis tight

repo: Aero-matlab


text:
function view_qva(lab1,valueMap)
    % view_maxz plots the maxz LTA values calculated
    % with maxzlta.m or other similar values as a color map
    % needs valueMap, gx, gy, stri
    %
    % define size of the plot etc.
    %
%     The Z-Value Map Window
%
%
%           This window displays seismicity rate changes
%           as z-values using a color code. Negative
%           z-values indicate an increase in the seismicity
%           rate, positive values a decrease.
%           Some of the menu-bar options are
%           described below:
%
%           Threshold: You can set the maximum size that
%             a volume is allowed to have in order to be
%             displayed in the map. Therefore, areas with
%             a low seismicity rate are not displayed.
%             edit the size (in km) and click the mouse
%             outside the edit window.
%         'FixAx: You can chose the minimum and maximum
%                  values of the color-legend used.
%         'Polygon: You can select earthquakes in a
%           polygon either by entering the coordinates or
%           defining the corners with the mouse
%
%         'Circle: Select earthquakes in a circular volume:
%                Ni, the number of selected earthquakes can
%                be edited in the upper right corner of the
%                window.
%           Refresh Window: Redraws the figure, erases
%                 selected events.
%
%           zoom: Selecting Axis -> zoom on allows you to
%                 zoom into a region. Click and drag with
%                 the left mouse button. type <help zoom>
%                 for details.
%           Aspect: select one of the aspect ratio options
%           Text: You can select text items by clicking. The
%                 selected text can be rotated, moved, you
%                 can change the font size etc.
%                 Double click on text allows editing it.

    report_this_filefun();
    ZG=ZmapGlobal.Data;
    % Find out if figure already exists
    %
    qmap=findobj('Type','Figure','-and','Name','q-detect-map');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(qmap)
        qmap = figure_w_normalized_uicontrolunits( ...
            'Name','q-detect-map',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        lab1 = 'day/night ratio';
        create_my_menu();

        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Info ',...
            'callback',@callbackfun_001);

        uicontrol('Units','normal',...
            'Position',[.92 .80 .08 .05],'String','set ni',...
            'callback',@callbackfun_006)

        set_nia = uicontrol('style','edit','value',ni,'string',num2str(ni));
        set(set_nia, 'Callback', @callbackfun_007);
        set(set_nia,'units','norm','pos',[.94 .85 .06 .05],'min',10,'max',10000);
        nilabel = uicontrol('style','text','units','norm','pos',[.90 .85 .04 .05]);
        set(nilabel,'string','ni:','background',[.7 .7 .7]);

        ZG.tresh_km = nan; re4 = valueMap;
        nilabel2 = uicontrol('style','text','units','norm','pos',[.60 .92 .25 .06]);
        set(nilabel2,'string','MinRad (in km):','background',color_fbg);
        set_ni2 = uicontrol('style','edit','value',ZG.tresh_km,'string',num2str(ZG.tresh_km),...
            'background','y');
        set(set_ni2, 'Callback', @callbackfun_008);
        set(set_ni2,'units','norm','pos',[.85 .92 .08 .06],'min',0.01,'max',10000);

        uicontrol('Units','normal',...
            'Position',[.95 .93 .05 .05],'String','Go ',...
            'callback',@callbackfun_009)

        colormap(cool)

    end   % This is the end of the figure setup

    % Now lets plot the color-map of the z-value
    %
    figure(qmap);
    delete(findobj(qmap,'Type','axes'));
    % delete(sizmap);
    reset(gca)
    cla
    set(gca,'NextPlot','replace')
    watchon;
    set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','SortMethod','childorder')

    rect = [0.18,  0.10, 0.7, 0.75];
    rect1 = rect;

    % find max and min of data for automatic scaling
    %
    ZG.maxc = max(valueMap(:));
    ZG.maxc = fix(ZG.maxc)+1;
    ZG.minc = min(valueMap(:));
    ZG.minc = fix(ZG.minc)-1;

    % set values gretaer ZG.tresh_km = nan
    %
    re4 = valueMap;
    re4(r > ZG.tresh_km) = nan;

    % plot image
    %
    orient landscape
    %set(gcf,'PaperPosition', [0.5 1 9.0 4.0])

    axes('position',rect)
    set(gca,'NextPlot','add')
    pco1 = pcolor(gx,gy,re4);

    axis([ min(gx) max(gx) min(gy) max(gy)])
    axis image
    set(gca,'NextPlot','add')

    shading(ZG.shading_style);

    fix_caxis.ApplyIfFrozen(gca);

    title([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
        'Color','r','FontWeight','bold')

    xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

    % plot overlay
    %
    set(gca,'NextPlot','add')
    zmap_update_displays();
    globalcatalog = ZG.primeCatalog;
    ploeq = plot(globalcatalog.X,globalcatalog.Y,'k.');
    set(ploeq,'Tag','eq_plot','MarkerSize',ZG.ms6,'Marker',ty,'Color',ZG.someColor,'Visible','on')

    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    h1 = gca;
    hzma = gca;

    % Create a colorbar
    %
    h5 = colorbar('horiz');
    set(h5,'Pos',[0.35 0.05 0.4 0.02],...
        'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

    rect = [0.00,  0.0, 1 1];
    axes('position',rect)
    axis('off')
    %  Text Object Creation
    txt1 = text(...
        'Units','normalized',...
        'Position',[ 0.33 0.07 0 ],...
        'HorizontalAlignment','right',...
        'FontSize',ZmapGlobal.Data.fontsz.s,....
        'FontWeight','bold',...
        'String',lab1);

    % Make the figure visible
    %
    set(gca,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    figure(qmap);
    axes(h1)
    watchoff(qmap)

    %% ui functions
    function create_my_menu()
        add_menu_divider();

        add_symbol_menu('eq_plot');

        options = uimenu('Label',' Select ');
        uimenu(options,'Label','Refresh ','MenuSelectedFcn',@callbackfun_002)
        uimenu(options,'Label','Select EQ in Circle','MenuSelectedFcn',@callbackfun_003)
        uimenu(options,'Label','Select EQ in Polygon ','MenuSelectedFcn',@callbackfun_004)

        op1 = uimenu('Label',' Maps ');
        uimenu(op1,'Label','day/night value map',...
            'MenuSelectedFcn',@callbackfun_005)

        add_display_menu(1);
    end

    %% callback functions

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        web(['file:' hodi '/help/quarry.htm']) ;
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        delete(findobj(qmap,'Type','axes'));
        view_qva(lab1,valueMap);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        circle;
        watchoff(qmap);
        hisgra(ZG.newt2.Date.Hour,'Hour',ZG.newt2.Name);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        stri = 'Polygon';
        h1 = gca;
        cufi = gcf;
        selectp;
        hisgra(ZG.newt2,'Hour');
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1 ='day/night ratio';
        valueMap = old;
        view_qva(lab1,valueMap);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2num(set_nia.String);
        'String';
        num2str(ni);
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        %FIXME callback does nothing
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ZG.tresh_km=str2double(set_ni2.String);
        set_ni2.String=num2str(ZG.tresh_km);
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ;
        pause(1);
        re4 =valueMap;
        view_bva(lab1,valueMap);
    end
end

author: CelsoReyes


text:
function xf = cfft2(x)

% calculate output size
in_sz = size(x);

% if both dimensions are odd
if all(mod(in_sz(1:2), 2) == 1)
    xf = fftshift(fftshift(fft2(x), 1), 2);
else
    out_sz = in_sz;
    out_sz(1:2) = out_sz(1:2) + mod(out_sz(1:2)+1,2);

    % allocate
    xf = complex(zeros(out_sz, 'single'));

    xf(1:in_sz(1),1:in_sz(2),:,:) = fftshift(fftshift(fft2(x), 1), 2);

    if out_sz(1) ~= in_sz(1)
        xf(end,:,:,:) = conj(fliplr(xf(1,:,:,:)));
    end
    if out_sz(2) ~= in_sz(2)
        xf(:,end,:,:) = conj(flipud(xf(:,1,:,:)));
    end
end

author: martin-danelljan

repo: Continuous-ConvOp


text:
%kkextract 'Extract Rectangular Region from Object'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kextract.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input object'
% Integer: wsize 'Width ', default: 1: 'Specifies the region width'
% Integer: hsize 'Height ', default: 1: 'Specifies the region height'
% Integer: dsize 'Depth ', default: 1: 'Specifies the region depth'
% Integer: tsize 'Time ', default: 1: 'Specifies the region size on time dimension'
% Integer: esize 'Elements', default: 1: 'Specifies the region size on elements dimension'
% Integer: woff 'Width ', default: 0: 'Specifies origin of region on width dimension'
% Integer: hoff 'Height ', default: 0: 'Specifies origin of region on height dimension'
% Integer: doff 'Depth ', default: 0: 'Specifies origin of region on depth dimension'
% Integer: toff 'Time ', default: 0: 'Specifies origin of region on time dimension'
% Integer: eoff 'Elements', default: 0: 'Specifies origin of region on elements dimension'
% OutputFile: o 'Output', required: 'Resulting output region'
%
% Example: o = kkextract(i, {'i','';'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1;'woff',0;'hoff',0;'doff',0;'toff',0;'eoff',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kextract - Extract Rectangular Region from Object
%
%  DESCRIPTION
%
% The "Extract Region" operator (kextract) extracts a region of data
% from the source object "Input" (i) and places it in the destination
% object "Output" (o).  The region size is specified by the Width (wsize),
% Height (hsize), Depth (dsize), Time (tsize), and Elements (esize) Region
% Size parameters.  The region will be extracted starting at the point defined
% by the Region Origin parameters, Width (woff), Height (hoff), Depth (doff),
% Time (toff), and Elements (eoff).  The datum at this starting point will be
% placed at the origin (0,0,0,0,0) in the destination object.
%
% The output data object size for each dimension is determined first by the
% Region Size parameters "if" they are specified, and then by the Region
% Origin parameters in combination with the original source data size.  In other
% words, if the size for a dimension is not selected on the graphical user
% interface (or explicitly specified when running kextract from the command
% line), all data on that dimension will be extracted, starting at the specified
% origin.   This will ensure the following behavior.  Given an input data file
% with dimensions (width=100, height=100, depth=10, time=1, elements = 1),
% and the following "specified" parameters:
%
% 	Region Size:  Width = 50, Height = 50
% 	Region Origin:  Width = 2, Height = 2, Depth = 2
%
% the output object size will be (width=50, height=50, depth=8, time=1,
% elements = 1).
%
% The extract operation is based on implicit indexing.  This means that if
% location or time data exist, the extract operation is not done in terms
% an interpretation of the location/time data values, but in terms of the
% implicit indexing of these data (which is specified by the width, height,
% depth, time, elements indices of the polymorphic data model).  This program
% can operate on any combination of value, mask, location and/or time data
% in the input object.  Map data is unaffected by the extract operation.
%
% If the "Region Origin" and/or the \fIRegion Size\fP fall outside of
% the size of the input object, destination object value, location, and time
% data will be padded with zero, and mask data will be padded with one.
%
% The subobject position attribute of a data object defines the coordinates
% at which that subobject is located in its parent data object.
% The "Record Region Origin" option (subpos) allows the extract operator to
% set the sub-object position attribute in the output object to the Region
% origin.  If this option is set to "Yes", the subobject position of the
% output object will be set to (woff, hoff, doff, toff, eoff).  If "Record
% Region Origin" is set to "No", the output subobject position is set to
% (0,0,0,0,0).  Programs like the insert operator (kinsert) can use this
% attribute for automatically inserting the subregion back into the original
% after processing.
%
% The subobject position attribute can also be explicitly set by using the
% Set Data Attributes operator (ksetdattr), and can be printed using the
% Print File Information operator (kfileinfo).
%
% The destination object data types are the same as the source object
% data types.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% This program has not yet been updated to completely support the
% polymorphic data model.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkextract(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkextract(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'wsize', 1;'hsize', 1;'dsize', 1;'tsize', 1;'esize', 1;'woff', 0;'hoff', 0;'doff', 0;'toff', 0;'eoff', 0;'o', '__output'};
maxval={0,2,2,2,2,2,0,0,0,0,0,0};
minval={0,2,2,2,2,2,0,0,0,0,0,0};
istoggle=[0,1,1,1,1,1,1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','Integer','Integer','Integer','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kextract"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function start_points = perform_lloyd_mesh(vertex,faces, start_points, options)

% perform_lloyd_mesh - perform lloyd relaxation to sample point on a mesh
%
%   start_points = perform_lloyd_mesh(vertex,faces, start_points, options);
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;

if size(vertex,1)>size(vertex,2)
    vertex = vertex';
end
if size(faces,1)>size(faces,2)
    faces = faces';
end

niter_lloyd = getoptions(options, 'niter_lloyd', 1);
if niter_lloyd>1
    for i=1:niter_lloyd
        start_points = perform_lloyd_mesh(vertex,faces, start_points, options)
        options.lambda = []; % enfore recomputing
    end
    return;
end

lambda = getoptions(options, 'lambda', []);
edges_id = getoptions(options, 'edges_id', []);
Q = getoptions(options, 'Q', []);
if isempty(lambda) || isempty(edges_id) || isempty(Q)
    % update voronoi
    [Q,DQ, ve, edges_id, lambda] = compute_voronoi_mesh(vertex,faces, start_points, options);
end

% compute distances for start
ne = length(lambda);
n = size(vertex,2);
% compute edge length
d = diff( reshape(vertex(:,edges_id'), [3 ne 2]), 1, 3 );
d = sqrt( sum(d.^2,1) );
d1 = (1-lambda).*d;
d2 = lambda.*d;
% compute initial values for FM
values = zeros(n,1); cnt = zeros(n,1);
for i=1:ne
    values(edges_id(1,i)) = values(edges_id(1,i)) + d1(i);
    cnt(edges_id(1,i)) = cnt(edges_id(1,i)) + 1;
    values(edges_id(2,i)) = values(edges_id(2,i)) + d1(i);
    cnt(edges_id(2,i)) = cnt(edges_id(2,i)) + 1;
end
values = values./cnt;
% perform FM from edge points
pts = unique( edges_id(:) );
options.values = values(pts);
options.values = [];
D = perform_fast_marching_mesh(vertex, faces, pts, options);

% seed new locations
nstart = length(start_points);
for i=1:nstart
    I = find(Q(:,1)==i);
    [v,k] = max(D(I));
    start_points(i) = I(k);
end
start_points = start_points(:);

author: alecjacobson

repo: gptoolbox


text:
function vol = tetrahedronVolume(vertices, varargin)
%TETRAHEDRONVOLUME Signed volume of a tetrahedron
%
%   VOL = tetrahedronVolume(TETRA)
%   Comptues the siged volume of the tetrahedron TETRA defined by a 4-by-4
%   array representing the polyhedron vertices.
%
%   Example
%     vi = [0 0 0;1 0 0;0 1 0;0 0 1];
%     tetrahedronVolume(vi)
%     ans =
%         0.1667
%
%     [V F] = createTetrahedron;
%     tetrahedronVolume(V)
%     ans =
%         -.3333
%
%   See also
%   meshes3d, createTetrahedron, meshVolume
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-04-05,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

if nargin == 2
    tetras = varargin{1};
    nTetras = size(tetras, 1);
    vol = zeros(nTetras, 1);
    for i = 1:nTetras
        tetra = tetras(i,:);
        vol(i) = det(bsxfun(@minus, vertices(tetra(2:4),:), vertices(tetra(1),:))) / 6;
    end
    return;
end

% control on inputs
if nargin == 4
    vertices = [vertices ; varargin{1} ; varargin{2} ; varargin{3}];
end

if size(vertices, 1) < 4
    error('Input vertex array requires at least 4 vertices');
end

% compute volume of tetrahedron, using first vertex as origin
vol = det(vertices(2:4,:) - vertices([1 1 1],:)) / 6;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [data, headerText] = load_ascii(filename,delimiter,header,nlines,offset)

% This function loads data from a tab delimited or csv ascii file similar
% to dlmread/importdata. If the file has fixed width columns and an offset
% is desired, it uses this information to quickly scan to desired area
% (making it significantly faster in those scenarios). Useful when working
% with very large ascii files that cannot be entirely loaded into memory at
% one time. If a file is loaded in with no offset, load_ascii is comparable
% to dlmread (faster/slower depends on delimiter). Importdata is many times
% slower than both.
%
% Inputs:
%
% filename: absolute or relative file path
% header: number of lines for header (defaults to 0)
% nlines: number of lines to read in (defaults to inf)
% offset: number of lines to skip (defaults to 0) (assumes data is
% written with a fixed field width. If this is not the
% case, offset option will not work.
% delimiter: delimiter in ascii file (defaults to tab)
%
% outputs:
% data: matrix of data in file
% headerText: text of header as specified by header input
%
% Example usage:
% Consider scenario where we have a very large ascii file (fixed width tab
% delimited, 24 header lines). We are only interested in a section of 2e6
% lines, 20e6 past the end of the header.
% %syntax
% %[data,headerText] = load_ascii(filename,delimiter,header,nlines,offset)
% tic;data = load_ascii(filename,'\t',24,2e6,20e6);toc;
% tic;data2 = dlmread(filename,'\t',[20e6+24 0 12e6+23 2]);toc;
% Elapsed time is 3.838641 seconds.
% Elapsed time is 21.537717 seconds.

% initialize some things
if nargin<3, header = 0;end
if nargin<4, nlines = inf;end
if nargin<5, offset = 0;end
if nargin<2, delimiter = '\t';end

fid = fopen(filename);
headerCell = cell(header,1);
for i=1:header
    headerCell{i} = fgets(fid);
end
headerText = char(headerCell);

% First line of Data...
fline = fgets(fid);

% determine number of columns
delKey = double(sprintf(delimiter));
ncols = sum(fline == delKey)+1;

% determine number of bytes per line
bytesPerLine = length(fline);

% skip offset
fseek(fid,bytesPerLine*(offset-1),0);

% read data and keep native format
d = fread(fid,bytesPerLine*nlines,'char=>char');

%convert format to matrix
format = ['%f' delimiter];
if ~strcmp(delimiter,'\t'),format = repmat(['%f' delimiter],1,ncols);format = format(1:end-1);end
data = sscanf(d,format,[ncols,nlines])';

fclose(fid);
end


text:
function info = cnn_imagenet_evaluate(varargin)
% CNN_IMAGENET_EVALUATE   Evauate MatConvNet models on ImageNet

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.dataDir = fullfile('data', 'ILSVRC2012') ;
opts.expDir = fullfile('data', 'imagenet12-eval-vgg-f') ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.modelPath = fullfile('data', 'models', 'imagenet-vgg-f.mat') ;
opts.networkType = [] ;
opts.lite = false ;
opts.numFetchThreads = 12 ;
opts.train.batchSize = 128 ;
opts.train.numEpochs = 1 ;
opts.train.gpus = [] ;
opts.train.prefetch = true ;
opts.train.expDir = opts.expDir ;

opts = vl_argparse(opts, varargin) ;
display(opts);

% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------

if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

net = load(opts.modelPath) ;
if isfield(net, 'net') ;
  net = net.net ;
end
% Cannot use isa('dagnn.DagNN') because it is not an object yet
isDag = isfield(net, 'params') ;

if isDag
  opts.networkType = 'dagnn' ;
  net = dagnn.DagNN.loadobj(net) ;
  trainfn = @cnn_train_dag ;

  % Drop existing loss layers
  drop = arrayfun(@(x) isa(x.block,'dagnn.Loss'), net.layers) ;
  for n = {net.layers(drop).name}
    net.removeLayer(n) ;
  end

  % Extract raw predictions from softmax
  sftmx = arrayfun(@(x) isa(x.block,'dagnn.SoftMax'), net.layers) ;
  predVar = 'prediction' ;
  for n = {net.layers(sftmx).name}
    % check if output
    l = net.getLayerIndex(n) ;
    v = net.getVarIndex(net.layers(l).outputs{1}) ;
    if net.vars(v).fanout == 0
      % remove this layer and update prediction variable
      predVar = net.layers(l).inputs{1} ;
      net.removeLayer(n) ;
    end
  end

  % Add custom objective and loss layers on top of raw predictions
  net.addLayer('objective', dagnn.Loss('loss', 'softmaxlog'), ...
               {predVar,'label'}, 'objective') ;
  net.addLayer('top1err', dagnn.Loss('loss', 'classerror'), ...
               {predVar,'label'}, 'top1err') ;
  net.addLayer('top5err', dagnn.Loss('loss', 'topkerror', ...
                                     'opts', {'topK',5}), ...
               {predVar,'label'}, 'top5err') ;

  % Make sure that the input is called 'input'
  v = net.getVarIndex('data') ;
  if ~isnan(v)
    net.renameVar('data', 'input') ;
  end

  % Swtich to test mode
  net.mode = 'test' ;
else
  opts.networkType = 'simplenn' ;
  net = vl_simplenn_tidy(net) ;
  trainfn = @cnn_train ;
  net.layers{end}.type = 'softmaxloss' ; % softmax -> softmaxloss
end

% Synchronize label indexes used in IMDB with the ones used in NET
imdb = cnn_imagenet_sync_labels(imdb, net);

% Run evaluation
[net, info] = trainfn(net, imdb, getBatchFn(opts, net.meta), ...
                      opts.train, ...
                      'train', NaN, ...
                      'val', find(imdb.images.set==2)) ;

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;
bopts = meta.normalization ;
bopts.numThreads = opts.numFetchThreads ;

% Most networks are trained by resizing images to 256 pixels and then
% cropping a slightly smaller subarea. Reproduce this effect (center
% crop) for a more accurate evaluation (it also avoids resizing the
% images if these have been pre-processed to be of this size, which
% accelerates everything).

bopts.border = 256 - meta.normalization.imageSize ;

switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getSimpleNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
im = cnn_imagenet_get_batch(images, opts, ...
                            'prefetch', nargout == 0) ;
labels = imdb.images.label(batch) ;

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
im = cnn_imagenet_get_batch(images, opts, ...
                            'prefetch', nargout == 0) ;
if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  inputs = {'input', im, 'label', imdb.images.label(batch)} ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
classdef TestTonemapDrago
    %TestTonemapDrago

    methods (Static)
        function test_1
            hdr = cv.imread(fullfile(mexopencv.root(),'test','memorial.hdr'), 'Flags',-1);
            hdr = cv.resize(hdr, 0.5, 0.5);

            obj = cv.TonemapDrago();
            obj.Gamma = 2.2;
            obj.Saturation = 1.0;
            obj.Bias = 0.85;

            ldr = obj.process(hdr);
            validateattributes(ldr, {'single'}, ...
                {'size',size(hdr), 'nonnegative'});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [x_struct, z_struct, u_struct, p_struct, ...
          x_bounds, z_bounds, u_bounds, p_bounds, ...
          x_order] = vars(varsfh, userdata)

vh = ocl.VarHandler(userdata);
varsfh(vh);

x_struct = vh.x_struct;
z_struct = vh.z_struct;
u_struct = vh.u_struct;
p_struct = vh.p_struct;

x_bounds = vh.x_bounds;
z_bounds = vh.z_bounds;
u_bounds = vh.u_bounds;
p_bounds = vh.p_bounds;

x_order = vh.x_order;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
%Tristan Ursell
%Relative Noise Transform
%(c) November 2012
%
%Iout=relnoise(Iin,sz,sigma);
%Iout=relnoise(Iin,sz,sigma,'field');
%[Iout,Ivar]=relnoise(Iin,sz,sigma,...);
%[Iout,Ivar,Imean]=relnoise(Iin,sz,sigma,...);
%
%Iin = the input image, of any numerical class.
%
%sz = (3 < sz < min(size(Iin))) is the size of the filter block used to
%calculate means and variances.  This value must be odd.
%
%sigma (sigma > 0) is the weighting parameter that defines the standard
%deviation relative to the filter block's standard deviation around which
%the center pixel will be Gaussian weighted. Setting sigma = 1 weights the
%current pixel using the STD of the current filter block. Lower values
%bring the current pixel closer to the mean, while high values are more
%tolerant of variations.  As sigma -> Inf, Iout = Iin.
%
%The field 'plot' will create an output plot comparing this transform to
%the original image, a Gaussian blur with STD = sz/2, and median filter
%with block size equal to sz.  At first glance, this filter appears similar
%to a median transform, but it does a better job of preserving local
%intensity extrema.  Comparison with the median filter requires the Image
%Processing Toolbox, but the rest of the script does not.
%
%The field 'disk' or 'square' will choose between using a disk or square
%filter block shape, where sz is the disk diameter or square side length.
%The default is square.
%
%The field 'custom' may be followed by a user-defined logical matrix or
%strel, e.g. relnoise(Iin,sz,sigma,'custom',strel('line',30,45)).  In this
%case 'sz' will be unused.
%
%Iout is the transformed output image.
%
%Ivar is the variance of the pixel intensities in the filter block at every
%point in the image -- essentially the spatially varying variance of the
%image.
%
%Imean is the mean smoothed image using the filter block, equivalent to a
%convolution averaging filter with the specified neighborhood.
%
%see also: wiener2  filter2
%
%Example:
%
%Iin=imread('spot_test.tif');
%
%Iout=relnoise(Iin,3,0.5,'square','plot');
% %OR
%Iout=relnoise(Iin,3,0.5,'custom',strel('line',30,0),'plot');
%
%figure;
%subplot(1,2,1)
%imagesc(Iout-double(Iin))
%title('What was removed from the original image.')
%axis equal tight
%box on
%
%subplot(1,2,2)
%imagesc(abs(fftshift(fft(Iout-double(Iin)))))
%title('FFT of difference between original and filtered images.')
%axis equal tight
%box on
%

function [varargout]=relnoise(Iin,sz,sigma,varargin)

%Nans
Inan=isnan(Iin);
if sum(Inan(:))>0
    error('Input matrix contains NaNs.')
end

%convert type
Iin=double(Iin);

%check filter size
if or(sz<1,sz>min(size(Iin)))
    error('The filter size is out of bounds.')
end

if mod(sz,2)~=1
    error('The filter size must be an odd integer.')
end

%parse field input
f1=find(strcmp('plot',varargin),1);
f2=find(strcmp('disk',varargin),1);
f3=find(strcmp('custom',varargin),1);

%choose plot option
if ~isempty(f1)
    plotq=1;
else
    plotq=0;
end

%choose filter type
if ~isempty(f3)
    hood_temp=varargin{f3+1};
    if strcmp(class(hood_temp),'strel')
        hood=hood_temp.getnhood;
    else
        hood=hood_temp;
    end

    sz=round(mean(size(hood)));
elseif ~isempty(f2)
    %disk filter block
    Xdisk = ones(sz,1)*(-(sz-1)/2:(sz-1)/2);
    Ydisk = (-(sz-1)/2:(sz-1)/2)'*ones(1,sz);
    Zdisk = sqrt(Xdisk.^2 + Ydisk.^2);

    hood=zeros(sz,sz);
    hood(Zdisk<=(sz-1)/2)=1;
else
    %square filter block
    hood=ones(sz,sz);
end

%convert hood class
hood=single(hood);

%calcualte means and variances
hood_sz=sum(hood(:));

%perform convolultion and normalization
Imean0=conv2(Iin,hood,'same')/hood_sz;
Inorm=conv2(ones(size(Iin)),hood,'same')/hood_sz;
Imean=Imean0./Inorm;
Ivar=conv2(Iin.^2,hood,'same')./Inorm*1/hood_sz-Imean.^2;

%compute weight matrix
W=exp(-(Iin-Imean).^2./(2*sigma^2*Ivar));

%correct for zero variance pixels
W(Ivar==0)=0;

%compute output image
if sigma==0
    Iout=Imean;
else
    Iout=Iin.*W+(1-W).*Imean;
end

%handle outputs
if nargout==1
    varargout{1}=Iout;
elseif nargout==2
    varargout{1}=Iout;
    varargout{2}=Ivar;
elseif nargout==3
    varargout{1}=Iout;
    varargout{2}=Ivar;
    varargout{3}=Imean;
elseif nargout==0
else
    error('Incorrect number of output arguments.')
end

%plot comparisons
if plotq==1

    figure;
    subplot(2,2,1)
    imagesc(Iin)
    xlabel('X')
    ylabel('Y')
    box on
    axis equal tight
    title('Original Image')

    subplot(2,2,2)
    imagesc(Iout)
    xlabel('X')
    ylabel('Y')
    box on
    axis equal tight
    title(['Relative Noise Reduction (this filter), size = ' num2str(sz) ', sigma = ' num2str(sigma)])

    subplot(2,2,3)
    %look for image processing toolbox
    boxes=ver;
    gotit=strfind([boxes.Name],'Image Processing Toolbox');
    if ~isempty(gotit)
        Iout2=medfilt2(Iin,[sz,sz],'symmetric');
    else
        disp('Sorry, you do not have the Image Processing Toolbox.')
        disp('The medfilt2 comparison image cannot be generated.')
        disp('Disabling `plot` will stop this message.')
        Iout2=Imean;
    end
    imagesc(Iout2)
    xlabel('X')
    ylabel('Y')
    box on
    axis equal tight

    if ~isempty(gotit)
        title(['Median Filter of size ' num2str(sz)])
    else
        title(['Mean Filter of size ' num2str(sz)])
    end

    %construct Gaussian filter (without Image Processing Toolbox)
    sz2=round(3/2*sz);
    Xgauss = ones(sz2,1)*(-(sz2-1)/2:(sz2-1)/2);
    Ygauss = (-(sz2-1)/2:(sz2-1)/2)'*ones(1,sz2);
    Zgauss = exp(-(Xgauss.^2+Ygauss.^2)/(2*(sz/2)^2));
    Zgauss = Zgauss/sum(Zgauss(:));
    Iout3=conv2(Iin,Zgauss,'same');

    subplot(2,2,4)
    imagesc(Iout3)
    xlabel('X')
    ylabel('Y')
    box on
    axis equal tight
    title(['Gaussian Blur with STD = ' num2str(sz/2)])
    pause(0.1)
end


text:
function roi = dtiRoiErode(roi, cropAxes)
%
% roi = dtiRoiErode(roi, [cropAxes=[1 1 1]])
%
% Erodes the roi by n points along any of the three axes. E.g.,
% cropAxes = [0 1 1] will erode by 1 point along Y and Z,
% cropAxes = [3 2 2] will erode by 3 points along X, and 2 points
% along Y and Z.
%
% HISTORY:
% 2009.08.19 RFD wrote it.

if(~exist('cropAxes','var')||isempty(cropAxes))
    cropAxes = [1 1 1];
end
cropAxes = cropAxes(:)';
if(numel(cropAxes)~=3), error('numel(cropAxes) ~= 3!'); end

cropAxes(cropAxes>0) = cropAxes(cropAxes>0)+1;
stEl = ones(cropAxes+1);

if(isstruct(roi))
    for(ii=1:numel(roi))
        [roiImg, imgXform, bb] = dtiRoiToImg(roi(ii));
        roiImg = imerode(roiImg, stEl);
        roi(ii).coords = dtiRoiCoordsFromImg(roiImg, imgXform, bb);
    end
else
    [roiImg, imgXform, bb] = dtiRoiToImg(roi);
    roiImg = imerode(roiImg, stEl);
    roi = dtiRoiCoordsFromImg(roiImg, imgXform, bb);
end

return;


text:
function plotFlow(u, v, imgOriginal, rSize, scale)
% Creates a quiver plot that displays the optical flow vectors on the
% original first frame (if provided). See the MATLAB Function Reference for
% "quiver" for more info.
%
% Usage:
% plotFlow(u, v, imgOriginal, rSize, scale)
%
% u and v are the horizontal and vertical optical flow vectors,
% respectively. imgOriginal, if supplied, is the first frame on which the
% flow vectors would be plotted. use an empty matrix '[]' for no image.
% rSize is the size of the region in which one vector is visible. scale
% over-rules the auto scaling.
%
% Author: Mohd Kharbat at Cranfield Defence and Security
% mkharbat(at)ieee(dot)org , http://mohd.kharbat.com
% Published under a Creative Commons Attribution-Non-Commercial-Share Alike
% 3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
% October 2008
% Rev: Jan 2009

figure();

if nargin>2
    if sum(sum(imgOriginal))~=0
        imshow(imgOriginal,[0 255]);
        hold on;
    end
end
if nargin<4
    rSize=5;
end
if nargin<5
    scale=3;
end

% Enhance the quiver plot visually by showing one vector per region
for i=1:size(u,1)
    for j=1:size(u,2)
        if floor(i/rSize)~=i/rSize || floor(j/rSize)~=j/rSize
            u(i,j)=0;
            v(i,j)=0;
        end
    end
end
quiver(u, v, scale, 'color', 'b', 'linewidth', 2);
set(gca,'YDir','reverse');


text:
function [x,y,z] = hilbert3(n)
% Hilbert 3D curve.
%
% function [x,y,z] = hilbert3(n) gives the vector coordinates of points
% in n-th order Hilbert curve of area 1.
%
% Example: plot the 3-rd order curve
%
% [x,y,z] = hilbert3(3); plot3(x,y,z)
%
% Inputs:
%   n  number of subdivision iterations
% Outputs:
%   x  list of x coordinates
%   y  list of y coordinates
%   z  list of z coordinates
%
% Example:
%   n=3;
%   [x,y,z] = hilbert3(n);
%   P = [x;y;z]';
%   [V,F] = wire_mesh(P,[1:size(P,1)-1;2:size(P,1)]','Thickness',2^(1-n)/6);
%   [V,F] = loop(V,F,2);
%   tsurf(F,V);
%

%   Copyright (c) by Ivan Martynov
%   Inspired by function [x,y] = hilbert(n) made by Federico Forte
%   Date: September 17, 2009

if nargin ~= 1
    n = 2;
end

if n <= 0
    x = 0;
    y = 0;
    z = 0;
else
    [xo,yo,zo] = hilbert3(n-1);
    x = .5*[.5+zo .5+yo -.5+yo -.5-xo -.5-xo -.5-yo .5-yo .5+zo];
    y = .5*[.5+xo .5+zo .5+zo .5+yo -.5+yo -.5-zo -.5-zo -.5-xo];
    z = .5*[.5+yo -.5+xo -.5+xo .5-zo .5-zo -.5+xo -.5+xo .5-yo];
end

%{
Copyright (c) 2009, Ivan Martynov
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * 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.
%}


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Returns vector of muscle activation forces for a single Lag.Pt.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Fm = give_3_Element_Muscle_Activation(v,LF,LFO,SK,a,b,Fmax,current_time,xPt,xLag)

% current_time: current time in simulation (s)
% xLag: vector of all x-Lagrangian Pts

% 3 ELEMENT HILL MUSCLE MODEL: contractile element (CE), series element
%                              (SE), parallel element (PE)
%
% NOTE:
%       (1) F_tot = F_PE + F_SE; F_CE = F_SE
%       (2) L = L_PE = L_CE + L_SE
%
% Fm = F_PE + F_CE = F_PE + FSE
%    = k_PE*(Lf - LR) + a_f * Fmax * F1(Lf) * F2(Vf)
%    = k_PE*(Lf - LR) + a_f * Fmax *exp( -( (Q-1)/SK )^2 ) * (1/P0)*(b*P0-a*v)/(v+b); Q = LF/LFO
%

% a_f: coefficient that triggers contraction (traveling square or Gaussian wave?)

% Length-Tension Model Parameters (F1)
% Fmax: maximum isometric force produced at the optimum length of the muscle fibers
% LF:   length of the muscle fibers
% LFO:  length at which the muscle fibers exert their maximum tension
% SK:   constant specific for each muscle where SK > 0.28.

% Hill Model (Force-Velocity) Parameters (F2)
% P0:   maximum load w/ NO contraction
% a:
% b:
% v:    velocity of muscle expansion/contraction

% Length Tension Model Parameters %
Q = LF/LFO;
F1 = exp( - ( (Q-1)/SK )^2 );

% Hill Model %
P0 = Fmax;   %Same as Fmax
F2 = (1/P0)*(b*P0-a*v)/(v+b);

% Get Activation Coefficient %
af_Val = give_Traveling_Triggering_Coefficient(current_time,xLag,xPt);

% Actually Compute Muscle Force %
for i=1:length(xLag)
   Fm = af_Val*Fmax*F1*F2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Returns value of Activation Trigger at specific time
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function af_Val = give_Traveling_Triggering_Coefficient(current_time,xLag,xPt)

% current_time: current_time in simulation
% xLag:         x-Lagrangian pts (both bottom and top of tube)
% xPt: x-Pt of interest

t = current_time;               % current time
freq = 10;                      % frequency of traveling wave down tube
Ltube = xLag(end/2) - xLag(1);  % length of heart tube
buff = 0.25;                    % percent-buff on each end of heart tube
L_AR = Ltube*(1-2*buff);        % length of the actual activation region of tube
SqWidth = L_AR/10;              % width of the square traveling activation wave
v = (L_AR-SqWidth) * freq;      % traveling wave velocity
k = 2*pi*freq/v;                % Wave-number for traveling wave

t = rem(t,1/freq);              % Gives remainder after "modular arithmetic" ("fmod" in C++)

xL = xLag(1) + (buff)*Ltube;
xR = xLag(1) + (buff)*Ltube + L_AR;
    x = xPt;
    if ( ( x >= xL ) && ( x <= xR ) )
        af_Val = 0.75*( sin( 1/(L_AR)*( 2*pi*freq*t + k*x) - pi/2 ) +1 );
    else
        af_Val = 0.0;
    end

xL = xLag(1) + (buff)*Ltube + (v*t);
xR = xLag(1) + (buff)*Ltube + SqWidth + (v*t);
xM = (xL+xR)/2;
c = SqWidth/1.5;
    x = xPt;
    if ( ( x >= xL ) && ( x <= xR ) )
        af_Val = 1;                         % Traveling Square Wave
        %af_Val = exp( -(x-xM)^2 / (2*c)^2 ); % Traveling Gaussian Wave
    else
        af_Val = 0.0;
    end

author: nickabattista


text:
function  status = writeascii(filename, l)

% status = writeascii(filename, l)
% write each line of list l into an ASCII file filename

status = 0;
fid = fopen(filename, 'w');

if fid ~= -1
	for i=1:length(l)
		line = get(l,i);
		fprintf(fid, ['%s \n'], line);
	end
	fclose(fid);
else
	status = -1;
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [r, R_r, R_v, R_a] = rpredict(r, v, a, dt)

% RPREDICT Position prediction.
%   RPREDICT(R,V,DT) performs the time update R = R + V*DT.
%
%   RPREDICT(R,V,A,DT) considers R = R + V*DT + 1/2*A*DT instead.
%
%   [R,R_r,R_v,R_a] = ... returns Jacobian matrices wrt position R,
%   velocity V and acceleration A.
%
%   See also VPREDICT, QPREDICT.

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

if nargin == 3
    dt  = a;
    r   = r + v*dt;
    R_r = eye(length(r));
    R_v = dt*R_r;
    R_a = zeros(3);

elseif nargin == 4
    r   = r + v*dt + .5*a*dt^2;
    R_r = eye(length(r));
    R_v = R_r*dt;
    R_a = 0.5*R_r*dt^2;

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:
function value = c8_le_l2 ( x, y )

%*****************************************************************************80
%
%% C8_LE_L2 := X <= Y for complex values, and the L2 norm.
%
%  Definition:
%
%    The L2 norm can be defined here as:
%
%      C8_NORM_L2(X) = sqrt ( ( real (X) )**2 + ( imag (X) )**2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex X, Y, the values to be compared.
%
%    Output, logical VALUE, is TRUE if X <= Y.
%
  if ( ( real ( x ) )^2 + ( imag ( x ) )^2 ...
    <= ( real ( y ) )^2 + ( imag ( y ) )^2 )
    value = 1;
  else
    value = 0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [d_gospa, x_to_y_assignment, decomposed_cost] = ...
    GOSPA(x_mat, y_mat, p, c, alpha)
% AUTHOR: Abu Sajana Rahmathullah
% DATE OF CREATION: 7 August, 2017
%
%  [d_gospa, x_to_y_assignment] = GOSPA(x_mat, y_mat, p, c, alpha)
% computes the generalized optimal sub-pattern assignment metric (GOSPA)
% metric between the two finite sets, x_mat and y_mat for the given
% parameters c, p and alpha. Note that this implementation is based on
% auction algortihm, implementation of which is also available in this
% repository. For details about the metric, check
% https://arxiv.org/abs/1601.05585. Also visit https://youtu.be/M79GTTytvCM
% for a 15-min presentation about the metric.
%
% INPUT:
%   x_mat, y_mat: Input sets represented as real matrices, where the
%                 columns represent the vectors in the set
%   p           : 1<=p<infty, exponent
%   c           : c>0, cutoff distance
%   alpha       : 0<alpha<=2, factor for the cardinality penalty.
%                 Recommended value 2 => Penalty on missed & false targets
%
% OUTPUT:
%   d_gospa          : Scalar, GOSPA distance between x_mat and y_mat
%   x_to_y_assignment: Integer vector, of length same as the number of
%                      columns of x_mat. The i^th entry in the vector
%                      denotes the column index of the vector in y_mat,
%                      that is assigned to the vector in the i^th column of
%                      x_mat. Note that these indices are based on the
%                      permutation. Therefore, if #columns in x_mat <=
%                      #columns in y_mat, the entries in this vector will
%                      be between 1 and #columns in y_mat. Otherwise, the
%                      entries will be between 0 and #columns in y_mat,
%                      where 0 indicates that the corresponding columns in
%                      x_mat are unassigned.
%   decomposed_cost : Struct that returns the decomposition of the GOSPA
%                     metric for alpha=2 into 3 components:
%                          'localisation', 'missed', 'false'.
%                     Note that
%                     d_gospa = (decomposed_cost.localisation +
%                                decomposed_cost.missed       +
%                                decomposed_cost.false)^(1/p)
%
% Note: Euclidean base distance between the vectors in x_mat and y_mat is
% used in this function. One can change the function 'computeBaseDistance'
% in this function for other choices.

% check that the input parameters are within the valid range

n_ouput_arg=nargout;

checkInput();

nx = size(x_mat, 2); % no of points in x_mat
ny = size(y_mat, 2); % no of points in y_mat

% compute cost matrix
cost_mat = zeros(nx, ny);
for ix = 1:nx
    for iy = 1:ny
        cost_mat(ix, iy) ...
            = min(computeBaseDistance(x_mat(:, ix), y_mat(:, iy)), c);
    end
end

% intialise output values
decomposed_cost     = struct( ...
    'localisation', 0, ...
    'missed',       0, ...
    'false',        0);

x_to_y_assignment   = [];
opt_cost            = 0;

dummy_cost = (c^p) / alpha; % penalty for the cardinality mismatch

% below, cost is negated to make it compatible with auction algorithm
if nx == 0 % when x_mat is empty, all entries in y_mat are false
    opt_cost              = -ny * dummy_cost;
    decomposed_cost.false = opt_cost;
else
    if ny == 0 % when y_mat is empty, all entries in x_mat are missed
        opt_cost               = -nx * dummy_cost;

        if(alpha==2)
            decomposed_cost.missed = opt_cost;
        end
    else % when both x_mat and y_mat are non-empty, use auction algorithm
        cost_mat = -(cost_mat.^p);
        [x_to_y_assignment, y_to_x_assignment, ~] ...
            = auctionAlgortihm(cost_mat, 10*(nx * ny));
        % use the assignments to compute the cost
        for ind = 1:nx
            if x_to_y_assignment(ind) ~= 0
                opt_cost = opt_cost + cost_mat(ind,x_to_y_assignment(ind));

                if(alpha==2)

                    decomposed_cost.localisation = ...
                        decomposed_cost.localisation ...
                        + cost_mat(ind,x_to_y_assignment(ind)) ...
                        .* double(cost_mat(ind,x_to_y_assignment(ind)) > -c^p);

                    decomposed_cost.missed      = decomposed_cost.missed ...
                        - dummy_cost ...
                        .* double(cost_mat(ind,x_to_y_assignment(ind)) == -c^p);

                    decomposed_cost.false       = ...
                        decomposed_cost.false ...
                        - dummy_cost ...
                        .* double(cost_mat(ind,x_to_y_assignment(ind)) == -c^p);
                end
            else
                opt_cost               = opt_cost - dummy_cost;
                if(alpha==2)
                    decomposed_cost.missed = decomposed_cost.missed - dummy_cost;
                end
            end
        end
        opt_cost = opt_cost - sum(y_to_x_assignment == 0) * dummy_cost;
        if(alpha==2)
            decomposed_cost.false = decomposed_cost.false ...
                - sum(y_to_x_assignment == 0) * dummy_cost;
        end
    end
end

% final output
d_gospa                      = (-opt_cost)^(1/p);
decomposed_cost.localisation = (-decomposed_cost.localisation);
decomposed_cost.missed       = (-decomposed_cost.missed);
decomposed_cost.false        = (-decomposed_cost.false);

 function checkInput()
        if size(x_mat, 1) ~= size(y_mat, 1)
            error('The number of rows in x_mat & y_mat should be equal.');
        end
        if ~((p >= 1) && (p < inf))
            error('The value of exponent p should be within [1,inf).');
        end
        if ~(c>0)
            error('The value of base distance c should be larger than 0.');
        end

        if ~((alpha > 0) && (alpha <= 2))
            error('The value of alpha should be within (0,2].');
        end
        if alpha ~= 2 && n_ouput_arg==3
            warning(['decomposed_cost is not valid for alpha = ' ...
                num2str(alpha)]);
        end
    end
end

function db = computeBaseDistance(x_vec, y_vec)
db = sum((x_vec - y_vec).^2)^(1/2);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function res = vl_simplenn(net, x, dzdy, res, varargin)
%VL_SIMPLENN  Evaluate a SimpleNN network.
%   RES = VL_SIMPLENN(NET, X) evaluates the convnet NET on data X.
%   RES = VL_SIMPLENN(NET, X, DZDY) evaluates the convnent NET and its
%   derivative on data X and output derivative DZDY (foward+bacwkard pass).
%   RES = VL_SIMPLENN(NET, X, [], RES) evaluates the NET on X reusing the
%   structure RES.
%   RES = VL_SIMPLENN(NET, X, DZDY, RES) evaluates the NET on X and its
%   derivatives reusing the structure RES.
%
%   This function process networks using the SimpleNN wrapper
%   format. Such networks are 'simple' in the sense that they consist
%   of a linear sequence of computational layers. You can use the
%   `dagnn.DagNN` wrapper for more complex topologies, or write your
%   own wrapper around MatConvNet computational blocks for even
%   greater flexibility.
%
%   The format of the network structure NET and of the result
%   structure RES are described in some detail below. Most networks
%   expect the input data X to be standardized, for example by
%   rescaling the input image(s) and subtracting a mean. Doing so is
%   left to the user, but information on how to do this is usually
%   contained in the `net.meta` field of the NET structure (see
%   below).
%
%   The NET structure needs to be updated as new features are
%   introduced in MatConvNet; use the `VL_SIMPLENN_TIDY()` function
%   to make an old network current, as well as to cleanup and check
%   the structure of an existing network.
%
%   Networks can run either on the CPU or GPU. Use VL_SIMPLENN_MOVE()
%   to move the network parameters between these devices.
%
%   To print or obtain summary of the network structure, use the
%   VL_SIMPLENN_DISPLAY() function.
%
%   VL_SIMPLENN(NET, X, DZDY, RES, 'OPT', VAL, ...) takes the following
%   options:
%
%   `Mode`:: `'normal'`
%      Specifies the mode of operation. It can be either `'normal'` or
%      `'test'`. In test mode, dropout and batch-normalization are
%      bypassed. Note that, when a network is deployed, it may be
%      preferable to *remove* such blocks altogether.
%
%   `ConserveMemory`:: `false`
%      Aggressively delete intermediate results. This in practice has
%      a very small performance hit and allows training much larger
%      models. However, it can be useful to disable it for
%      debugging. Keeps the values in `res(1)` (input) and `res(end)`
%      (output) with the outputs of `loss` and `softmaxloss` layers.
%      It is also possible to preserve individual layer outputs
%      by setting `net.layers{...}.precious` to `true`.
%      For back-propagation, keeps only the derivatives with respect to
%      weights.
%
%   `CuDNN`:: `true`
%      Use CuDNN when available.
%
%   `Accumulate`:: `false`
%      Accumulate gradients in back-propagation instead of rewriting
%      them. This is useful to break the computation in sub-batches.
%      The gradients are accumulated to the provided RES structure
%      (i.e. to call VL_SIMPLENN(NET, X, DZDY, RES, ...).
%
%   `BackPropDepth`:: `inf`
%      Limit the back-propagation to top-N layers.
%
%   `SkipForward`:: `false`
%      Reuse the output values from the provided RES structure and compute
%      only the derivatives (backward pass).
%
%   ## The result format
%
%   SimpleNN returns the result of its calculations in the RES
%   structure array. RES(1) contains the input to the network, while
%   RES(2), RES(3), ... contain the output of each layer, from first
%   to last. Each entry has the following fields:
%
%   - `res(i+1).x`: the output of layer `i`. Hence `res(1).x` is the
%     network input.
%
%   - `res(i+1).aux`: any auxiliary output data of layer i. For example,
%     dropout uses this field to store the dropout mask.
%
%   - `res(i+1).dzdx`: the derivative of the network output relative
%     to the output of layer `i`. In particular `res(1).dzdx` is the
%     derivative of the network output with respect to the network
%     input.
%
%   - `res(i+1).dzdw`: a cell array containing the derivatives of the
%     network output relative to the parameters of layer `i`. It can
%     be a cell array for multiple parameters.
%
%   ## The network format
%
%   The network is represented by the NET structure, which contains
%   two fields:
%
%   - `net.layers` is a cell array with the CNN layers.
%
%   - `net.meta` is a grab-bag of auxiliary application-dependent
%     information, including for example details on how to normalize
%     input data, the class names for a classifiers, or details of
%     the learning algorithm. The content of this field is ignored by
%     VL_SIMPLENN().
%
%   SimpleNN is aware of the following layers:
%
%   Convolution layer::
%     The convolution layer wraps VL_NNCONV(). It has fields:
%
%     - `layer.type` contains the string `'conv'`.
%     - `layer.weights` is a cell array with filters and biases.
%     - `layer.stride` is the sampling stride (e.g. 1).
%     - `layer.pad` is the padding (e.g. 0).
%     - `layer.dilate` is the dilation factor (e.g. 1).
%
%   Convolution transpose layer::
%     The convolution transpose layer wraps VL_NNCONVT(). It has fields:
%
%     - `layer.type` contains the string `'convt'`.
%     - `layer.weights` is a cell array with filters and biases.
%     - `layer.upsample` is the upsampling factor (e.g. 1).
%     - `layer.crop` is the amount of output cropping (e.g. 0).
%
%   Max pooling layer::
%     The max pooling layer wraps VL_NNPOOL(). It has fields:
%
%     - `layer.type` contains the string `'pool'`.
%     - `layer.method` is the pooling method (either 'max' or 'avg').
%     - `layer.pool` is the pooling size (e.g. 3).
%     - `layer.stride` is the sampling stride (usually 1).
%     - `layer.pad` is the padding (usually 0).
%
%   Normalization (LRN) layer::
%     The normalization layer wraps VL_NNNORMALIZE(). It has fields:
%
%     - `layer.type` contains the string `'normalize'` or `'lrn'`.
%     - `layer.param` contains the normalization parameters (see VL_NNNORMALIZE()).
%
%   Spatial normalization layer::
%     The spatial normalization layer wraps VL_NNSPNORM(). It has fields:
%
%     - `layer.type` contains the string `'spnorm'`.
%     - `layer.param` contains the normalization parameters (see VL_NNSPNORM()).
%
%   Batch normalization layer::
%     This layer wraps VL_NNBNORM(). It has fields:
%
%     - `layer.type` contains the string `'bnorm'`.
%     - `layer.weights` contains is a cell-array with, multiplier and
%       biases, and moments parameters
%
%     Note that moments are used only in `'test'` mode to bypass batch
%     normalization.
%
%   ReLU and Sigmoid layers::
%     The ReLU layer wraps VL_NNRELU(). It has fields:
%
%     - `layer.type` contains the string `'relu'`.
%     - `layer.leak` is the leak factor (e.g. 0).
%
%     The sigmoid layer is the same, but for the sigmoid function,
%     with `relu` replaced by `sigmoid` and no leak factor.
%
%   Dropout layer::
%     The dropout layer wraps VL_NNDROPOUT(). It has fields:
%
%     - `layer.type` contains the string `'dropout'`.
%     - `layer.rate` is the dropout rate (e.g. 0.5).
%
%     Note that the block is bypassed in `test` mode.
%
%   Softmax layer::
%     The softmax layer wraps VL_NNSOFTMAX(). It has fields
%
%     - `layer.type` contains the string`'softmax'`.
%
%   Log-loss layer and softmax-log-loss::
%     The log-loss layer wraps VL_NNLOSS(). It has fields:
%
%     - `layer.type` contains `'loss'`.
%     - `layer.class` contains the ground-truth class labels.
%
%     The softmax-log-loss layer wraps VL_NNSOFTMAXLOSS() instead. it
%     has the same parameters, but `type` contains the `'softmaxloss'`
%     string.
%
%   P-dist layer::
%     The p-dist layer wraps VL_NNPDIST(). It has fields:
%
%     - `layer.type` contains the string  `'pdist'`.
%     - `layer.p` is the P parameter of the P-distance (e.g. 2).
%     - `layer.noRoot` it tells whether to raise the distance to
%     the P-th power (e.g. `false`).
%     - `layer.epsilon` is the regularization parameter for the derivatives.
%
%   Custom layer::
%     This can be used to specify custom layers.
%
%     - `layer.type` contains the string `'custom'`.
%     - `layer.forward` is  a function handle computing the block.
%     - `layer.backward` is a function handle computing the block derivative.
%
%     The first function is called as
%
%          res(i+1) = layer.forward(layer, res(i), res(i+1))
%
%     where RES is the structure array specified before. The second function is
%     called as
%
%          res(i) = layer.backward(layer, res(i), res(i+1))
%
%     Note that the `layer` structure can contain additional custom
%     fields if needed.
%
%   See also: dagnn.DagNN, VL_SIMPLENN_TIDY(),
%   VL_SIMPLENN_DISPLAY(), VL_SIMPLENN_MOVE().

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.conserveMemory = false ;
opts.sync = false ;
opts.mode = 'normal' ;
opts.accumulate = false ;
opts.cudnn = true ;
opts.backPropDepth = +inf ;
opts.skipForward = false ;
opts.parameterServer = [] ;
opts.holdOn = false ;
opts = vl_argparse(opts, varargin);

n = numel(net.layers) ;
assert(opts.backPropDepth > 0, 'Invalid `backPropDepth` value (!>0)');
backPropLim = max(n - opts.backPropDepth + 1, 1);

if (nargin <= 2) || isempty(dzdy)
  doder = false ;
  if opts.skipForward
    error('simplenn:skipForwardNoBackwPass', ...
      '`skipForward` valid only when backward pass is computed.');
  end
else
  doder = true ;
end

if opts.cudnn
  cudnn = {'CuDNN'} ;
  bnormCudnn = {'NoCuDNN'} ; % ours seems slighty faster
else
  cudnn = {'NoCuDNN'} ;
  bnormCudnn = {'NoCuDNN'} ;
end

switch lower(opts.mode)
  case 'normal'
    testMode = false ;
  case 'test'
    testMode = true ;
  otherwise
    error('Unknown mode ''%s''.', opts. mode) ;
end

gpuMode = isa(x, 'gpuArray') ;

if nargin <= 3 || isempty(res)
  if opts.skipForward
    error('simplenn:skipForwardEmptyRes', ...
    'RES structure must be provided for `skipForward`.');
  end
  res = struct(...
    'x', cell(1,n+1), ...
    'dzdx', cell(1,n+1), ...
    'dzdw', cell(1,n+1), ...
    'aux', cell(1,n+1), ...
    'stats', cell(1,n+1), ...
    'time', num2cell(zeros(1,n+1)), ...
    'backwardTime', num2cell(zeros(1,n+1))) ;
end

if ~opts.skipForward
  res(1).x = x ;
end

% -------------------------------------------------------------------------
%                                                              Forward pass
% -------------------------------------------------------------------------

for i=1:n
  if opts.skipForward, break; end;
  l = net.layers{i} ;
  res(i).time = tic ;
  switch l.type
    case 'conv'
      res(i+1).x = vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, ...
        'pad', l.pad, ...
        'stride', l.stride, ...
        'dilate', l.dilate, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case 'convt'
      res(i+1).x = vl_nnconvt(res(i).x, l.weights{1}, l.weights{2}, ...
        'crop', l.crop, ...
        'upsample', l.upsample, ...
        'numGroups', l.numGroups, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case 'pool'
      res(i+1).x = vl_nnpool(res(i).x, l.pool, ...
        'pad', l.pad, 'stride', l.stride, ...
        'method', l.method, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case {'normalize', 'lrn'}
      res(i+1).x = vl_nnnormalize(res(i).x, l.param) ;

    case 'softmax'
      res(i+1).x = vl_nnsoftmax(res(i).x) ;

    case 'loss'
      res(i+1).x = vl_nnloss(res(i).x, l.class) ;

    case 'softmaxloss'
      res(i+1).x = vl_nnsoftmaxloss(res(i).x, l.class) ;

    case 'relu'
      if l.leak > 0, leak = {'leak', l.leak} ; else leak = {} ; end
      res(i+1).x = vl_nnrelu(res(i).x,[],leak{:}) ;

    case 'sigmoid'
      res(i+1).x = vl_nnsigmoid(res(i).x) ;

    case 'noffset'
      res(i+1).x = vl_nnnoffset(res(i).x, l.param) ;

    case 'spnorm'
      res(i+1).x = vl_nnspnorm(res(i).x, l.param) ;

    case 'dropout'
      if testMode
        res(i+1).x = res(i).x ;
      else
        [res(i+1).x, res(i+1).aux] = vl_nndropout(res(i).x, 'rate', l.rate) ;
      end

    case 'bnorm'
      if testMode
        res(i+1).x = vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, ...
                                'moments', l.weights{3}, ...
                                'epsilon', l.epsilon, ...
                                bnormCudnn{:}) ;
      else
        res(i+1).x = vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, ...
                                'epsilon', l.epsilon, ...
                                bnormCudnn{:}) ;
      end

    case 'pdist'
      res(i+1).x = vl_nnpdist(res(i).x, l.class, l.p, ...
        'noRoot', l.noRoot, ...
        'epsilon', l.epsilon, ...
        'aggregate', l.aggregate, ...
        'instanceWeights', l.instanceWeights) ;

    case 'custom'
      res(i+1) = l.forward(l, res(i), res(i+1)) ;

    otherwise
      error('Unknown layer type ''%s''.', l.type) ;
  end

  % optionally forget intermediate results
  needsBProp = doder && i >= backPropLim;
  forget = opts.conserveMemory && ~needsBProp ;
  if i > 1
    lp = net.layers{i-1} ;
    % forget RELU input, even for BPROP
    forget = forget && (~needsBProp || (strcmp(l.type, 'relu') && ~lp.precious)) ;
    forget = forget && ~(strcmp(lp.type, 'loss') || strcmp(lp.type, 'softmaxloss')) ;
    forget = forget && ~lp.precious ;
  end
  if forget
    res(i).x = [] ;
  end

  if gpuMode && opts.sync
    wait(gpuDevice) ;
  end
  res(i).time = toc(res(i).time) ;
end

% -------------------------------------------------------------------------
%                                                             Backward pass
% -------------------------------------------------------------------------

if doder
  res(n+1).dzdx = dzdy ;
  for i=n:-1:backPropLim
    l = net.layers{i} ;
    res(i).backwardTime = tic ;
    switch l.type

      case 'conv'
        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
          vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx, ...
          'pad', l.pad, ...
          'stride', l.stride, ...
          'dilate', l.dilate, ...
          l.opts{:}, ...
          cudnn{:}) ;

      case 'convt'
        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
          vl_nnconvt(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx, ...
          'crop', l.crop, ...
          'upsample', l.upsample, ...
          'numGroups', l.numGroups, ...
          l.opts{:}, ...
          cudnn{:}) ;

      case 'pool'
        res(i).dzdx = vl_nnpool(res(i).x, l.pool, res(i+1).dzdx, ...
                                'pad', l.pad, 'stride', l.stride, ...
                                'method', l.method, ...
                                l.opts{:}, ...
                                cudnn{:}) ;

      case {'normalize', 'lrn'}
        res(i).dzdx = vl_nnnormalize(res(i).x, l.param, res(i+1).dzdx) ;

      case 'softmax'
        res(i).dzdx = vl_nnsoftmax(res(i).x, res(i+1).dzdx) ;

      case 'loss'
        res(i).dzdx = vl_nnloss(res(i).x, l.class, res(i+1).dzdx) ;

      case 'softmaxloss'
        res(i).dzdx = vl_nnsoftmaxloss(res(i).x, l.class, res(i+1).dzdx) ;

      case 'relu'
        if l.leak > 0, leak = {'leak', l.leak} ; else leak = {} ; end
        if ~isempty(res(i).x)
          res(i).dzdx = vl_nnrelu(res(i).x, res(i+1).dzdx, leak{:}) ;
        else
          % if res(i).x is empty, it has been optimized away, so we use this
          % hack (which works only for ReLU):
          res(i).dzdx = vl_nnrelu(res(i+1).x, res(i+1).dzdx, leak{:}) ;
        end

      case 'sigmoid'
        res(i).dzdx = vl_nnsigmoid(res(i).x, res(i+1).dzdx) ;

      case 'noffset'
        res(i).dzdx = vl_nnnoffset(res(i).x, l.param, res(i+1).dzdx) ;

      case 'spnorm'
        res(i).dzdx = vl_nnspnorm(res(i).x, l.param, res(i+1).dzdx) ;

      case 'dropout'
        if testMode
          res(i).dzdx = res(i+1).dzdx ;
        else
          res(i).dzdx = vl_nndropout(res(i).x, res(i+1).dzdx, ...
                                     'mask', res(i+1).aux) ;
        end

      case 'bnorm'
        [res(i).dzdx, dzdw{1}, dzdw{2}, dzdw{3}] = ...
          vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx, ...
                     'epsilon', l.epsilon, ...
                     bnormCudnn{:}) ;
        % multiply the moments update by the number of images in the batch
        % this is required to make the update additive for subbatches
        % and will eventually be normalized away
        dzdw{3} = dzdw{3} * size(res(i).x,4) ;

      case 'pdist'
        res(i).dzdx = vl_nnpdist(res(i).x, l.class, ...
          l.p, res(i+1).dzdx, ...
          'noRoot', l.noRoot, ...
          'epsilon', l.epsilon, ...
          'aggregate', l.aggregate, ...
          'instanceWeights', l.instanceWeights) ;

      case 'custom'
        res(i) = l.backward(l, res(i), res(i+1)) ;

    end % layers

    switch l.type
      case {'conv', 'convt', 'bnorm'}
        if ~opts.accumulate
          res(i).dzdw = dzdw ;
        else
          for j=1:numel(dzdw)
            res(i).dzdw{j} = res(i).dzdw{j} + dzdw{j} ;
          end
        end
        dzdw = [] ;
        if ~isempty(opts.parameterServer) && ~opts.holdOn
          for j = 1:numel(res(i).dzdw)
            opts.parameterServer.push(sprintf('l%d_%d',i,j),res(i).dzdw{j}) ;
            res(i).dzdw{j} = [] ;
          end
        end
    end
    if opts.conserveMemory && ~net.layers{i}.precious && i ~= n
      res(i+1).dzdx = [] ;
      res(i+1).x = [] ;
    end
    if gpuMode && opts.sync
      wait(gpuDevice) ;
    end
    res(i).backwardTime = toc(res(i).backwardTime) ;
  end
  if i > 1 && i == backPropLim && opts.conserveMemory && ~net.layers{i}.precious
    res(i).dzdx = [] ;
    res(i).x = [] ;
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet