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

%% Syntax
% |[Ep,Vp,Cp]=pillowHex(E,V,C,shrinkFactor);|

%% Description
%

%%
clear; close all; clc;

%%
% Plot settings
fontSize=25;
faceAlpha1=0.25;
edgeWidth=2;
markerSize=35;
cMap=gjet(6);

%% Examples
%

%% Example: Pillowing a hexahedral element

%%
% Creating an example hexahedral element
V=[0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1;]; %nodes
E=1:8; %Element

%%

shrinkFactor=0.5;
[Ep,Vp]=pillowHex(E,V,[],shrinkFactor);

%%
% Visualize results

[F]=element2patch(E);  %Patch data for plotting
[Fp]=element2patch(Ep);  %Patch data for plotting

cFigure;

subplot(1,2,1); hold on;
title('Original hex element','FontSize',fontSize);
gpatch(F,V,'gw','g',0.5,edgeWidth);
% patchNormPlot(F,V,0.25);
plotV(V,'k.','MarkerSize',markerSize);
colormap(cMap);
axisGeom;
axis off;
camlight headlight;

subplot(1,2,2); hold on;
title('Pillowed hex element','FontSize',fontSize);
gpatch(Fp,Vp,'bw','b',0.5,edgeWidth);
% patchNormPlot(Fp,Vp,0.25);
plotV(Vp,'k.','MarkerSize',markerSize);
colormap(cMap);
axisGeom;
axis off;
camlight headlight;

drawnow;

%% Example: Pillowing a set of hexahedral elements

%%
% Creating an example hexahedral element set
boxDim=[5 6 7];
boxEl=[4 5 6];

[meshStruct]=hexMeshBox(boxDim,boxEl);

E=meshStruct.E;
V=meshStruct.V;
C=(1:1:size(E,1))';

%%
%
shrinkFactor=0.5;
[Ep,Vp,Cp]=pillowHex(E,V,C,shrinkFactor);

%%
% Visualize results
[F,CF]=element2patch(E,C);  %Patch data for plotting
[Fp,CFp]=element2patch(Ep,Cp);  %Patch data for plotting

cFigure;

subplot(1,2,1); hold on;
title('Original hex element','FontSize',fontSize);
gpatch(F,V,CF,'k',0.5);
% patchNormPlot(F,V,0.25);
plotV(V,'k.','MarkerSize',markerSize);
colormap(cMap);
axisGeom;
axis off;
camlight headlight;

subplot(1,2,2); hold on;
title('Pillowed hex element','FontSize',fontSize);
gpatch(Fp,Vp,CFp,'k',0.5);
% patchNormPlot(Fp,Vp,0.25);
plotV(Vp,'k.','MarkerSize',markerSize);
colormap(cMap);
axisGeom;
axis off;
camlight headlight;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function a = reshape(a,varargin)
%RESHAPE      Reshape for interval vectors/matrices
%
%   r = reshape(a,vector)  or  r = reshape(a,n1,n2,...)
%
% functionality as Matlab function reshape
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  improved speed
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 06/16/13     S.M. Rump  performance improvement
%

  if a.complex
    a.mid = reshape(a.mid,varargin{:});
    if ~isequal(a.rad,0)
      a.rad = reshape(a.rad,varargin{:});
    end
  else
    a.inf = reshape(a.inf,varargin{:});
    a.sup = reshape(a.sup,varargin{:});
  end

repo: OpenMAS


text:
function [y, dsdx, dsdp] = VBA_sigmoid(x, varargin)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [y, dsdx, dsdp] = VBA_sigmoid(x, [name, value, ...])
% Apply a sigmoid transformation to x
%
% By default, the canonical sigmoid is used:
%
%   y = 1 / (1 + exp(- x))
%
% However, this function can be parametrized using name/value pairs of
% arguments to get, using all options:
%
%   y = offset + scale * 1 / (1 + exp(- slope * (x - center))
%
% IN:
%   - x: values to be transformed
%   - optional key/value pairs or structure (or both) that parametrize the
%     sigmoid transformation:
%       * 'slope'     : inverse-temperature parameter (default = 1)
%       * 'center'    : absissa of the inflexion point (default = 0)
%       * 'scale'     : multiplicative gain of the transformation (default = 1)
%       * 'offset'    : additive gain (default = 0)
%       * 'lapseRate' : shotcut for offset = lapseRate and scale = 1 - 2 * lapseRate
%                       this option is incompatible with offset or scale
%   - other optional key/value pairs or structure (or both)
%       * 'inverse' : reversed transformation, ie. x = VBA_sigmoid(y, opt)
%                       (default = false)
%       * 'finite'  : boundaries that enforce precision to be finite and
%                     derivative to stay numerically non zero (default = 1e-9)
%                     Set to 0 to deactivate
%       * 'derivatives' : cell array of parameter names wrt which dsdp must be
%                         computed (in given order).
%
%
%   Note that if the 'inverse' flag is set to true, the derivatives will
%   not be computed and dsdx = dsdp = [] will be returned instead.
%
% OUT:
%   - y: transformed values, with the same dimension as x.
%   - dsdx: derivative of the transformation wrt. x, taken at each point x.
%     Has the same dimension as x.
%   - dsdp: derivative of the transformation wrt. the parameters specified
%     as arguments. Derivatives are for parameters taken in alphabetical
%     order and aggregated along the first dimension. For example, calling
%     VBA_sigmoid([x1; x2], 'slope', 2, 'center', 3) will return dsdp as:
%     [ d_s/d_center(x1) d_s/d_center(x2) ;
%       d_s/d_slope(x1)  d_s/d_slope(x2)  ]
%     If x is multidimensional, size(dsdp) = [nb_params, size(x,1), size(x,2), ...]
%
% /////////////////////////////////////////////////////////////////////////

%% Globals
% =========================================================================
% truncature for finite sigmoid
epsilon = 1e-9;

%% Shrtcut
% =========================================================================
% quick version!
if nargin == 1 && nargout == 1
    y = epsilon + (1 - 2 * epsilon) ./ (1 + exp (- x));
    return
end

%% Parse arguments
% =========================================================================

% define inputParser
% -------------------------------------------------------------------------
persistent parser;

if isempty (parser)
    parser = getParser ();
end

% parse arguments
% -------------------------------------------------------------------------

parser.parse (varargin{:});
params = parser.Results ;

% apply shortcuts if needed
% -------------------------------------------------------------------------

if ~ ismember ('lapseRate', parser.UsingDefaults)
    % lapseRate is mutually ex
    if any (~ ismember ({'offset', 'scale'}, parser.UsingDefaults))
        error('*** VBA_sigmoid: you can not specify lapseRate and offset or scale at the same time');
    end
    params.offset = params.lapseRate;
    params.scale = 1 - 2 * params.lapseRate;
end

%% Compute transformation
% =========================================================================

%% Inversed case
if params.inverse
    % check that values are valid
    if ~ VBA_isInRange(x, [0 1])
        error('*** VBA_sigmoid: inverse sigmoid inputs must be between 0 and 1');
    end

    % inverse sigmoid transformation
    lx = params.scale * (x - params.offset) .^-1  - 1;
    y = params.center - params.slope ^-1 * log (lx) ;

    % skip derivatives, they are generally not used in inverse case
    dsdx = [];
    dsdp = [];

%% Normal case
else
    % evaluate sigmoid
    % ---------------------------------------------------------------------
    sx = 1 ./ (1 + exp(- params.slope * (x - params.center)));
    y = params.offset + params.scale * sx;

    % ensure finite precision ('finite' flag)
    % ---------------------------------------------------------------------
    if params.finite > 0
        minY = params.offset + epsilon;
        y = max(y,minY);
        maxY = params.offset + params.scale - epsilon;
        y = min(y,maxY);
    end

    % compute derivatives with respect to value
    % ---------------------------------------------------------------------

    % skip if not not required
    if nargout < 2
        return
    end

    % actual computation
    dsdx = params.slope * (y - params.offset) .* (1 - (y - params.offset) ./ params.scale);

    % compute derivatives with respect to parameters
    % ---------------------------------------------------------------------

    % skip if not not required
    if nargout < 3
        return
    end

    % concatenate derivatives for all parameters
    dims = size(x);

    dsdp = cat (2, ...
        - VBA_vec (dsdx), ... % d_center
        1 - 2 * VBA_vec (sx), ... d_lapseRate
        ones(numel(sx),1), ... d_offset
        VBA_vec (sx), ... d_scale
        ((VBA_vec (x) - VBA_vec(params.center)) / params.slope) .* VBA_vec (dsdx) ... d_slope
        );

    % keep only those passed s parameter
    derivables = {'center','lapseRate','offset','scale','slope'};
    if isempty(params.derivatives)
         params.derivatives = setdiff(derivables, parser.UsingDefaults);
    end
    dIdx = cellfun(@(l) find(strcmp (derivables, l)), params.derivatives);

    dsdp = dsdp(:,dIdx);

    % set derived parameter as first dimension
    if all(dims == 1)
        dsdp = dsdp';
    else
        dims(dims==1) = [];
        dsdp = reshape(dsdp', [size(dsdp,2) dims]);
    end

end

end

function parser = getParser()

parser = inputParser;
parser.PartialMatching = false;
parser.KeepUnmatched = true;

% define parameters
% -------------------------------------------------------------------------

% flags
parser.addParameter ('inverse', false, @islogical);
parser.addParameter ('finite', true, @(z) VBA_isInRange(z, [0 1e-2]));
parser.addParameter ('derivatives', {}, @iscellstr);

% x transfomations
parser.addParameter ('slope', 1, @isnumeric);
parser.addParameter ('center', 0, @isnumeric);

% sig transformation
parser.addParameter ('offset', 0, @isnumeric);
parser.addParameter ('scale', 1, @(z) VBA_isInRange(z, [eps Inf]));

% shortcut for lapse rate model
parser.addParameter ('lapseRate', 0, @(z) VBA_isInRange(z, [0 0.5]));
end

author: MBB-team

repo: VBA-toolbox


text:
% Demonstration of the degree elevation algorithm.
%

crv = nrbtestcrv;

% plot the control points
plot(crv.coefs(1,:),crv.coefs(2,:),'bo')
title('Degree elevation of test curve by 1');
hold on;
plot(crv.coefs(1,:),crv.coefs(2,:),'b--');

% draw nurbs curve
nrbplot(crv,48);

% degree elevate the curve by 1
icrv = nrbdegelev(crv, 1);

% insert new knots and plot new control points
plot(icrv.coefs(1,:),icrv.coefs(2,:),'ro')
plot(icrv.coefs(1,:),icrv.coefs(2,:),'r--');

hold off;

author: Sable

repo: mcbench-benchmarks


text:
function [str] = strengths_und(CIJ)
%STRENGTHS_UND        Strength
%
%   str = strengths_und(CIJ);
%
%   Node strength is the sum of weights of links connected to the node.
%
%   Input:      CIJ,    undirected weighted connection matrix
%
%   Output:     str,    node strength
%
%
%   Olaf Sporns, Indiana University, 2002/2006/2008

% compute strengths
str = sum(CIJ);        % strength

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% compute X*Y on GPU in batches to avoid memory errors
function out = gpuBlockXY(X, Y)

% this code is for a very tall X, and small Y

g = gpuDevice;
bytesGPU = g.AvailableMemory;
infoX = whos('X');
infoY = whos('Y');
[hX, wX] = size(X);
[hY, wY] = size(Y); % hY should be equal to wX
% the dimensions of out will be hX-by-wY
bytesPerRowX = infoX.bytes/hX;
bytesPerRowY = infoY.bytes/hY;
% theoretically, to keep the tmp result, X, and Y in GPU memory
% bytesRequired = bytesPerRowY*batchSize + bytesPerRowX*batchSize + infoY.bytes;
batchSize = floor((bytesGPU-infoY.bytes)/(bytesPerRowX+bytesPerRowY));
batchSize = floor(batchSize/8); % just to be on the safe side
% also, interestingly with smaller batches the code runs slightly faster
nBatches = ceil(hX/batchSize);

startIdx = 1:batchSize:hX;
endIdx = startIdx + batchSize-1;
endIdx = min(endIdx, hX);

out = zeros([hX, wY], 'like', X);
gpuY = gpuArray(Y);
for iBatch = 1:nBatches
%     fprintf('Batch %d/%d\n', iBatch, nBatches);
    batchX = gpuArray(X(startIdx(iBatch):endIdx(iBatch), :));
    tmp = batchX*gpuY;

% MK - this is an annoying line of code, but it helps Matlab to do memory
% management on the GPU properly. Otherwise I get OUT OF MEMORY errors
% maybe it is just my GPU is a bit unstable
    wait(g);
    out(startIdx(iBatch):endIdx(iBatch), :) = gather(tmp);
end

return

author: cortex-lab


text:
function [] = showTSDF(sdf, grid)
sdf(or(sdf <= grid.disp_range(1), sdf >= grid.disp_range(2))) = nan;

h = slice(grid.x, grid.y, grid.z, sdf, [], [], grid.z);

set(h, 'EdgeColor','none', 'FaceColor','interp')
set(h, 'EdgeColor','none')
alpha(.2)
axis equal
end

author: ChirikjianLab

repo: Marching-Primitives


text:
function C = PruneTree(C)

% Logic:
% Start with a clique, scan through its neighbors. If you find a neighbor
% such that it is a superset of the clique you started with, then you know
% that you can prune the tree. For instance, let's take the following
% clique tree:
% ABE -- AB ---AD
% Let's say we started with AB. We scan through its neighbors and find that
% AB is a subset of ABE. So we cut off the edges connected to AB and add an
% edge between ABE and all of AB's other neighbors. This maintains the
% running intersection property and gives us a more compact clique tree
% which looks like: ABE -- AD.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

toRemove = [];

for i=1:length(C.nodes)

    if ismember(i,toRemove), continue, end;
    neighborsI = find(C.edges(i,:));

    for c = 1: length(neighborsI),

        j = neighborsI(c);
        assert(i ~= j);

        if ismember(j,toRemove), continue, end;

        if (sum(ismember(C.nodes{i}, C.nodes{j})) == length(C.nodes{i}))

            for nk = neighborsI

                % find neighbors and connect with that.
                if length(intersect(C.nodes{i}, C.nodes{nk})) == length(C.nodes{i})
                    C.edges(setdiff(neighborsI,[nk]),nk) = 1;
                    C.edges(nk,setdiff(neighborsI,[nk])) = 1;
                    break;
                end
            end

            % kill the edges for the clique that is to be removed.
            C.edges(i,:) = 0;
            C.edges(:,i) = 0;
            toRemove = [i toRemove];

        end
    end
end

toKeep = setdiff(1:length(C.nodes),toRemove);

C.nodes(toRemove) = [];

if isfield(C, 'edges')
    C.edges = C.edges(toKeep,toKeep);
else
    C.edges = [];
end

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [Kd] = arap_linear_block(varargin)
  % ARAP_LINEAR_BLOCK constructs a block of the matrix which constructs the
  % linear terms of a given arap energy. When treating rotations as knowns
  % (arranged in a column) then this constructs Kd of K such that the linear
  % portion of the energy is as a column:
  %   K * R = [Kx Z  ... Ky Z  ...
  %            Z  Kx ... Z  Ky ...
  %            ... ]
  % These blocks are also used to build the "covariance scatter matrices". Here
  % we want to build a scatter matrix that multiplies against positions
  % (treated as known) producing covariance matrices to fit each rotation.
  % Notice that in the case of the RHS of the poisson solve the rotations are
  % known and the positions unknown, and vice versa for rotation fitting. These
  % linear block just relate the rotations to the positions, linearly in each.
  %
  % Kd = arap_linear_block(V,F,d)
  % Kd = arap_linear_block(V,F,d,'ParameterName','ParameterValue,...)
  %
  % Inputs:
  %   V  #V by dim list of initial domain positions
  %   F  #F by #simplex size list of triangle indices into V
  %   d  coordinate of linear constructor to build
  %   Optional:
  %     'Energy'
  %       followed by a string specifying which arap energy definition to use.
  %       One of the following:
  %         'spokes'  "As-rigid-as-possible Surface Modeling" by [Sorkine and
  %           Alexa 2007], rotations defined at vertices affecting incident
  %           edges, default
  %         'elements'  "A local-global approach to mesh parameterization" by
  %           [Liu et al.  2010] or "A simple geometric model for elastic
  %           deformation" by [Chao et al.  2010], rotations defined at
  %           elements (triangles or tets)
  %         'spokes-and-rims'  Adapted version of "As-rigid-as-possible Surface
  %           Modeling" by [Sorkine and Alexa 2007] presented in section 4.2 of
  %           or "A simple geometric model for elastic deformation" by [Chao et
  %           al.  2010], rotations defined at vertices affecting incident
  %           edges and opposite edges
  % Outputs:
  %   Kd  #V by #V/#F block of the linear constructor matrix corresponding to
  %     coordinate d
  %
  % See also: arap, arap_rhs, covariance_scatter_matrix
  %

  % default is Sorkine and Alexa style local rigidity energy
  energy = 'spokes';
  V = varargin{1};
  F = varargin{2};
  d = varargin{3};
  % number of vertices
  n = size(V,1);
  % number of elements
  m = size(F,1);
  % simplex size
  simplex_size = size(F,2);
  assert(simplex_size == 3 || simplex_size == 4);
  % number of dimensions
  dim = size(V,2);
  assert(d <= dim);

  ii = 4;
  while(ii <= nargin)
    switch varargin{ii}
    case 'Energy'
      ii = ii + 1;
      assert(ii<=nargin);
      energy = varargin{ii};
    otherwise
      error(['Unsupported parameter: ' varargin{ii}]);
    end
    ii = ii + 1;
  end

  switch energy
  case 'spokes'
    Kd = spokes_linear_block(V,F,d);
  case 'spokes-and-rims'
    Kd = spokes_and_rims_linear_block(V,F,d);
  case 'elements'
    Kd = elements_linear_block(V,F,d);
  otherwise
    error(['Unsupported energy type: ' energy]);
  end

  function K = spokes_linear_block(V,F,d)
    % Computes a matrix K such that V'* K * R computes
    %  ∑ wij * 0.5 * (V(i,d)-V(j,d)) * (Ri + Rj)
    % j∈N(i)
    %
    % Inputs:
    %   V  #V by dim list of coordinates
    %   F  #F by 3 list of triangle indices into V
    %   d  index into columns of V
    % Output:
    %   K  #V by #V matrix
    %
    E = edges(F);
    % Build upper part of adjacency matrix where instead of a 1 for edge from i
    % to j we have the difference of position in dimension d
    A = sparse(E(:,1),E(:,2),V(E(:,1),d)-V(E(:,2),d),size(V,1),size(V,1));
    % anti-symmetric, or considers direction of edges
    A = A-A';
    % Multiply with cotangent weights (don't worry about diagonal begin wrong
    % since in A it's all zeros
    L = cotmatrix(V,F);
    K = L.*A;
    % correct the diagonal (notice that the sign is positive
    K = K + diag(sum(K,2));
    K = 0.5*K;
  end

  function K = spokes_and_rims_linear_block(V,F,d)
    % Computes a matrix K such that V' * K * R computes
    %  ∑    -2*(cot(aij) + cot(bij) * (V(i,d)-V(j,d)) * (Ri + Rj) +
    %       -2*cot(aij) * (V(i,d)-V(j,d)) * Raij +
    %       -2*cot(bij) * (V(i,d)-V(j,d)) * Rbij
    % j∈N(i)
    %
    % where:          vj
    %              /  |  \
    %             /   |   \
    %            /    |    \
    %           /     |     \
    %          aij    |    bij
    %           \     |     /
    %            \    |    /
    %             \fij|gij/
    %              \  |  /
    %                 vi
    %
    % Inputs:
    %   V  #V by dim list of coordinates
    %   F  #F by 3 list of triangle indices into V
    %   d  index into columns of V
    % Output:
    %   K  #V by #F matrix
    %
    if simplex_size == 3
      % triangles
      C = cotangent(V,F);
      i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
      I = [i1;i2;i2;i3;i3;i1;i1;i2;i3];
      J = [i2;i1;i3;i2;i1;i3;i1;i2;i3];
      v = [ ...
         C(:,3).*(V(i1,d)-V(i2,d)) + C(:,2).*(V(i1,d)-V(i3,d)); ...
        -C(:,3).*(V(i1,d)-V(i2,d)) + C(:,1).*(V(i2,d)-V(i3,d)); ...
         C(:,1).*(V(i2,d)-V(i3,d)) + C(:,3).*(V(i2,d)-V(i1,d)); ...
        -C(:,1).*(V(i2,d)-V(i3,d)) + C(:,2).*(V(i3,d)-V(i1,d)); ...
         C(:,2).*(V(i3,d)-V(i1,d)) + C(:,1).*(V(i3,d)-V(i2,d)); ...
        -C(:,2).*(V(i3,d)-V(i1,d)) + C(:,3).*(V(i1,d)-V(i2,d)); ...
         ... % diagonal
        C(:,3).*(V(i1,d)-V(i2,d)) - C(:,2).*(V(i3,d)-V(i1,d)); ...
        C(:,1).*(V(i2,d)-V(i3,d)) - C(:,3).*(V(i1,d)-V(i2,d)); ...
        C(:,2).*(V(i3,d)-V(i1,d)) - C(:,1).*(V(i2,d)-V(i3,d)); ...
        ];
      % construct and divide by 3 so laplacian can be used as is
      K = sparse(I,J,v,n,n)/3;
    elseif simplex_size == 4
      % tetrahedra
      assert(false)
    end
  end

  function K = elements_linear_block(V,F,d)
    % Computes a matrix K such that V' * K * R computes
    %  ∑    -2*cot(aij) * (V(i,d)-V(j,d)) * Raij +
    %       -2*cot(bij) * (V(i,d)-V(j,d)) * Rbij
    % j∈N(i)
    %
    % where:        vj
    %              / |\
    %             /  | \
    %            aij | bij
    %             \  | /
    %              \ |/
    %               vi
    %
    % Inputs:
    %   V  #V by dim list of coordinates
    %   F  #F by 3 list of triangle indices into V
    %   d  index into columns of V
    % Output:
    %   K  #V by #F matrix
    %
    if simplex_size == 3
      % triangles, #T by 3
      C = cotangent(V,F);
      i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
      I = [i2;i3;i3;i1;i1;i2];
      J = repmat((1:m)',3*2,1);
      v = [ ...
         C(:,1).*(V(i2,d)-V(i3,d)); ...
        -C(:,1).*(V(i2,d)-V(i3,d)); ...
         C(:,2).*(V(i3,d)-V(i1,d)); ...
        -C(:,2).*(V(i3,d)-V(i1,d)); ...
         C(:,3).*(V(i1,d)-V(i2,d)); ...
        -C(:,3).*(V(i1,d)-V(i2,d))];
      K = sparse(I,J,v,n,m);
    elseif simplex_size == 4
      % tetrahedra, #T by 6
      C = cotangent(V,F);
      i1 = F(:,1); i2 = F(:,2); i3 = F(:,3); i4 = F(:,4);
      I = [i2;i3;i3;i1;i1;i2;i4;i1;i4;i2;i4;i3];
      J = repmat((1:m)',6*2,1);
      v = [ ...
         C(:,1).*(V(i2,d)-V(i3,d)); ...
        -C(:,1).*(V(i2,d)-V(i3,d)); ...
         C(:,2).*(V(i3,d)-V(i1,d)); ...
        -C(:,2).*(V(i3,d)-V(i1,d)); ...
         C(:,3).*(V(i1,d)-V(i2,d)); ...
        -C(:,3).*(V(i1,d)-V(i2,d)); ...
         C(:,4).*(V(i4,d)-V(i1,d)); ...
        -C(:,4).*(V(i4,d)-V(i1,d)); ...
         C(:,5).*(V(i4,d)-V(i2,d)); ...
        -C(:,5).*(V(i4,d)-V(i2,d)); ...
         C(:,6).*(V(i4,d)-V(i3,d)); ...
        -C(:,6).*(V(i4,d)-V(i3,d))];
      K = sparse(I,J,v,n,m);
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
function D = gamma_dist()

    D.gamma_expectation = @gamma_expectation;
    function ratio = gamma_expectation(alpha, beta)

        ratio = alpha ./ beta;

    end

    D.gamma_expectation_log = @gamma_expectation_log;
    function val = gamma_expectation_log(alpha, beta)

        val = psi(alpha) - log(beta);

    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [exp_index, exp_id] = load_em_history_data(emHistoryFile)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    %% Import data from text file.
    % Script for importing data from the following text file:

    %% Initialize variables.
    delimiter = ',';
    startRow = 1;

    %% Format string for each line of text:
    %   column1: double (%d)
    %	column2: double (%f)
    %   column3: double (%f)
    %	column4: double (%f)
    %   column5: double (%f)
    %	column6: double (%f)
    %   column7: double (%f)
    % For more information, see the TEXTSCAN documentation.
    formatSpec = '%d%d%[^\n\r]';

    %% Open the text file.
    fileID = fopen(emHistoryFile,'r');

    %% Read columns of data according to format string.
    % This call is based on the structure of the file used to generate this
    % code. If an error occurs for a different file, try regenerating the code
    % from the Import Tool.
    dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);

    %% Close the text file.
    fclose(fileID);

    %% Allocate imported array to column variable names
    exp_index = dataArray{:, 1};
    exp_id = dataArray{:, 2};

    %% Clear temporary variables
    clearvars emHistoryFile delimiter startRow formatSpec fileID dataArray ans;
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [rfAmp,rfPhase,rfFreq,rfCoil,rfTime]=rfTanhTan(p)
%create a Tanh/Tan adiabatic inversion rf pulse starting from tStart and ending at tEnd
%tStart rf start time
%tEnd rf end time
%dt rf sample time
%rfPhase rf phase
%rfFreq rf off-res freq

tStart=p.tStart;
tEnd=p.tEnd;
dt=p.dt;
MaxB1=p.MaxB1; % Maxium B1
TBP=p.TBP; % Time bandwidth product
rfPhase=p.rfPhase;
rfCoil=p.CoilID;
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

tEnd=tEnd-tStart;
tStart=0; % time scale shift to 0

Zeta=10;
Kappa=atan(20);
A=TBP*pi/(tEnd-tStart);

rfTime=linspace(tStart,tEnd,ceil((tEnd-tStart)/dt)+1);

rfTime1=rfTime(rfTime<(tEnd-tStart)/2);
rfAmp1=MaxB1.*tanh((2*Zeta.*rfTime1)/(tEnd-tStart)); % rf amplitude modulation
rfFreq1=A.*(tan(Kappa*(1-2*rfTime1/(tEnd-tStart)))/tan(Kappa))/(2*pi); % rf frequency modulation
rfTime2=(tEnd-tStart)-rfTime(rfTime>=(tEnd-tStart)/2);
rfAmp2=MaxB1.*tanh((2*Zeta.*rfTime2)/(tEnd-tStart)); % rf amplitude modulation
rfFreq2=-A.*(tan(Kappa*(1-2*rfTime2/(tEnd-tStart)))/tan(Kappa))/(2*pi); % rf frequency modulation

rfAmp=[rfAmp1 rfAmp2];
rfFreq=[rfFreq1 rfFreq2];
rfPhase=(rfPhase)*ones(size(rfTime)); % rf Phase

rfTime=rfTime+p.tStart; % time scale shift back
rfCoil=(rfCoil)*ones(size(rfTime));
rfAmp(1)=0;
rfAmp(end)=0;
rfFreq(1)=0;
rfFreq(end)=0;
rfPhase(1)=0;
rfPhase(end)=0;

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    rfAmp=repmat(rfAmp,[1 Duplicates]);
    rfFreq=repmat(rfFreq,[1 Duplicates]);
    rfPhase=repmat(rfPhase,[1 Duplicates]);
    rfCoil=repmat(rfCoil,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(rfTime) 1]);
    rfTime=repmat(rfTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function tests = SE2Test
    tests = functiontests(localfunctions);
end

% we will assume that the primitives rotx,trotx, etc. all work

function constructor_test(tc)

    verifyClass(tc, SE2(), 'SE2');

    %% null
    tc.verifyEqual(SE2().double, eye(3,3));

    %% translation only
    t = [1 2];
    tc.verifyEqual(SE2(t).double, transl2(t));
    tc.verifyEqual(SE2(t').double, transl2(t));
    tc.verifyEqual(SE2(t(1), t(2)).double, transl2(t));

    %% R
    R = rot2(-pi/2);
    tc.verifyEqual(SE2(R).double, r2t(R));

    %% R,t

    tc.verifyEqual(SE2(R,t).double, transl2(t)*r2t(R));

    %% T
    T = transl2(1, 2)*trot2(0.3);
    tc.verifyEqual(SE2(T).double, T);

    %% x,y,theta
    T = transl2(1, 2)*trot2(0.3);
    tc.verifyEqual(SE2(1, 2, 0.3).double, T);
    tc.verifyEqual(SE2([1, 2, 0.3]).double, T);
    tc.verifyEqual(SE2([1, 2], 0.3).double, T);

    %% T
    T = rt2tr(R,t);
    tc.verifyEqual(SE2(T).double, T);

    %% copy constructor
    TT = SE2(T);
    tc.verifyEqual(SE2(TT).double, T);

    %% vectorised versions

    T1 = transl2(1,2) * trot2(0.3);
    T2 = transl2(1,-2) * trot2(-0.4);

    TT = cat(3, T1, T2, T1, T2);

    tt = SE2(TT);
    tc.verifyEqual(length(tt), size(TT, 3) );
    tc.verifyEqual(tt.T, TT);

end

function concat_test(tc)
    x = SE2();
    xx = [x x x x];

    tc.verifyClass(xx, 'SE2');
    tc.verifySize(xx, [1 4]);
end

function staticconstructors_test(tc)

    %% exponential
    tc.verifyEqual(SE2.exp( skew(0.3) ).R, rot2(0.3), 'AbsTol', 1e-10  );

    %% exponential
    tc.verifyEqual(SE2.exp(zeros(3,3)).T, eye(3,3), 'AbsTol', 1e-10  );
    t = [1 2];
    tc.verifyEqual(SE2.exp(skewa([t 0])).T, transl2(t), 'AbsTol', 1e-10  );
end

function isa_test(tc)

    verifyTrue(tc, SE2.isa(trot2(0)) );
    verifyFalse(tc, SE2.isa(1) )
end

function resulttype_test(tc)

    t = SE2();
    verifyClass(tc, t, 'SE2');

    verifyClass(tc, t*t, 'SE2');

    verifyClass(tc, t/t, 'SE2');

    verifyClass(tc, inv(t), 'SE2');
    end

function inverse_test(tc)

    T1 = transl2(1, 2) * trot2(0.3);
    TT1 = SE2(T1);

    % test inverse
    tc.verifyEqual(double(TT1.inv()), inv(T1), 'AbsTol', 1e-10  );

    tc.verifyEqual(double(TT1*TT1.inv()), eye(3,3), 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT1.inv()*TT1), eye(3,3), 'AbsTol', 1e-10  );

    % vector case
    tc.verifyEqual(double(TT1.inv()*TT1), eye(3,3), 'AbsTol', 1e-10  );
end

function Rt_test(tc)

    TT1 = SE2.rand
    T1 = TT1.double; R1 = t2r(T1); t1 = transl2(T1);

    tc.verifyEqual(TT1.T, T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(TT1.R, R1, 'AbsTol', 1e-10  );
    tc.verifyEqual(TT1.t, t1, 'AbsTol', 1e-10  );

    tc.verifyEqual(TT1.transl, t1', 'AbsTol', 1e-10  );
    TT = [TT1 TT1 TT1];
    tc.verifyEqual(TT.transl, [t1 t1 t1]', 'AbsTol', 1e-10  );
end

function arith_test(tc)

    R1 = rpy2r( randn(1,3) );  t1 = randn(3,1); T1 = rt2tr(R1, t1);
    R2 = rpy2r( randn(1,3) );  t2 = randn(3,1); T2 = rt2tr(R2, t2);

    TT1 = SE2.rand; T1 = TT1.double;
    TT2 = SE2.rand; T2 = TT2.double;

    I = SE2();

    %% SE2 * SE2 product
    % scalar x scalar

    tc.verifyEqual(double(TT1*TT2), T1*T2, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT2*TT1), T2*T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT1*I), T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT2*I), T2, 'AbsTol', 1e-10  );

    % vector x vector
    tc.verifyEqual([TT1 TT1 TT2] * [TT2 TT1 TT1], [TT1*TT2 TT1*TT1 TT2*TT1]);

    % scalar x vector
    tc.verifyEqual(TT1 * [TT2 TT1], [TT1*TT2 TT1*TT1]);

    % vector x scalar
    tc.verifyEqual([TT1 TT2]*TT2, [TT1*TT2 TT2*TT2]);

    %% SE2 * vector product
    vx = [1 0]'; vy = [0 1]';

    % scalar x scalar

    tc.verifyEqual(TT1*vy, h2e( T1*e2h(vy) ), 'AbsTol', 1e-10);

    % vector x vector
    tc.verifyEqual([TT1 TT2] * [vx vy], [h2e(T1*e2h(vx)) h2e(T2*e2h(vy))], 'AbsTol', 1e-10);

    % scalar x vector
    tc.verifyEqual(TT1 * [vx vy], h2e( T1*e2h([vx vy]) ), 'AbsTol', 1e-10);

    % vector x scalar
    tc.verifyEqual([TT1 TT2 TT1] * vy, [h2e(T1*e2h(vy)) h2e(T2*e2h(vy)) h2e(T1*e2h(vy))], 'AbsTol', 1e-10);

end

function function_tests(tc)

    % log
    T = SE2.exp([2 3 0.5]);
    tc.verifyEqual(log(T), [0 -0.5 2; 0.5 0 3; 0 0 0], 'AbsTol', 1e-10  );

end

function conversions_test(tc)

    %%  SE2                     convert to SE2 class

    TT = SE2(1, 2, 0.3);

    verifyClass(tc, TT.SE3, 'SE3');
    tc.verifyEqual(double(TT.SE3), transl(1, 2, 0) * trotz(0.3), 'AbsTol', 1e-10 );

    %% xyt
    tc.verifyEqual(TT.xyt(), [1 2, 0.3]', 'AbsTol', 1e-10);

    %% Twist
    T = SE2.exp([2 3 0.5]);
    t = T.Twist();
    verifyInstanceOf(tc, t, 'Twist');
    tc.verifyEqual(t.v, [2 3]', 'AbsTol', 1e-10);
    tc.verifyEqual(t.w, 0.5, 'AbsTol', 1e-10);

    %% Lie stuff
    th = 0.3;
    RR = SO2(th);
    tc.verifyEqual(RR.log, skew(th), 'AbsTol', 1e-10 );

end

function interp_test(tc)
    TT = SE2.rand
    I = SE2;

    z = interp(I, TT, 0);
    tc.verifyClass(z, 'SE2')

    tc.verifyEqual(double(interp(I, TT, 0)),   double(I), 'AbsTol', 1e-10 );
    tc.verifyEqual(double(interp(I, TT, 1)),   double(TT), 'AbsTol', 1e-4 );
    tc.verifyEqual(double(interp(I, TT, 0.5)), double(trinterp2(TT.T, 0.5)), 'AbsTol', 1e-10  );

end

function miscellany_test(tc)

    TT = SE2(1, 2, 0.3);

    tc.verifyEqual(dim(TT), 3);

    tc.verifyEqual(isSE(TT), true );

    tc.verifyClass(TT.new, 'SE2');

    tc.verifyClass(SE2.convert(TT), 'SE2');
    tc.verifyClass(SE2.convert(TT.T), 'SE2');
    z = SE2.convert(TT);
    tc.verifyEqual(double(z), double(TT));

    z = SE2.convert(TT.T);
    tc.verifyEqual(double(z), TT.T);

end

function display_test(tc)

    T1 = SE2.rand;
    T2 = SE2.rand

    T1.print
    trprint2(T1)   % old style syntax

    T1.plot

    T1.print
    trprint2(T1)   % old style syntax

    T1.plot
    trplot2(T1)   % old style syntax

    T1.animate
    T1.animate(T2)
    tranimate2(T1)   % old style syntax
    tranimate2(T1, T2)   % old style syntax
    tranimate2(T1)   % old style syntax
    tranimate2(T1, T2)   % old style syntax
end

author: petercorke

repo: spatialmath-matlab


text:
function [mmHg] = ftH2O2mmHg(ftH2O)
% Convert pressure from feet of water column at 4 degrees to millimeters of
% mercury at 0 degrees C
% Chad Greene 2012
mmHg = ftH2O*22.4198;


text:
function obj = removeConcurrentAnalysisFiles(obj)
	% DESCRIPTION
	% Biafra Ahanonu
	% started: 2014.07.31 - branch from calciumImagingAnalysis 2020.05.07 [15:47:29]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	obj.foldersToAnalyze
	if obj.guiEnabled==1
		if isempty(obj.foldersToAnalyze)
			scnsize = get(0,'ScreenSize');
			[fileIdxArray, ~] = listdlg('ListString',obj.fileIDNameArray,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','which folders to analyze?');
		else
			fileIdxArray = obj.foldersToAnalyze;
		end
	else
		if isempty(obj.foldersToAnalyze)
			fileIdxArray = 1:length(obj.fileIDNameArray);
		else
			fileIdxArray = obj.foldersToAnalyze;
		end
	end

	% Find all concurrent analysis files and delete them
	nFolders = length(fileIdxArray);
	nFoldersTotal = length(obj.inputFolders);
	for folderNo = 1:nFolders
		i = fileIdxArray(folderNo);
		try
			% disp([num2str(i) ' | ' obj.inputFolders{i}])
			fileList = getFileList(obj.inputFolders{i},obj.concurrentAnalysisFilename);
			if ~isempty(fileList)
				fprintf('%d/%d (%d/%d) will delete analysis file: %s\n',i,nFoldersTotal,folderNo,nFolders,fileList{1})
			end
		catch err
			disp(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			disp(repmat('@',1,7))
		end
	end
	userInput = input('Continue (1 = yes, 0 = no)? ');
	if userInput==1
		for folderNo = 1:nFolders
			i = fileIdxArray(folderNo);
			try
				% disp([num2str(i) ' | ' obj.inputFolders{i}])
				fprintf('%d/%d (%d/%d): %s\n',i,nFoldersTotal,folderNo,nFolders,obj.inputFolders{i})
				fileList = getFileList(obj.inputFolders{i},obj.concurrentAnalysisFilename);
				if ~isempty(fileList)
					fprintf('Deleting analysis file: %s\n',fileList{1})
					delete(fileList{1});
				end
			catch err
				disp(repmat('@',1,7))
				disp(getReport(err,'extended','hyperlinks','on'));
				disp(repmat('@',1,7))
			end
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function node_num = sphere_grid_q9_node_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% SPHERE_GRID_Q9_NODE_NUM counts nodes in a Q9 sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.
%
%    Output, integer NODE_NUM, the number of nodes in the grid.
%
  node_num = 4 * nelemx * nelemy - 2 * nelemx + 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = haar ( n, x )

%*****************************************************************************80
%
%% HAAR performs a Haar transform.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2011
%
%  Author:
%
%    Ken Beauchamp
%
%  Reference:
%
%    Ken Beauchamp,
%    Walsh functions and their applications,
%    Academic Press, 1975,
%    ISBN: 0-12-084050-2,
%    LC: QA404.5.B33.
%
%  Parameters:
%
%    Input, integer N, the number of items in X.
%    N must be a power of 2.
%
%    Input, real X(N), the data to be transformed.
%
%    Output, real X(N), the transformed data.
%
  k = i4_log_2 ( n );

  for i = 1 : k

    l = k + 1 - i;
    l2 = 2^( l - 1 );

    y(1:2*l2) = x(1:2*l2);

    for j = 1 : l2
       l3 = l2 + j;
       jj = 2 * j - 1;
       x(j) = y(jj) + y(jj+1);
       x(l3) = y(jj) - y(jj+1);
    end

  end

  return
end


text:
%########################################################################
%
%	- PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : Alpha RA 1.0
%
%########################################################################
%
%	ground_truth_stats.m:
%
% Description:
%
%   A static MATLAB class to compute ground truth statistics
%

classdef ground_truth_stats

   methods(Static)
       function [pearson, rmse, snr, snr_std, bpm] = get(ppg,waveform, fs)
           y=ppg;
           L=length(y);
           NFFT = 2^nextpow2(L);
           f = fs/2*linspace(0,1,NFFT/2+1);

           [spec_ppg]=spectrogram(y,256,240,f,fs,'yaxis');

           [peak_values, peak_frequencies] = max(abs(spec_ppg));
           bpm_ppg=f(peak_frequencies)*60;

           y=waveform;
           L=length(y);
           NFFT = 2^nextpow2(L);
           f = fs/2*linspace(0,1,NFFT/2+1);

           [spec,freqs,tmp, pxx]=spectrogram(y,256,240,f,fs,'yaxis');

           [peak_values, peak_frequencies] = max(abs(spec));
           bpm=f(peak_frequencies)*60;

           num_elements=size(bpm_ppg,2);
           if num_elements>size(bpm,2)
               num_elements=size(bpm,2);
           end

           [pearson]=corr(bpm_ppg(1,1:num_elements)',bpm(1,1:num_elements)');
           [rmse]=sum(sqrt((bpm_ppg(1,1:num_elements)-bpm(1,1:num_elements)).^2))/num_elements;

           %snr
           hr_f=bpm_ppg/60;
           for i=1:num_elements
               GTMask1 = (freqs >= hr_f(1,i)-0.1)&(freqs <= hr_f(1,i)+0.1);
               GTMask2 = (freqs >= hr_f(1,i)*2-0.2)&(freqs <= hr_f(1,i)*2+0.2);
               %SPower = sum(abs(spec(GTMask1|GTMask2,i)));

               %SPower = sum(abs(spec(GTMask1,i)));
               signal_power = sum(pxx(GTMask1,i));
               FMask2 = (freqs >= 0.5)&(freqs <= 4);
               %total_power = sum(abs(spec(FMask2,i)));
               total_power = sum(pxx(FMask2,i));
               SNR(i) = pow2db(signal_power/(total_power-signal_power));
           end

           [snr] = mean(SNR);
           [snr_std] =std(SNR);
       end
   end
end

author: partofthestars

repo: PPGI-Toolbox


text:
%THIS FUMCTION ILLUSTRATES HOW TO USE THE USER-DEFINED ALGORITHM
function [AH,XH] = user_alg5(Y,r,X)
%
% The example of implementing the user-defined algorithm (this is the Lee-Seung algorithm based on the KL divergence)
%
% INPUTS:
% Y - mixed signals (matrix of size [m by T])
% r - number of estimated signals
% X - true source signals
%
% OUTPUTS
% AH - estimated mixing matrix (matrix of size [m by r])
% XH - estimated source signals (matrix of size [r by T])
%
% #########################################################################
% Initialization
[m,T]=size(Y);
Y(Y <=0) = eps; % this enforces the positive value in the data
AH=rand(m,r);
XH=rand(r,T);

IterNo = 1000; % number of alternating steps

% Iterations
for k = 1:IterNo
    XH = XH.*(AH'*(Y./(AH*XH + eps)));
    AH = AH.*((Y./(AH*XH + eps))*XH')./repmat(sum(XH,2)',m,1);
    AH = AH*diag(1./(sum(AH,1) + eps));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [command] = generate_command(drone_type, autopilot_version, time)
% GENERATE_COMMAND - function that generates commands for a fixed_wing or a
% quadcopter drone. This function is used to test the autopilots of the
% drones, the commanded values here are hard-coded.

if drone_type == "fixed_wing"
        % command set
        command = [0 100 0.26 0]';
elseif drone_type == "quadcopter" || drone_type == "point_mass"
        flag = mod(floor(time/3),2);

        switch autopilot_version
            case 1 % attitude control
                if flag == 0
                    command = [0 0 0 0.5]';
                else
                    command = [0 0 0 -0.5]';
                end
            case 2 % velocity controller
                if flag == 0
                    command = [0 0 6 0]';
                else
                    command = [0 0 -6 0]';
                end
            case 3 % acceleration controller
                command = [0.5 0 0 0]';
        end

end

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


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

%%
clear; close all; clc;

%% Syntax
% |[F,V,C]=grid2patch(X,Y,Z,C,perdiocOpt);|

%% Description
% This function generates patch data for gridded surface coordinates
% defined by X,Y,Z (e.g. as produced by functions like meshgrid or
% ndgrid). Optional inputs include the color data C (a color for each
% vertex), perdiocOpt (a 1x2 logic) which sets periodicity of the surface
% in the first or second direction.
% Similar to MATLAB's |surf2patch| function (which does not include
% periodic connectivity).

%% Examples

%%
% Plot settings
cMap=gjet(250);
fontSize=15;

%% Example: Converting gridded surface data to patch data
[X,Y]=ndgrid(linspace(0,2*pi,25));
Z=sin(X)+sin(Y);

%%
% Creating patch data for the cylinger
[F,V,C]=grid2patch(X,Y,Z,Z);

%%
% Visualizing outputs

cFigure; hold on;
gpatch(F,V,C,'k');
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: periodicity in direction 1

n=12;
t=linspace(-pi,pi,n+1);
t=t(1:end-1);
v1=[cos(t(:)) sin(t(:)) zeros(size(t(:))) ];
v2=3*v1;

%Get coordinate matrices
m=5;
X=linspacen(v1(:,1),v2(:,1),m);
Y=linspacen(v1(:,2),v2(:,2),m);
Z=linspacen(v1(:,3),v2(:,3),m);

%Collect node set
VN=[X(:) Y(:) Z(:)];
numSteps=n-1;

%%
% Creating patch data for the cylinger
[F1,V1]=grid2patch(X,Y,Z,[],[0 0]); %Not periodic
[F2,V2]=grid2patch(X,Y,Z,[],[1 0]); %Periodic in first direction

%%
% Visualizing outputs

cFigure;
subplot(1,2,1); hold on;
title('Not periodic')
gpatch(F1,V1,'bw','k');
gpatch(patchBoundary(F1),V1,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

subplot(1,2,2); hold on;
title('Periodic in 1st direction')
gpatch(F2,V2,'bw','k');
gpatch(patchBoundary(F2),V2,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: periodicity in direction 2

%%
% Creating an example of a cylinder
t=linspace(-pi,pi,25);
t=t(1:end-1);
[T,Z]=meshgrid(t,linspace(0,1,5));

[X,Y]=pol2cart(T,ones(size(T)));

%%
% Creating patch data for the cylinger
[F1,V1,C1]=grid2patch(X,Y,Z,T,[0 0]); %Not periodic
[F2,V2,C2]=grid2patch(X,Y,Z,T,[0 1]); %Periodic in first direction

%%
% Visualizing outputs

cFigure;
subplot(1,2,1); hold on;
title('Not periodic')
gpatch(F1,V1,C1,'k');
gpatch(patchBoundary(F1),V1,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

subplot(1,2,2); hold on;
title('Periodic in 2nd direction')
gpatch(F2,V2,C2,'k');
gpatch(patchBoundary(F2),V2,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: periodicity in both directions

%%
% Creating an example of a cylinder
n=12; m=25;
t=linspace(-pi,pi,n+1);
t=t(1:end-1);
v=[cos(t(:)) zeros(size(t(:))) sin(t(:))];
v(:,1)=v(:,1)+3;

a=linspace(0,2*pi,m+1); a=a(1:end-1);
X=zeros(n,m);
Y=zeros(n,m);
Z=zeros(n,m);

for q=1:1:m
    R=euler2DCM([0 0 a(q)]);
    vr=v*R;
   X(:,q)=vr(:,1);
   Y(:,q)=vr(:,2);
   Z(:,q)=vr(:,3);
end

%%
% Creating patch data for the cylinger
[F1,V1,C1]=grid2patch(X,Y,Z,Z,[0 0]); %Not periodic
[F2,V2,C2]=grid2patch(X,Y,Z,Z,[0 1]); %Periodic in first direction
[F3,V3,C3]=grid2patch(X,Y,Z,Z,[1 0]); %Periodic in second direction
[F4,V4,C4]=grid2patch(X,Y,Z,Z,[1 1]); %Periodic both directions

%%
% Visualizing outputs

cFigure;
subplot(2,2,1); hold on;
title('Not periodic')
gpatch(F1,V1,C1,'k');
gpatch(patchBoundary(F1),V1,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

subplot(2,2,2); hold on;
title('Periodic in 1st direction')
gpatch(F2,V2,C2,'k');
gpatch(patchBoundary(F2),V2,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

subplot(2,2,3); hold on;
title('Periodic in 2nd direction')
gpatch(F3,V3,C3,'k');
gpatch(patchBoundary(F3),V3,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

subplot(2,2,4); hold on;
title('Periodic in 1st and 2nd direction')
gpatch(F4,V4,C4,'k');
gpatch(patchBoundary(F4),V4,'none','r',1,3);
axisGeom(gca,fontSize); camlight headlight;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON


text:
function [ap_cascade, ap_dp, speedup] = compare_runtime(model, prec, suffix)

error('Read the README in this file and then remove this line.');

%                        !!! README FIRST !!!
%
% To perform experiments using the methodology in the paper you MUST
% compile the SINGLE-THREADED version of the convolution routine.
% Otherwise you will be comparing a multi-threaded, BLAS accelerated
% version of the DP algorithm to a single-threaded, non-BLAS accelerated
% version of the cascade algorithm.  To compile the single-threaded
% convolution routine simply run this command in matlab:
%
% >> mex -O fconv.cc -o fconv
%
% When you're done, you'll probably want to switch back to the multi-threaded,
% BLAS accelerated version.  You can do this by running the command:
%
% >> mex -O fconvblas.cc -lmwblas -o fconv

globals;
if nargin < 3
  suffix = num2str(prec);
end
[boxes, times] = ...
  dotest(model, 'test', model.year, [model.year '_CASCADE_' suffix], prec);

ap_cascade = pascal_eval(model.class, boxes{1}, 'test', model.year, ['2007_CASCADE_' suffix]);
ap_dp = pascal_eval(model.class, boxes{2}, 'test', model.year, ['2007_DP_' suffix]);
speedup = mean(times(:,2)./times(:,1));

function [boxes, times] = dotest(model, testset, year, suffix, prec)

setVOCyear = year;
globals;
pascal_init;

pca = 5;
cls = model.class;
ids = textread(sprintf(VOCopts.imgsetpath, testset), '%s');

[prec, recall, thresh] = cascade_thresh(model, year, prec);
model.thresh = thresh;
pca = 5;
cscmodel = cascade_model(model, model.year, pca, model.thresh);

times = zeros(length(ids), 3);
turns = round(rand(length(ids),1));

% run detector in each image
for i = 1:length(ids);
  if strcmp('inriaperson', cls)
    % INRIA uses a mixutre of PNGs and JPGs, so we need to use the annotation
    % to locate the image.  The annotation is not generally available for PASCAL
    % test data (e.g., 2009 test), so this method can fail for PASCAL.
    rec = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    im = imread([VOCopts.datadir rec.imgname]);
  else
    im = imread(sprintf(VOCopts.imgpath, ids{i}));
  end

  % time feature computation
  th = tic();
  pyra = featpyramid(im, model);
  time_feat = toc(th);

  % randomize order to wash out who-went-first cache effects
  if turns(i)
    th = tic();
    [dets_CSC, boxes_CSC] = cascade_detect(pyra, cscmodel, cscmodel.thresh);
    time_CSC = toc(th);

    th = tic();
    [dets_DP, boxes_DP] = gdetect(pyra, model, model.thresh);
    time_DP = toc(th);
  else
    th = tic();
    [dets_DP, boxes_DP] = gdetect(pyra, model, model.thresh);
    time_DP = toc(th);

    th = tic();
    [dets_CSC, boxes_CSC] = cascade_detect(pyra, cscmodel, cscmodel.thresh);
    time_CSC = toc(th);
  end

  if ~isempty(boxes_CSC)
    [dets_CSC boxes_CSC] = clipboxes(im, dets_CSC, boxes_CSC);
    I = nms(dets_CSC, 0.5);
    boxes{1}{i} = dets_CSC(I,[1:4 end]);
    parts{1}{i} = boxes_CSC(I,:);
  else
    boxes{1}{i} = [];
    parts{1}{i} = [];
  end

  if ~isempty(boxes_DP)
    boxes_DP = reduceboxes(model, boxes_DP);
    [dets_DP boxes_DP] = clipboxes(im, dets_DP, boxes_DP);
    I = nms(dets_DP, 0.5);
    boxes{2}{i} = dets_DP(I,[1:4 end]);
    parts{2}{i} = boxes_DP(I,:);
  else
    boxes{2}{i} = [];
    parts{2}{i} = [];
  end

  times(i,:) = [time_CSC time_DP time_feat];
  ravg = mean(times(1:i,2)./times(1:i,1));
  fprintf('%s: testing: %s %s, %d/%d (%.3f)\n', ...
          cls, testset, VOCyear, i, length(ids), ravg);
end
save([cachedir cls '_boxes_' testset '_' suffix '.mat'], 'boxes', 'times');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function z = s2z4(S);

% Z = s2z4(S)
%
% Scattering to Impedance transformation
% only for N-by-4 matrix

%d = (1 - S(:, 1)).*(1 - S(:, 4)) - S(:,2).* S(:, 3);

for i = 1:size(S,1)
  d(i) = (1 - S(i, 1)).*(1 - S(i, 4)) - S(i,2).* S(i, 3);
  while abs(d(i)) < 1e-8
    fckindex = 1+round(rand*3);
    S(i, fckindex) = S(i, fckindex)*(1+rand*1e-8);
    d(i) = (1 - S(i, 1)).*(1 - S(i, 4)) - S(i,2).* S(i, 3);
  end;
end;

d = rot90(d, 3);

% at this point the I-S matrix should be non-singular

z(:, 1) = ((1 + S(:, 1)).*(1 - S(:, 4)) + S(:,2).* S(:, 3))./d;
z(:, 2) = 2* S(:,2)./d;
z(:, 3) = 2* S(:,3)./d;
z(:, 4) = ((1 - S(:, 1)).*(1 + S(:, 4)) + S(:,2).* S(:, 3))./d;

author: Sable


text:
% The COBRA Toolbox: testexportSetToGAMS
%
% Purpose:
%     - test exportSetToGAMS function
%
% Authors:
%     - Loic Marx, November 2018

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which(mfilename));
cd(fileDir);

% if nargin < 2
model = getDistributedModel('ecoli_core_model.mat');
fileName = 'refData_reactions.txt';

% read in the file
rxnsRef = textread(fileName, '%s', 'delimiter', '\n');

% test if the first and end cells are equal to /
assert(isequal(rxnsRef{1}, '/'))
assert(isequal(rxnsRef{end}, '/'))

% test if all the other cells correspond to the reaction names
for k = 2:length(model.rxns)-1
    assert(isequal(cellstr(rxnsRef{k}(2:end-1)), model.rxns(k-1)));
end

% run the function
exportSetToGAMS(model.rxns, 'Reactions.txt');

% compare to the ref data
refData = textread('Reactions.txt', '%s', 'delimiter', '\n');
assert(isequal(refData, rxnsRef))

% test for 1 input: exportSetToGAMS(model.rxns)
assert(verifyCobraFunctionError('exportSetToGAMS', 'inputs', {model.rxns}, 'testMessage', 'All inputs for the function exportSetToGAMS must be specified'))

% remove the test file
delete('Reactions.txt');

% change back to the current directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function [x,y]=limits(a)
% LIMITS returns min & max values of matrix; else scalar value.
%
%   [lo,hi]=LIMITS(a) returns LOw and HIgh values respectively.
%
%   lim=LIMITS(a) returns 1x2 result, where lim = [lo hi] values

% Copyright 2004-2010 The MathWorks, Inc.

if nargin~=1 | nargout>2 %bogus syntax
  error('usage: [lo,hi]=limits(a)')
end

siz=size(a);

if prod(siz)==1 %scalar
  result=a;                         % value
else %matrix
  result=[min(a(:)) max(a(:))];     % limits
end

if nargout==1 %composite result
  x=result;                         % 1x2 vector
elseif nargout==2 %separate results
  x=result(1);                      % two scalars
  y=result(2);
else %no result
  ans=result                        % display answer
end

repo: mcbench-benchmarks


function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,varargin)
% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
%                                       P1,P2,P3,P4,P5,P6,P7,P8)
% retcodes: 0, normal step.  5, largest step still improves too fast.
% 4,2 back and forth adjustment of stepsize didn't finish.  3, smallest
% stepsize still improves too slow.  6, no improvement found.  1, zero
% gradient.
%---------------------
% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
% Places where the number of P's need to be altered or the code could be returned to
% its old form are marked with ARGLIST comments.
%
% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
% update.
%
% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
% it's not psd.
%
%tailstr = ')';
%for i=nargin-6:-1:1
%   tailstr=[ ',P' num2str(i)  tailstr];
%end
%ANGLE = .03;
ANGLE = .005;
%THETA = .03;
THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
FCHANGE = 1000;
MINLAMB = 1e-9;
% fixed 7/15/94
% MINDX = .0001;
% MINDX = 1e-6;
MINDFAC = .01;
fcount=0;
lambda=1;
xhat=x0;
f=f0;
fhat=f0;
g = g0;
gnorm = norm(g);
%
if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
   retcode =1;
   dxnorm=0;
   % gradient convergence
else
   % with badg true, we don't try to match rate of improvement to directional
   % derivative.  We're satisfied just to get some improvement in f.
   %
   %if(badg)
   %   dx = -g*FCHANGE/(gnorm*gnorm);
   %  dxnorm = norm(dx);
   %  if dxnorm > 1e12
   %     disp('Bad, small gradient problem.')
   %     dx = dx*FCHANGE/dxnorm;
   %   end
   %else
   % Gauss-Newton step;
   %---------- Start of 7/19/93 mod ---------------
   %[v d] = eig(H0);
   %toc
   %d=max(1e-10,abs(diag(d)));
   %d=abs(diag(d));
   %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
%      toc
   dx = -H0*g;
%      toc
   dxnorm = norm(dx);
   if dxnorm > 1e12
      disp('Near-singular H problem.')
      dx = dx*FCHANGE/dxnorm;
   end
   dfhat = dx'*g0;
   %end
   %
   %
   if ~badg
      % test for alignment of dx with gradient and fix if necessary
      a = -dfhat/(gnorm*dxnorm);
      if a<ANGLE
         dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
         % suggested alternate code:  ---------------------
         dx = dx*dxnorm/norm(dx)    % This keeps scale invariant to the angle correction
         % ------------------------------------------------
         dfhat = dx'*g;
         % dxnorm = norm(dx);  % this line unnecessary with modification that keeps scale invariant
         disp(sprintf('Correct for low angle: %g',a))
      end
   end
   disp(sprintf('Predicted improvement: %18.9f',-dfhat/2))
   %
   % Have OK dx, now adjust length of step (lambda) until min and
   % max improvement rate criteria are met.
   done=0;
   factor=3;
   shrink=1;
   lambdaMin=0;
   lambdaMax=inf;
   lambdaPeak=0;
   fPeak=f0;
   lambdahat=0;
   while ~done
      if size(x0,2)>1
         dxtest=x0+dx'*lambda;
      else
         dxtest=x0+dx*lambda;
      end
      % home
      f = feval(fcn,dxtest,varargin{:});
      %ARGLIST
      %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
      % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
      disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f ))
      %debug
      %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda))
      if f<fhat
         fhat=f;
         xhat=dxtest;
         lambdahat = lambda;
      end
      fcount=fcount+1;
      shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ;
      growSignal = ~badg & ( (lambda > 0)  &  (f0-f > -(1-THETA)*dfhat*lambda) );
      if  shrinkSignal  &   ( (lambda>lambdaPeak) | (lambda<0) )
         if (lambda>0) & ((~shrink) | (lambda/factor <= lambdaPeak))
            shrink=1;
            factor=factor^.6;
            while lambda/factor <= lambdaPeak
               factor=factor^.6;
            end
            %if (abs(lambda)*(factor-1)*dxnorm < MINDX) | (abs(lambda)*(factor-1) < MINLAMB)
            if abs(factor-1)<MINDFAC
               if abs(lambda)<4
                  retcode=2;
               else
                  retcode=7;
               end
               done=1;
            end
         end
         if (lambda<lambdaMax) & (lambda>lambdaPeak)
            lambdaMax=lambda;
         end
         lambda=lambda/factor;
         if abs(lambda) < MINLAMB
            if (lambda > 0) & (f0 <= fhat)
               % try going against gradient, which may be inaccurate
               lambda = -lambda*factor^6
            else
               if lambda < 0
                  retcode = 6;
               else
                  retcode = 3;
               end
               done = 1;
            end
         end
      elseif  (growSignal & lambda>0) |  (shrinkSignal & ((lambda <= lambdaPeak) & (lambda>0)))
         if shrink
            shrink=0;
            factor = factor^.6;
            %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) | ( abs(lambda)*(factor-1)< MINLAMB)
            if abs(factor-1)<MINDFAC
               if abs(lambda)<4
                  retcode=4;
               else
                  retcode=7;
               end
               done=1;
            end
         end
         if ( f<fPeak ) & (lambda>0)
            fPeak=f;
            lambdaPeak=lambda;
            if lambdaMax<=lambdaPeak
               lambdaMax=lambdaPeak*factor*factor;
            end
         end
         lambda=lambda*factor;
         if abs(lambda) > 1e20;
            retcode = 5;
            done =1;
         end
      else
         done=1;
         if factor < 1.2
            retcode=7;
         else
            retcode=0;
         end
      end
   end
end
disp(sprintf('Norm of dx %10.5g', dxnorm))


text:
function test_issue922

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqanalysis

%test case for ft_freqanalysis rounding bug

load(dccnpath('/home/common/matlab/fieldtrip/data/test/issue922.mat'));
[freq] = ft_freqanalysis(cfg, test);

author: fieldtrip

repo: fieldtrip


text:
function [ a, b ] = p33_lim ( )

%*****************************************************************************80
%
%% P33_LIM returns the integration limits for problem 33.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the BEAM-EXAMPLE geometry and prints associated input files
%
% PYTHON INDEXING!!!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BeamCurve()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  32;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  32;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
N = 2*Nx;        % Number of Lagrangian Pts. (2x resolution of Eulerian grid)
a = 0.5;         % Length of beam (only in x-coordinate)
struct_name = 'BeamCurve'; % Name for .vertex, .spring, .beam, .target, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,Lx,a);

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis square;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .beam file! (TORSIONAL SPRINGS)
%k_Beam = 7.5e10; C = 0.0;
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .nonInv_beam file! (NON-INVARIANT))
k_Beam = 1e10;
print_Lagrangian_nonInv_Beams(xLag,yLag,k_Beam,struct_name)

% Prints .target file!
k_Target = 2e8;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', 2 );

    fprintf(target_fid, '%d %1.16e\n', 1-1, k_Target);
    fprintf(target_fid, '%d %1.16e\n', N-1, k_Target);

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (NON-INVARIANT) points to a file called rubberband.nonInv_beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_nonInv_Beams(xLag,yLag,k_Beam,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.nonInv_beam'], 'w');

    fprintf(beam_fid, '%d\n', N-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n',s-1-1, s-1, s+1-1, k_Beam, 0, 0);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n',s-1-1, s-1, 1-1,   k_Beam, 0, 0);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,Lx,a)

% The immsersed structure is a curved line %
ds = a/(N-1);

xLag(1) = -a/2;%Lx/2-a/2;
yLag(1) = 0.0;
for i=2:N

    xLag(i) = xLag(i-1) + ds;
    x = xLag(i);
    yLag(i) = -2*( (x)^2 - (a/2)^2 );

end

xLag = xLag + Lx/2;
yLag = yLag + Lx/2;

author: nickabattista


text:
function [windowsOut,scoreOut]=mvg_runObjectDetection(img,config)
% function [windowsOut,scoreOut]=mvg_runObjectDetection(img,config)
% returns a set of candidate windows which are prominent to contain
% objects over broad variety of classes.
%
% Inputs:
% img, imgRow*imgCol*3, double, is the input color image.
% config, struct, contains the parameter settings. See below for formatting
%                 details. If undefined, default configuration will be used.
%
% Outputs:
% windowsOut, numWindows*4, double, is a matrix containing the candidate
%                                   windows. Each row corresponds to one
%                                   window in format:
%                                   windowsOut(i,:)=[xmin,ymin,xmax,ymax];
% scoreOut, numWindows*1, double, is a vector of objectness scores associated
%                                 with each window in windowsOut. i:th entry
%                                 corresponds to i:th row in windowsOut.
%
% Example usage:
% >>img=imread('exampleImage.jpg');
% >>[windowsOut,scoreOut]=mvg_runObjectDetection(img);
% >>mvg_drawWindows(img,windowsOut((1:10),:));

% This program implements the method described in:
%
% Rahtu E. & Kannala J. & Blaschko M. B.
% Learning a Category Independent Object Detection Cascade.
% Proc. International Conference on Computer Vision (ICCV 2011).

% 2011 MVG, Oulu, Finland, Esa Rahtu and Juho Kannala
% 2011 VGG, Oxford, UK, Matthew Blaschko

%% Default configuration (for details, please refer to the paper above)
% The parameters you most likely need to change are config.NMS.numberOfOutputWindows that
% defines the number of output windows and config.NMS.trhNMSb, which defines the non-maxima
% suppression threshold (0.4-0.6 gives better recall with lower overlaps, 0.6-> gives better reall with high overlaps).
if ~exist('config','var') || isempty(config)
    %% Initial window sampling parameters
    % General
    config.InitialWindows.loadInitialWindows=false; % false->nothing is loaded, true->load initial windows from the file given in storage parameter below (overrides all other initial window settings)
    config.InitialWindows.windowTypes={'Prior','Superpix'}; % Methods for generating the initial windows
    config.InitialWindows.numInitialWindows=100000; % The number of initial windows returned
    % Superpixel windows
    config.InitialWindows.numberOfConnectedSuperpix=3; % The maximum number of connected superpixels used with the initial windows
    % Prior windows
    config.InitialWindows.windowPriorDistribution='ICCV_windowPriorDistribution.mat'; % mat-file or the structure containing the distribution of prior windows

    %% Feature parameters
    config.Features.loadFeatureScores=false; % false->nothing is loaded, true->load feature scores from the file given in storage parameter below (overrides all other initial feature settings)
    config.Features.featureTypes={'SS','WS','BE','BI'}; % Feature to be computed from initial windows
    config.Features.featureWeights=[1.685305e+00, 7.799898e-02, 3.020189e-01, -7.056292e-04]; % Relative feature weights (same ordering as with the features above)

    %% NMS parameters
    config.NMS.NMStype='NMSab'; % The type of Non-Maximum-Suppression used
    config.NMS.numberOfOutputWindows=1000; % Number of output windows retuned
    config.NMS.numberOfIntermediateWindows=10001; % Number of windows after NMSa algorithm
    config.NMS.trhNMSa=0; % Threshold for NMSa method
    config.NMS.trhNMSb=0.75; % Threshold for NMSb method

    %% Storage parameters (file names with full path, set empty if nothing will be stored or loaded)
    config.Storage.initialWindows=[];
    config.Storage.features=[];
    config.Storage.outputWindows=[];

    %% General parameters
    config.verbose=0; % 0->show nothing, 1->display progress
end

%% Initialize
% General
timeStamp=clock;
config.Features.verbose=config.verbose;
config.NMS.imageSize=[size(img,1),size(img,2)];

% Check loading paramets
if config.Features.loadFeatureScores && ~config.InitialWindows.loadInitialWindows
    error('If you load feature scores, also corresponding windows must be loaded. Otherwise windows and scores do not match!');
end

% Ensure superpixel representation exists if it is needed anywhere in the algorithm
if (sum(strcmp('Superpix',config.InitialWindows.windowTypes))>eps && ~config.InitialWindows.loadInitialWindows) || ((sum(strcmp('SS',config.Features.featureTypes))>eps || sum(strcmp('BI',config.Features.featureTypes))) && ~config.Features.loadFeatureScores)
    % Compute superpixels if they are not in config.superPixels
    if ~isfield(config,'superPixels')
        % Display
        if config.verbose>eps
            fprintf('\nComputing superpixels...\n');
            intermediateTimeStamp=clock;
        end
        % Compute superpixels
        config.superPixels=mvg_computeSuperpixels(img);
        % Display
        if config.verbose>eps
            fprintf('Superpixels computed! Time taken %1.2f sec.\n',etime(clock,intermediateTimeStamp));
        end
    end
    % Store superpixels for initial windows and feature scoring
    config.InitialWindows.superPixels=config.superPixels;
    config.Features.superPixels=config.superPixels;
end

%% Load or make initial windows
if config.InitialWindows.loadInitialWindows
    % Display
    if config.verbose>eps
        fprintf('\nLoading initial windows from:\n%s...\n',config.Storage.initialWindows);
    end

    % Load data
    load(config.Storage.initialWindows,'windows','configInitialWindows');

    % Reset initial window config
    config.InitialWindows=configInitialWindows;

else
    % Display
    if config.verbose>eps
        fprintf('\nCreating initial windows...\n');
        intermediateTimeStamp=clock;
    end
    % Make windows
    windows=mvg_makeInitialWindows(img,config.InitialWindows);
    % Display
    if config.verbose>eps
        fprintf('Initial windows done! Time taken %1.2f sec.\n',etime(clock,intermediateTimeStamp));
    end
    % Store initial windows
    if ~isempty(config.Storage.initialWindows)
        % Display
        if config.verbose>eps
            fprintf('Storing initial windows to:\n%s\n',config.Storage.initialWindows);
        end
        % Store initial windows
        configInitialWindows=config.InitialWindows;
        save(config.Storage.initialWindows,'windows','configInitialWindows');
    end
end

%% Load or compute feature scores
if config.Features.loadFeatureScores
    % Display
    if config.verbose>eps
        fprintf('\nLoading feature scores from:\n%s...\n',config.Storage.features);
    end

    % Load data
    load(config.Storage.features,'featureScores','configFeature');

    % Reset initial window config
    config.Features=configFeature;

else
    % Display
    if config.verbose>eps
        fprintf('\nComputing feature scores...\n');
        intermediateTimeStamp=clock;
    end
    % Run score computation
    featureScores=mvg_computeFeatureScores(img,windows,config.Features);
    % Display
    if config.verbose>eps
        fprintf('Features computed! Time taken %1.2f sec.\n',etime(clock,intermediateTimeStamp));
    end
    % Store feature scores
    if ~isempty(config.Storage.features)
        % Display
        if config.verbose>eps
            fprintf('Storing feature scores to:\n%s\n',config.Storage.features);
        end
        % Store feature scores
        configFeature=config.Features;
        save(config.Storage.features,'featureScores','configFeature');
    end
end

%% Make combined score
% Display
if config.verbose>eps
    fprintf('\nComputing combined scores...\n');
    intermediateTimeStamp=clock;
end
% Compute linear combination
combinedScores=(config.Features.featureWeights*featureScores')';
% Display
if config.verbose>eps
    fprintf('Combined score computed! Time taken %1.2f sec.\n',etime(clock,intermediateTimeStamp));
end

%% Run NMS to get indices to selected windows
% Display
if config.verbose>eps
    fprintf('\nSelecting final output windows...\n');
    intermediateTimeStamp=clock;
end
% Run NMS
selectedId=mvg_runWindowNMS(windows,combinedScores,config.NMS);
% Display
if config.verbose>eps
    fprintf('Final windows selected! Time taken %1.2f sec.\n',etime(clock,intermediateTimeStamp));
end

%% Take selected window coordinates and scores
windowsOut=windows(selectedId,:);
if nargout>1
    scoreOut=combinedScores(selectedId);
end
% Store output windows
if ~isempty(config.Storage.outputWindows)
    % Display
    if config.verbose>eps
        fprintf('Storing output windows to:\n%s\n',config.Storage.outputWindows);
    end
    % Store feauture scores
    save(config.Storage.outputWindows,'windowsOut','scoreOut','config');
end

%% Print done
if config.verbose>eps
    fprintf('Done!\n');
    fprintf('Total time taken: %1.2f sec.\n',etime(clock,timeStamp));
end

author: Cloud-CV

repo: object-proposals


text:
function fhandle = plotCylinder(VORT)

fhandle = figure

vortmin = -5;
vortmax = 5;
V2 = VORT;
% normalize values... not symmetric
minval = min(V2(:));
maxval = max(V2(:));
if(abs(minval)<5 && abs(maxval)<5)
    if(abs(minval)>abs(maxval))
        vortmax = maxval;
        vortmin = -maxval;
    else
        vortmin = minval;
        vortmax = -minval;
    end
end
V2(V2>vortmax) = vortmax;
V2(V2<vortmin) = vortmin;
imagesc(V2)
colormap(jet);
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 100 600 260])
axis equal
hold on

cvals = [-4 -2 -1 -.5 -.25 -.155];

[c1,h1] = contour(V2,cvals*vortmax/5,'--k','LineWidth',1.);
contour(V2,-cvals*vortmax/5,'-k','LineWidth',1.)

% Take all the info from the contourline output argument:
i0 = 1;
i2 = 1;
while i0 <  length(c1)
    i1 = i0+[1:c1(2,i0)];
    zLevel(i2) = c1(1,i0);
    hold on
    % And plot it with dashed lines:
    ph(i2) = plot(c1(1,i1),c1(2,i1),'k--','linewidth',1.);
    i0 = i1(end)+1;
    i2 = i2+1;
end
% Scrap the contourlines:
delete(h1)

t = (1:100)/100'*2*pi;
x = 49+25*sin(t);
y = 99+25*cos(t);
fill(x,y,[.3 .3 .3])
plot(x,y,'k','LineWidth',1.2)

set(gcf,'PaperPositionMode','auto')
set(gcf, 'Renderer', 'painters')

author: dynamicslab


text:
function [Et,Vt,Ft]=hexahedral_hexagon_beam(r,ne,hz,nz,XYZ_centre)

%%

%% CREATING 2D SECTION MESH

%Basic regular quad mesh
[X,Y]=meshgrid(linspace(0,1,ne+1));
Z=zeros(size(X));
V=[X(:) Y(:) Z(:)];
[F,V] = surf2patch(X,Y,Z);

%Creating 3 sheared quad meshes to construct hexagon

%V1
V1=V;
SH=eye(3,3); SH(1,2)=0; SH(2,1)=0.5;
V1=V1*SH;
S=eye(3,3); S(1,1)=1; S(2,2)=sqrt(3/4);
V1=V1*S;

%V2
a=pi/3;
R=[cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 0;];
V2=V1*R;

%V3
V3=V;
SH=eye(3,3); SH(1,2)=0; SH(2,1)=-0.5;
V3=V3*SH;
S=eye(3,3); S(1,1)=1; S(2,2)=sqrt(3/4);
V3=V3*S;
V3(:,1)=V3(:,1)+0.5; V3(:,2)=V3(:,2)+sqrt(3/4);

%Composing hexagon
Vu=[V1;V2;V3];
Fu=[F;F+size(V1,1);F+2.*size(V1,1)];

%Removing double points
[Fu,Vu,~,~,~,~]=unique_patch(Fu,Vu,[]);

%Scaling radius
[THETA,R] = cart2pol(Vu(:,1),Vu(:,2));
[Vu(:,1),Vu(:,2)] = pol2cart(THETA,r.*R);

%% CREATING 3D EXTRUDED MESH

F=Fu; V=Vu;

z_range=linspace(0,hz,nz+1);
Vt=repmat(V,numel(z_range),1);
Z_add=ones(size(V,1),1)*z_range; Z_add=Z_add(:);
Vt(:,3)=Vt(:,3)+Z_add;
Vt=Vt-ones(size(Vt,1),1)*mean(Vt,1); %centering around mean
Vt=Vt+ones(size(Vt,1),1)*XYZ_centre; %Translate to desired centre

Et=[];
for iz=1:1:numel(z_range)-1
    %       bottom                top
    Et=[Et; F+(size(V,1).*(iz-1)) F+(size(V,1).*iz)];
end

% f_order=[4 3 2 1 8 7 6 5];
% Et=Et(:,f_order);

[Ft]=hex2patch(Et);

end

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


text:
function slimgsetprep(srcfolder, dstpath, matsize, maxsec)
%SLIMGSETPREP organizes the images in a MATLAB friendly way
%
% $ Syntax $
%   - slimgsetprep(srcfolder, dstpath, matsize, maxsec)
%
% $ Argument $
%   - srcfolder:        the source file folder
%   - dstpath:          the destination file (without extension name)
%   - matsize:          the matrix size (row vector): [nrows, ncols]
%   - maxsec:           the maximum section size
%
% $ Description $
%   - slimgsetprep(srcfolder, dstpath, imgsiz, maxsec) prepares the matlab
%     files in destination for a set of images. The requirement is
%       1) In source folder: there are all the images and a DSDML
%          description file named dataset.xml
%       2) For destination, say abc, then there would be a core file named
%          abc.xml. And if the actual data is separately stored, there are
%          a series of files named abc.arr.<begin index>-<last index>. For
%          example, abc.arr.0001-2000, means that the 1st t0 2000th samples
%          are stored in the array file.
%       3) maxsec is the number of samples in each section. If maxsec is
%          specified, the images are stored in separate array files. If
%          maxsec is not specified or empty, the images are stored in core
%          file.
%       4) The core file is a MAT file with following variables:
%           'desc':         the DSDML data object descriptor
%           'sections':     the starting and ending indices of all sections
%           'matsize':      the row vector describing the size of images
%           'data':         the image array
%                           if imgs is numeric, it is the actual image
%                           array.
%                           if imgs is a cell array, it is the set of
%                           array filenames.
%
% $ History $
%   - Created by Dahua Lin, on Jul 26th, 2006
%

%% Parse and verify input arguments

if nargin < 3
    raise_lackinput('slimgsetprep', 3);
end

isseparate = false;
if nargin >= 4
    isseparate = ~isempty(maxsec);
end

matsize = matsize(:)';
if length(matsize) ~= 2
    error('sltoolbox:invalidarg', ...
        'imgsize should be a 2-element vector');
end

%% Read the dataset

descrfile = [srcfolder, '\dataset.xml'];
desc = dataset(descrfile);

N = desc.numsamples;
fns = desc.filenames;

%% Process data

if ~isseparate

    datasiz = [matsize, N];
    data = zeros(datasiz);

    for i = 1 : N
        img = imread([srcfolder, '\', fns{i}]);
        img = slimg2mat(img);
        data(:,:,i) = img;
    end

else

    sections = slpartition(N, 'maxblksize', maxsec);
    nsecs = length(sections.sinds);

    data = cell(nsecs, 1);

    % process filenames
    dstdir = fileparts(dstpath);
    numlen = length(num2str(N));
    if isempty(dstdir)  % local dst
        dstfnpat = [dstpath, sprintf('.arr.%%0%dd-%%0%dd', numlen, numlen)];
        dstfppre = [];
    else
        dstfnpat = [dstpath(length(dstdir)+2:end), sprintf('.arr.%%0%dd-%%0%dd', numlen, numlen)];
        dstfppre = [dstdir, '\'];
    end

    for i = 1 : nsecs
        si = sections.sinds(i);
        ei = sections.einds(i);

        curfn = sprintf(dstfnpat, si, ei);
        curfp = [dstfppre, curfn];
        data{i} = curfn;

        curN = ei - si + 1;
        arr = zeros([matsize, curN]);

        for j = 1 : curN
            curidx = si + j - 1;
            img = imread([srcfolder, '\', fns{curidx}]);
            img = slimg2mat(img);
            arr(:,:,j) = img;
        end

        slwritearray(arr, curfp);
    end

end

%% Output core

corefile = [dstpath, '.mat'];
save(corefile, 'desc', 'sections', 'matsize', 'data', '-v6');

author: lmthang


text:
function quality_test_halton ( )

%*****************************************************************************80
%
%% QUALITY_TEST_HALTON calls the QUALITY routines for points in the unit hypercube.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUALITY_TEST_HALTON\n' );
  fprintf ( 1, '  Test the MATLAB QUALITY library\n' );
  fprintf ( 1, '  on points in the unit hypercube.' );

  ns = 100000;
  seed_init = 123456789;
  input_filename = 'halton_02_00100.txt';

  [ dim_num, n ] = r8mat_header_read ( input_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The pointset was read from "%s"\n', input_filename );
  fprintf ( 1, '  The sample routine is SAMPLE_HYPERCUBE_UNIFORM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM =  %12d\n', dim_num   );
  fprintf ( 1, '  The number of points N =	   %12d\n', n         );
  fprintf ( 1, '  The number of sample points NS = %12d\n', ns        );
  fprintf ( 1, '  The random number SEED_INIT =    %12d\n', seed_init );
  fprintf ( 1, '\n' );

  z = r8mat_data_read ( input_filename, dim_num, n );

  r8mat_transpose_print_some ( dim_num, n, z, 1, 1, 5, 5, ...
    '  5x5 portion of data read from file:' );
%
%  For 2 dimensional datasets, some quality measures work from the Delaunay triangulation.
%  Compute that here.
%
  if ( dim_num == 2 )

    triangle_node = delaunay ( z(1,1:n), z(2,1:n) );

    triangle_node = triangle_node';
    [ dummy, triangle_num ] = size ( triangle_node );

  else

    triangle_node = [];
    triangle_num = 0;

  end

  if ( dim_num == 2 )
    quality_test005 ( n, z, triangle_num, triangle_node );
    quality_test006 ( n, z, triangle_num, triangle_node );
  end
  quality_test007 ( dim_num, n, z );
  quality_test01 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test02 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test03 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test04 ( dim_num, n, z );
  quality_test05 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test06 ( dim_num, n, z );
  quality_test07 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test08 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  if ( dim_num == 2 )
    quality_test083 ( n, z, triangle_num, triangle_node );
  end
  quality_test085 ( dim_num, n, z );
  quality_test09 ( dim_num, n, z );
  quality_test10 ( dim_num, n, z, ns, 'sample_hypercube_uniform', seed_init );
  quality_test11 ( dim_num, n, z );

  return
end


text:
function SE = functionComputeSE_DL_poweralloc(rho,signal,interference,prelogFactor)
%Compute the SE in Theorem 4.6 using the formulation in (7.1) for a given
%power allocation scheme.
%
%INPUT:
%rho          = K x L matrix where element (k,j) is the downlink transmit
%               power allocated to UE k in cell j
%signal       = K x L matrix where element (k,j,n) is a_jk in (7.2)
%interference = K x L x K x L matrix where (l,i,jk,n) is b_lijk in (7.3)
%prelogFactor = Prelog factor
%
%OUTPUT:
%SE = K x L matrix where element (k,j) is the downlink SE of UE k in cell j
%     using the power allocation given as input
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Extract number of UEs
K = size(signal,1);

%Extract number of cells
L = size(signal,2);

%Prepare to save results
SE = zeros(K,L);

%% Go through all cells
for j = 1:L

    %Go through all UEs in cell j
    for k = 1:K

        %Compute the SE in Theorem 4.6 using the formulation in (7.1)
        SE(k,j) = prelogFactor*log2(1+(rho(k,j)*signal(k,j)) / (sum(sum(rho.*interference(:,:,k,j))) + 1));

    end

end

author: emilbjornson

repo: massivemimobook


text:
function y = isPalindrome(inputWord)
    % Return 1 if input string is a palindrome; 0 otherwise.

    % Gets the length of the input word.
    inputLength = length(inputWord);

    % Number of pairs to compare.
    lastIndex = floor( inputLength / 2);
    y=1;
    % for every pair of letters in the word check if they match
    for i = 1 : lastIndex
        if inputWord(i) ~= inputWord(end + (1 - i))
            y = 0;
        end
    end

end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function b = r8but_to_r8ge ( n, mu, a )

%*****************************************************************************80
%
%% R8BUT_TO_R8GE copies a R8BUT matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8BUT storage format is for a banded upper triangular matrix.
%
%    To save storage, only the diagonal and upper triangle of A is stored,
%    in a compact diagonal format that preserves columns.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrices.
%    N must be positive.
%
%    Input, integer MU, the upper bandwidth of A1.
%    MU must be nonnegative, and no greater than N-1.
%
%    Input, real A(MU+1,N), the R8BUT matrix.
%
%    Output, real B(N,N), the R8GE matrix.
%
  for i = 1 : n
    for j = 1 : n
      if ( i <= j & j <= i+mu )
        b(i,j) = a(mu+1+i-j,j);
      else
        b(i,j) = 0.0;
      end
    end
  end

  return
end


text:
function mu = kseeds(X, k)
% Perform kmeans++ seeding
% Input:
%   X: d x n data matrix
%   k: number of seeds
% Output:
%   mu: d x k seeds
% Written by Mo Chen (sth4nth@gmail.com).
n = size(X,2);
D = inf(1,n);
mu = X(:,ceil(n*rand));
for i = 2:k
    D = min(D,sum((X-mu(:,i-1)).^2,1));
    mu(:,i) = X(:,randp(D));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function inv = cauchy_inv(x,m,s)
%---------------------------------------------------
% PURPOSE:
% Estimates the Inverse of the Cumulative Distribution Function of the
% Cauchy-Lorentz Distribution for a series of x values, with m location
% parameter and s scale parameter
%---------------------------------------------------
% USAGE:
% inv = cauchy_inv(x,m,s)
% where:
% x = (n x 1) or (1 x n) vector
% m = (n x 1) or (1 x n) vector
% s = (n x 1) or (1 x n) vector
%---------------------------------------------------
% RETURNS:
% cdf = (n x 1) or (1 x n) vector containing the probability for each
% element of x with corresponding location and scale parameters
%---------------------------------------------------
% Author:
% Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
% Date: 06/2010
%---------------------------------------------------

if nargin == 0
    error('Data, Location Parameter, Scale Parameter')
end

inv = m + s.*tan(pi*(x-1/2));

end


text:
function sl_customization(cm)

% Copyright 1990-2013 The MathWorks, Inc.
% $Revision: 1.1.8.1 $

    cm.addCustomMenuFcn('Simulink:ContextMenu',@getMyMenuItems);
    cm.addCustomFilterFcn('sig_call_back',@myFilter);
end

%% Define the custom menu function.
function schemaFcns = getMyMenuItems(callbackInfo)
  schemaFcns = {@getItem1};
end

%% create filter for menu to appear only for Simscape Probe blocks
function state = myFilter(callbackInfo)
  if strmatch(get_param(gcb,'maskdescription'),'Simscape Probe','exact')
    state = 'Enabled';
  else
    state = 'Hidden';
  end
end

%% Define the schema function for first menu item.
function schema = getItem1(callbackInfo)
  schema = sl_container_schema;
  if strmatch(get_param(gcb,'maskdescription'),'Simscape Probe','exact')
    schema.state = 'Enabled';
    schema.childrenFcns = {@units, @outPort @tagType};
    schema.label = 'Configure Probe';
  else
    schema.state = 'Hidden';
  end
end
%% Define Sub-menu entries
function schema = units(callbackInfo)
  schema = sl_action_schema;
  schema.label = 'Set Units';
  schema.callback = @setUnits;
end

function schema = outPort(callbackInfo)
  schema = sl_action_schema;
  schema.label = 'Add Output Port';
  schema.callback = @addOutputPort;
end

function schema = tagType(callbackInfo)
  schema = sl_toggle_schema;
  schema.label = 'Global';
  if strcmp(get_param([gcb, '/Goto'], 'TagVisibility'),'global') == 1
    schema.checked = 'checked';
  else
    schema.checked = 'unchecked';
  end
  schema.callback = @setTag;
end

%% Function callbacks from menu options
function setUnits(callbackInfo)

open_system([gcb, '/PS-S']) % Open PS-S dialog box to set units

end

function addOutputPort(callbackInfo)
p=get_param([gcb, '/PS-S'],'Position'); % get position of PS-S block as a reference
hOutblock = add_block('built-in/Outport', [gcb '/Out'],...
            'MakeNameUnique', 'on','Position', ...
            [p(1)+150 p(2)+100 p(3)+175 p(4)+100]);
add_line(gcb, ['PS-S', '/1'], ['Out', '/1'],'autorouting','on');
end

function setTag(callbackInfo)
tv=get_param([gcb, '/Goto'],'TagVisibility');
  if  regexp(tv,'global') == 1
     set_param([gcb, '/Goto'],'TagVisibility', 'local');
  elseif regexp(tv,'local') == 1
     set_param([gcb, '/Goto'],'TagVisibility', 'global');
  end
end


text:
function [Phi,S,Lambda,M] = extract_eigen_functions_new(shape,k)

    [M,S]=laplacian([shape.X shape.Y shape.Z],shape.TRIV);
    M = diag(sum(M,2));
    [Phi,Lambda] = eigs(-S,M,k,1e-5);
    Lambda = diag(Lambda);
    [Lambda,idx] = sort(Lambda,'descend');
    Phi = Phi(:,idx);
%     PhiI = Phi'*M;
    Lambda = abs(Lambda); %added this to return positive eigen values.

end

author: OshriHalimi


text:
function f = feature_extract( imdist ,scale)
    f=[];
    imdist=rgb2gray(imdist);
    weight=[0.2 0.8];
    for i=1:scale
        im=imdist;
        fun0=@(x)secal(x);
        emat=blkproc(im,[8 8] ,fun0);

        sort_t = sort(emat(:),'ascend');
        len = length(sort_t);
        t=sort_t(ceil(len*weight(1)):ceil(len*weight(2)));
        mu= mean(t);
        ske=skewness(sort_t);

        f1=[ mu  ske];

        im=imdist;
        fun1=@(x)fecal(x);
        im=double(im);
        femat=blkproc(im,[8 8],fun1);

        sort_t = sort(femat(:),'ascend');
        len = length(sort_t);
        t=sort_t(ceil(len*weight(1)):ceil(len*weight(2)));
        mu= mean(t);
        ske=skewness(sort_t);

        f2=[ mu  ske];

        f=[f f1 f2] ;

        imdist = imresize(imdist,0.5);
    end
end

author: thomas-koehler

repo: SupER


text:
function [Ka, pKa] = acidDissociationConstant(metAbbr, Alberty2006, metAbbrAlbertyAbbr, temp, is, chi)
% Acid dissociation constant for the different metabolite species that make up a reactant
%
% USAGE:
%
%    [Ka, pKa] = acidDissociationConstant(metAbbr, Alberty2006, metAbbrAlbertyAbbr, temp, is, chi)
%
% INPUTS:
%    metAbbr:               reconstruction reactant abbreviation
%    Alberty2006:           Basic data on the metabolite species that make
%                           up a reactant, compiled by Robert A. Alberty,
%                           Massachusetts Institute of Technology.
%                           In Print: `Robert A. Alberty, Biochemical Thermodynamics:
%                           Applications of Mathematica. John Wiley & Sons, 2006. p391-395`
%                           Online: BasicBioChemData3.nb
%                           http://library.wolfram.com/infocenter/MathSource/5704/
%    metAbbrAlbertyAbbr:    mapping from model metabolite primary key to
%                           primary key of reactants in `Alberty2006`
%
% OPTIONAL INPUT
%    temp:                  temperature (default 298.15 K)
%    is:                    ionic strength (default 0 M)
%    chi:                   electrical potential (default 0)
%
% OUTPUTS:
%    Ka:                    apparent equilibrium constants
%    pKa:                   :math:`-log_{10}(Ka)`

if ~exist('temp','var')
    temp=298.15;
end
if ~exist('is','var')
    is=0;
end
if ~exist('chi','var')
    chi=0;
end

%find the alberty abbreviation for this metabolite Abbreviation
albertyAbbr=metAbbrAlbertyAbbr(strcmp(metAbbr,metAbbrAlbertyAbbr(:,2)),3);

%make a list of alberty abbreviations
klt=size(Alberty2006,2);
allAlbertyAbbr=cell(klt,1);
for k=1:klt
    allAlbertyAbbr{k}=Alberty2006(k).abbreviation;
end
%index for matching data for the alberty abbreviation
n=find(strcmp(albertyAbbr,allAlbertyAbbr));

%find the number of species within pseudoisomer group
p=max(find(~isnan(Alberty2006(n).basicData(:,1))));

%no Legendre transformation for pH or electrical potential
Legendre     = 0;
LegendreCHI  = 0;
pHr = 7; %dummy, has no effect

[dGf0,dHf0,mf,aveHbound,aveZi,lambda,gpfnsp]=calcdGHT(Alberty2006(n).basicData(1:p,1),Alberty2006(n).basicData(1:p,2),Alberty2006(n).basicData(1:p,3),Alberty2006(n).basicData(1:p,4),pHr,is,temp,chi,Legendre,LegendreCHI);

Ka=zeros(p-1);
pKa=zeros(p-1);

gasConstant = 8.314472/1000; % kJ K-1 mol-1

%RTalpha p 49 Alberty 2003
%where alpha is the Debye-Huckel Constant
gibbscoeff = (9.20483*temp)/10^3 - (1.284668*temp^2)/10^5 + (4.95199*temp^3)/10^8;

% A <--> A- + H+
% Take into account the ionic strength effect on the Gibbs energy of a
% hydrogen ion
zi=1;
hydrogenIonistermG = -(gibbscoeff*(zi.^2)*is^0.5)/(1 + 1.6*is^0.5);

for x=1:p-1

    Ka(x)  =  exp((gpfnsp(x+1,1)-gpfnsp(x,1)-hydrogenIonistermG)/(gasConstant*temp));
    if 1
        pKa(x) =  -(gpfnsp(x+1,1)-gpfnsp(x,1)-hydrogenIonistermG)/(gasConstant*temp*log(10));
    else
        %p244 Alberty 2003 uses a slightly less accurate gas constant so slight
        %difference
        gpfnsp=Alberty2006(n).basicData(1:p,1);
        pKa(x) =  -(gpfnsp(x+1,1)-gpfnsp(x,1))/(log(10)* 8.31451 * .29815);
    end
end


text:
function d = Mtv(solver,I,i)
% forward operator
%
% Input
%  I - intensities
%  alpha - scaling factors
%
% Output
%  d - result
%

% extend specimen symmetry
lss = numProper(solver.SS);
I_ext = repmat(I.',lss,1)./lss;

% compute Fourier coefficients
nfsftmex('set_f', solver.nfft_r(i), I_ext);
nfsftmex('adjoint', solver.nfft_r(i));
fhat = nfsftmex('get_f_hat_linear', solver.nfft_r(i));

% convolution with kernel function
fhat = 4*pi*fhat .* solver.A;

% evaluate Fourier series at pole figure points g h_i
nfsftmex('set_f_hat_linear', solver.nfft_gh(i), fhat);
nfsftmex('trafo', solver.nfft_gh(i));
d = real(nfsftmex('get_f', solver.nfft_gh(i)));
d = reshape(d,length(solver.c),[]) * solver.refl{i}.';

end

author: mtex-toolbox

repo: mtex

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