text:
% test_itk_imfilter.m
%
% Script to test the filters provided by itk_imfilter.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011-2014 University of Oxford
% Version: 0.4.5
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%% itk::BinaryThinningImageFilter3D (skel) filter

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd.data = nrrd.data > 150;

% plot data
close all
imagesc(nrrd.data(:, :, 3))

% skeletonise (0.04 s)
tic
im = itk_imfilter('skel', nrrd);
toc

figure
imagesc(im(:, :, 3))

%% itk::DanielssonDistanceMapImageFilter (dandist) filter

% very simple 3x3 matrix with only central pixel = 1
im = zeros(3);
im(2, 2) = 1;
nrrd = scimat_im2scimat(im, [.5 .5 .5], [0 0 0]);
[im2, v, w] = itk_imfilter('dandist', nrrd);

% vector is given in voxel coordinates as (r, c, s) instead of (x, y, z)
w()

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd2 = nrrd;
nrrd2.data(:) = 0;
nrrd2.data(:, 20, :) = 1;
nrrd2.data(:, 45, :) = 2;

% plot data
close all
imagesc(nrrd2.data(:, :, 3))

% compute distance (0.11 s)
tic
[im, v, w] = itk_imfilter('dandist', nrrd2);
toc

% distance map
figure
imagesc(im(:, :, 3))

% voronoi diagram
figure
imagesc(v(:, :, 3))

%% itk::SignedMaurerDistanceMapImageFilter (maudist) filter

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd2 = nrrd;
nrrd2.data(:) = 0;
nrrd2.data(:, 20:23, :) = 1;
nrrd2.data = uint16(nrrd2.data);

% plot data
close all
imagesc(nrrd2.data(:, :, 3))

% compute distance (0.11 s)
tic
im = itk_imfilter('maudist', nrrd2.data);
toc

% distance map
figure
imagesc(im(:, :, 3))

%% itk:::BinaryDilateImageFilter (bwdilate) filter

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd.data = nrrd.data > 150;

% plot data
close all
imagesc(nrrd.data(:, :, 3))

% dilate (0.02 s)
tic
im = itk_imfilter('bwdilate', nrrd, 4);
toc

figure
imagesc(im(:, :, 3))

%% itk:::BinaryErodeImageFilter (bwerode) filter

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd.data = nrrd.data > 150;

% plot data
close all
imagesc(nrrd.data(:, :, 3))

% erode (0.02 s)
tic
im = itk_imfilter('bwerode', nrrd, 1);
toc

figure
imagesc(im(:, :, 3))

%% itk::AnisotropicDiffusionVesselEnhancementImageFilter (advess) filter

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd.data = single(nrrd.data);

% plot data
close all
imagesc(nrrd.data(:, :, 3))

% user-provided parameters
sigmaMin = nrrd.axis(1).spacing * 5;
sigmaMax = nrrd.axis(1).spacing * 10;
sigmaSteps = 15;
isSigmaStepLog = false;
iterations = 30;
wStrength = 24;
sensitivity = 4.0;
timeStep = 1e-3;
epsilon = 1e-2;

% smooth along vessels (11.4 s)
tic
im = itk_imfilter('advess', nrrd, sigmaMin, sigmaMax, sigmaSteps, ...
    isSigmaStepLog, iterations, wStrength, sensitivity, timeStep, epsilon);
toc

figure
imagesc(im(:, :, 3))

figure
imagesc(nrrd.data(:, 2:end, 3) - im(:, 2:end, 3))

%% itk::MultiScaleHessianSmoothed3DToVesselnessMeasureImageFilter (hesves)

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');
nrrd.data = single(nrrd.data);

% plot data
close all
imagesc(nrrd.data(:, :, 5))

% user-provided parameters
sigmaMin = 1;
sigmaMax = 8;
sigmaSteps = 4;
isSigmaStepLog = false;

% compute vesselness measure
tic
im = itk_imfilter('hesves', nrrd.data, sigmaMin, sigmaMax, sigmaSteps, ...
    isSigmaStepLog);
toc

figure
imagesc(im(:, :, 5))

%% itk::MedianImageFilter (median)

% load test data
nrrd = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');

% default median filtering (no filtering)
im2 = itk_imfilter('median', nrrd);

err = im2 - nrrd.data;
% expected result = 0
any(err(:) ~= 0)

% median filtering only along rows, line of length 5+1+5=11
im2 = itk_imfilter('median', nrrd, [5 0 0]);

% plot difference
subplot(2, 1, 1)
imagesc(nrrd.data(:, :, 4))
subplot(2, 1, 2)
imagesc(im2(:, :, 4))

% median filtering only along columns, line of length 20+1+20=41
im2 = itk_imfilter('median', nrrd, [0, 20, 0]);

% plot difference
subplot(2, 1, 1)
imagesc(nrrd.data(:, :, 4))
subplot(2, 1, 2)
imagesc(im2(:, :, 4))

% compare speed and results with Matlab's implementation
tic
im2 = itk_imfilter('median', nrrd, [10, 10, 10]); % 5.9 sec
toc
tic
im3 = medfilt3(nrrd.data, [21, 21, 21]); % 36.5 sec
toc

% difference: expected result = 0
any(double(im2(:)) - double(im3(:)))

% plot both results
subplot(2, 1, 1)
imagesc(im2(:, :, 4))
subplot(2, 1, 2)
imagesc(im3(:, :, 4))

%% itk::MRFImageFilter (mrf)

% load test data
scimat = scimat_load('../../cpp/src/third-party/IJ-Vessel_Enhancement_Diffusion.1/CroppedWholeLungCTScan.mhd');

% plot image
hold off
subplot(2, 2, 1)
imagesc(scimat.data(:,:,4))

% initial segmentation
[thr, q, obj, seg] = gmthr_seg(double(scimat.data), 2);

% quality of the separation between classes (expected: q = 0.9020)
q

% plot segmented image using the Gaussian mixture model
subplot(2, 2, 2)
imagesc(seg(:,:,4))

% compute neighbourhood weights that decrease with Euclidean distance
[gr, gc, gs] = ndgrid(-3:3, -3:3, -2:2);
weights = 1./sqrt(gr.^2 + gc.^2 + gs.^2);
weights(isinf(weights)) = 0.0;

% segment the image using the Markov Random Filter algorithm
seg = itk_imfilter('mrf', scimat.data, obj.mu', weights);

% plot segmented image using the Gaussian mixture model
subplot(2, 2, 3)
imagesc(seg(:,:,4))

% repeat segmentation, now with significant smoothing
seg = itk_imfilter('mrf', scimat.data, obj.mu', weights, 2);

% plot segmented image using the Gaussian mixture model
subplot(2, 2, 4)
imagesc(seg(:,:,4))

%% itk::VotingBinaryIterativeHoleFillingImageFilter (voteholefill)

% create small toy image with a small hole
im = zeros(15, 15, 'uint8');
im(3:13, 3:13) = 1;
im(7:8, 7) = 0;

% plot input image
hold off
imagesc(im)

% run filter to fill in holes
RAD = [2 2];
THR = 2;
im2 = itk_imfilter('voteholefill', im, RAD, THR);

% plot output image
imagesc(im2)

%% itk::CannyEdgeDetectionImageFilter

% load rat heart MRI
scimat = scimat_load('data/008-mri-downsampled-4.mha');

% normalise the spacing, otherwise the filter produces NaN voxels
inc = min([scimat.axis.spacing]);
scimat.axis(1).spacing = scimat.axis(1).spacing / inc;
scimat.axis(2).spacing = scimat.axis(2).spacing / inc;
scimat.axis(3).spacing = scimat.axis(3).spacing / inc;

% plot one of the intermediate slices
subplot(2, 2, 1)
hold off
imagesc(scimat.data(:, :, 112));
colormap(gray)

% filter parameters

% standard deviation of the Gaussian filtering is roughly 5 times the
% the voxel size in each dimension
sigma = [scimat.axis.spacing] * 2;

% upper threshold
uppthr = 450;

% lower threshold
lowthr = uppthr/2;

% maximum error
maxerr = 0.01*ones(1, 3);

% run Canny edge filter
[im, im2] = itk_imfilter('canny', scimat, sigma.^2, uppthr, lowthr, maxerr);

% plot edge detection result
subplot(2, 2, 2)
imagesc(im2(:, :, 112))
subplot(2, 2, 3)
imagesc(im(:, :, 112))

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function dist = distance(x1,y1,x2,y2)
%This function calculates the distance between any two cartesian
%coordinates.
%   Copyright 2009-2010 The MathWorks, Inc.
dist=sqrt((x1-x2)^2 + (y1-y2)^2);

author: Mesywang

repo: Motion-Planning-Algorithms


text:
clear all;
clc;
clf;
% coordinates and connectivities
 x=[0 2 2 0];
 y=[0 0.5 1 1];

a = 3;
%mesh number in one side
n=2^(a-1);%(nxn)
% numberElements: number of elements
numberElements=n^2;
% numberNodes: number of nodes
numberNodes=(n+1)^2;
x1=linspace(x(1),x(2),n+1);
x2=linspace(x(4),x(3),n+1);
y1=linspace(y(1),y(4),n+1);
y2=linspace(y(2),y(3),n+1);
tempY=zeros(n+1,n+1);
tempX=zeros(n+1,n+1);
for i=1:n+1
    tempY(i,:)=linspace(y1(i),y2(i),n+1);
    tempX(:,i)=linspace(x1(i),x2(i),n+1);
end
ct=1;
nodeCoordinates=zeros(numberNodes,2);
for i=1:n+1
    for j=1:n+1
        nodeCoordinates(ct,:)=[tempX(i,j) tempY(i,j)];
        ct=ct+1;
    end
end

elementNodes=zeros(numberElements,4);
for i=1:n
    for j=1:n
        elementNodes(j+n*(i-1),:)=...
            [j+(n+1)*(i-1)...
            j+1+(n+1)*(i-1)...
            j+n+2+(n+1)*(i-1)...
            j+n+1+(n+1)*(i-1)];
    end
end
drawingMesh(nodeCoordinates,elementNodes,'Q4','k-');

repo: Aero-matlab


text:
function S=int2rgbm(R,G,B)
%INT2RGBM Convert intensity arrays to RGB movie.
%   INT2RGBM(R,G,B) converts separate 3-D intensity movies R, G,
%   and B, representing the red, green, and blue color planes of
%   an RGB video, respectively, into a 4-D RGB movie array.

% Copyright 2003 The MathWorks, Inc.
% $Revision: 1.1 $ $Date: 2004/12/17 02:43:25 $

error(nargchk(3,3,nargin));
sz=size(R);
sz=[sz(1:2) 1 sz(3)];
S = cat(3,reshape(R,sz),reshape(G,sz),reshape(B,sz));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of an example serial arm
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   20/11/2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.name= 'Serial example';

robot.DH.theta= '[ ]';
robot.DH.d='[]';
robot.DH.a='[]';
robot.DH.alpha= '[]';

robot.J=[];

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

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

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

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-165) deg2rad(165); %Axis 1, minimum, maximum
                deg2rad(-110) deg2rad(110); %Axis 2, minimum, maximum
                deg2rad(-110) deg2rad(70); %Axis 3
                deg2rad(-160) deg2rad(160); %Axis 4:
                deg2rad(-120) deg2rad(120); %Axis 5
                deg2rad(-400) deg2rad(400)]; %Axis 6:

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

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

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

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

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

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text:
function tr_xform=cosmo_fmri_convert_xform(fmt, xform)
% convert xform code between numeric and string in fmri dataset
%
% tr_xform=cosmo_fmri_convert_xform(fmt, xform)
%
% Inputs:
%   fmt              one of:
%                    - 'nii' (nifti)
%                    - 'afni' (afni)
%                    - 'bv' (BrainVoyager)
%   xform            either
%                    - a string, in which it must be one of
%                       'scanner_anat','aligned_anat','talairach',
%                       'mni_152', 'unknown'
%                    - a number, in which case it must be the view type for
%                       nii (sform or qform) or afni (the first elment of
%                       SCENE_TYPE). Any number in bv corresponds with
%                       'unknown'.
%
% Output:
%   tr_xform         if the input is a string, then the output is the
%                    corresponding view type (numeric); and vice versa
%                    if the viewtype cannot be found, the string 'unknown'
%                    or the corresponding view type is returned
%
% Examples
%     % Number to string:
%
%     % - NIFTI
%     cosmo_fmri_convert_xform('nii',1)
%     %|| 'scanner_anat'
%     cosmo_fmri_convert_xform('nii',2)
%     %|| 'aligned_anat'
%     cosmo_fmri_convert_xform('nii',3)
%     %|| 'talairach'
%     cosmo_fmri_convert_xform('nii',4)
%     %|| 'mni_152'
%     % any unknown number returns 'scanner_anat'
%     cosmo_fmri_convert_xform('nii',0)
%     %|| 'scanner_anat'
%     cosmo_fmri_convert_xform('nii',NaN)
%     %|| 'scanner_anat'
%
%     % - AFNI
%     cosmo_fmri_convert_xform('afni',0)
%     %|| 'scanner_anat'
%     cosmo_fmri_convert_xform('afni',1)
%     %|| 'aligned_anat'
%     cosmo_fmri_convert_xform('afni',2)
%     %|| 'talairach'
%     % any unknown number returns 'scanner_anat'
%     cosmo_fmri_convert_xform('afni',10)
%     %|| 'scanner_anat'
%     cosmo_fmri_convert_xform('afni',NaN)
%     %|| 'scanner_anat'
%
%     % - BV
%     % any number returns 'talairach'
%     cosmo_fmri_convert_xform('bv',0)
%     %|| 'talairach'
%     cosmo_fmri_convert_xform('bv',NaN)
%     %|| 'talairach'
%
%     % String to number:
%
%     % - NIFTI
%     cosmo_fmri_convert_xform('nii','scanner_anat')
%     %|| 1
%     cosmo_fmri_convert_xform('nii','aligned_anat')
%     %|| 2
%     cosmo_fmri_convert_xform('nii','talairach')
%     %|| 3
%     cosmo_fmri_convert_xform('nii','mni_152')
%     %|| 4
%     cosmo_fmri_convert_xform('nii','unknown')
%     %|| 1
%
%     % - AFNI
%     cosmo_fmri_convert_xform('afni','scanner_anat')
%     %|| 0
%     cosmo_fmri_convert_xform('afni','aligned_anat')
%     %|| 1
%     cosmo_fmri_convert_xform('afni','talairach')
%     %|| 2
%     % treat as talairach
%     cosmo_fmri_convert_xform('afni','mni_152')
%     %|| 2
%     % unkown
%     cosmo_fmri_convert_xform('afni','unknown')
%     %|| 0
%
%     % - BV
%     % all are unknown in BV, because BV does not support this
%     cosmo_fmri_convert_xform('bv','scanner_anat')
%     %|| 0
%     cosmo_fmri_convert_xform('bv','aligned_anat')
%     %|| 0
%     cosmo_fmri_convert_xform('bv','talairach')
%     %|| 0
%     cosmo_fmri_convert_xform('bv','mni_152')
%     %|| 0
%     cosmo_fmri_convert_xform('bv','unknown')
%     %|| 0
%
%     % any other input gives an error
%     cosmo_fmri_convert_xform('unknown','talairach')
%     %|| error('unsupported format ''unknown''')
%     cosmo_fmri_convert_xform('unknown',1)
%     %|| error('unsupported format ''unknown''')
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [x,y]=get_mapping(fmt);

    if isnumeric(xform)
        % convert numeric to string
        i=find(y==xform,1,'first');
        if isempty(i)
            tr_xform=get_default(x,y);
        else
            tr_xform=x{i};
        end

    elseif ischar(xform)
        % convert string to numeric
        i=find(cosmo_match(x,xform),1,'first');
        if isempty(i)
            [unused,tr_xform]=get_default(x,y);
        else
            assert(numel(i)==1);
            tr_xform=y(i);
        end
    else
        error('unsupported input');
    end

function [xd,yd]=get_default(x,y)
    i=find(isnan(y));
    assert(numel(i)==1);
    xd=x{i};

    msk=cosmo_match(x,xd);
    msk(i)=false;

    j=find(msk);
    assert(numel(j)==1);
    yd=y(j);

function [x,y]=get_mapping(fmt)
    if ~ischar(fmt)
        error('first argument must be a string');
    end

    % mapping from codes to indices
    switch fmt
        case 'nii'
            codes={'scanner_anat',1;...
                    'aligned_anat',2;...
                    'talairach',3;...
                    'mni_152',4;...
                    'scanner_anat',NaN}; % default

        case 'afni'
            codes={'scanner_anat',0;...
                   'aligned_anat',1;... % ACPC
                   'talairach',2;...
                   'mni_152',2;...
                   'scanner_anat',NaN}; % default

        case 'bv'
            codes={'talairach',0;...
                    'talairach',NaN}; %

        otherwise
            error('unsupported format ''%s''', fmt);
    end

    x=codes(:,1);
    y=cell2mat(codes(:,2));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function y = HighPassModulationFilter(x, para)

[nFr, dim] = size(x);

switch para.filter_type
    case 'CMN'
        y = CMN(x);

    case 'CMN_online'
        y = CMN_online(x);

    case 'IIR'
        % Hd = FilterDesignChebchev();
        Hd = FilterDesignButterworth2_1();
        % Hd = FilterDesignCLS0_5();

        y = filter(Hd, [x;x;x;x]);
        y= y(end-nFr+1:end,:);

    otherwise
        fprintf('Unknown filter type\n');
end

end

author: singaxiong

repo: SignalGraph


text:
function y = softmax(x)
    M = bsxfun(@minus, x, max(x, [], 1));
    numerator = exp(M);
    denominator = sum(numerator) + 0.00001;
    y = bsxfun(@rdivide, numerator, denominator) + 0.00001;
end

author: jimmy-ren


text:
function [y]=tt_exp(x, eps, varargin)
%Computation of the pointwise exponential in TT format
%   [Y]=TT_EXP(X,EPS,VARARGIN) This function computes pointwise exponential
%   using scaling and squaring method of the TT-vector X. EPS is the
%   accuracy parameter, varargins:
%       N is the number of summand for the local Taylor
%           series (N=10 by default, usually enough).
%       RMAX is the TT-rank bound.

% nrm = norm(x);
if (isa(x, 'tt_tensor'))
    nrm = tt_max_abs(x);
else
    nrm = tt_max_abs(qtttucker_to_linqtt(x, eps));
end;
n0 = floor(max(log2(nrm), 0))+1;
x = x./(2^n0);
if (isa(x, 'tt_tensor'))
    ons = tt_ones(x.n);
else
    ons = [];
    for i=1:(x.core.d)
        curons = tt_ones(x.tuck{i}.n);
        curons = qtt_tucker(curons, x.tuck{i}.d, eps);
        ons = kron(ons, curons);
    end;
end;
y = ons;

N = 10;
hdm = 'svd';
rmax = Inf;
epst = eps;
while (length(varargin)==1)
    varargin = varargin{1};
end;
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'n'
            N=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        case 'epst'
            epst=varargin{i+1};
        case 'hdm'
            hdm=varargin{i+1};

        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

for k=(N-1):-1:1
% for k=1:N-1
    y=ons+(y.*x)/k;
    y=round(y,epst,rmax);
end

for k=1:n0
    if (strcmp(hdm, 'svd'))
        y=round(y.*y,eps*(0.5^(n0-k)),rmax);
        fprintf('squaring %d\n', k);
    else
        if (isa(x, 'tt_tensor'))
%             y = mvk3(diag(y), y, eps, 'nswp', 20, 'kickrank', 2);
            y = tt_mvk4(diag(y), y, eps, 'nswp', 20);
        else
%             y = mvrk(diag(y), y, eps, 'nswp', 20, 'kickrank', 2);
            y = mvrk2(diag(y), y, eps, 'nswp', 20);
        end;
    end;
end

end

author: oseledets

repo: TT-Toolbox


text:
function [ r, seed ] = r8mat_uniform_01 ( m, n, seed )

%*****************************************************************************80
%
%% R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.
%
%    Output, real R(M,N), an array of random values between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8MAT_UNIFORM_01 - Fatal error!' );
  end

  for j = 1 : n
    for i = 1 : m

      seed = floor ( seed );

      seed = mod ( seed, 2147483647 );

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      r(i,j) = seed * 4.656612875E-10;

    end
  end

author: johannesgerer

repo: jburkardt-m


text:
classdef EigenValueAndVectorComputerFactory < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            switch cParams.type
                case 'SYMBOLIC'
                    obj = EigenValueAndVectorComputerSymbolic(cParams);
                case 'PRECOMPUTED'
                    obj = EigenValueAndVectorComputerPrecomputed(cParams);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = interp2(varargin)
% bsarray/interp2: 2-D interpolation (table lookup)
% usage: ZI = interp2(B,XI,YI);
%    or: ZI = interp2(B,XI,YI,EXTRAPVAL);
%
% arguments:
%   B - bsarray object having tensorOrder 2. The positions of the
%       x-coordinates of the underlying data array are assumed to be
%       X = s(2).*(1:Nx), where Nx is the number of data points in the X
%       dimension (i.e., second element of get(B,'dataSize')), and s
%       is the element spacing (i.e., get(B,'elementSpacing')). The
%       positions of the y-coordinates of the underlying data array are
%       assumed to be Y = s(1).*(1:Ny).
%   XI - x-coordinates of points at which to interpolate B.
%   YI - y-coordinates of points at which to interpolate B. YI must be the
%       same size as XI.
%
%   EXTRAPVAL - value to return for points in XI and YI that are outside
%       the range of X and Y, respectively. Default EXTRAPVAL = NaN.
%
%   ZI - the values of the underlying bsarray B evaluated at the points in
%       the array XI.
%

% author: Nathan D. Cahill
% email: ndcahill@gmail.com
% date: 18 April 2008

% parse input arguments
[b,xi,yi,extrapval] = parseInputs(varargin{:});

% get flag to determine if basis functions in each dimension are centred or
% shifted
m = double(get(b,'centred'));
mx = m(2); my = m(1);

% get number of data elements and coefficients, determine amount of padding
% that has been done to create coefficients
nData = get(b,'dataSize'); nDatax = nData(2); nDatay = nData(1);
n = get(b,'coeffsSize'); nx = n(2); ny = n(1);
padNum = (n-nData-1)/2; padNumx = padNum(2); padNumy = padNum(1);

% get the spacing between elements, and then construct vectors of the
% locations of the data and of the BSpline coefficients
h = get(b,'elementSpacing'); hx = h(2); hy = h(1);
xCol = hx.*((1-padNumx):(nDatax+padNumx+1))';
xDataCol = xCol((1+padNumx):(end-padNumx));
yCol = hy.*((1-padNumy):(nDatay+padNumy+1))';
yDataCol = yCol((1+padNumy):(end-padNumy));

% turn evaluation points into a column vector, but retain original size so
% output can be returned in same size as input
siz_xi = size(xi);
siz_zi = siz_xi;

% grab the BSpline coefficients
cMat = get(b,'coeffs');

% initialize some variables for use in interpolation
numelXi = numel(xi);
ziMat = zeros(numelXi,1);
p = 1:numelXi;

% Find indices of subintervals, x(k) <= u < x(k+1),
% or u < x(1) or u >= x(m-1).
kx = min(max(1+floor((xi(:)-xCol(1))/hx),1+padNumx),nx-padNumx) + 1-mx;
ky = min(max(1+floor((yi(:)-yCol(1))/hy),1+padNumy),ny-padNumy) + 1-my;
sx = (xi(:) - xCol(kx))/hx;
sy = (yi(:) - yCol(ky))/hy;

% perform interpolation
d = get(b,'degree'); dx = d(2); dy = d(1);
xflag = (~mx && mod(dx,2));
yflag = (~my && mod(dy,2));
for j=1:ceil((dx+1)/2) % loop over BSpline degree in x dimension
    Bx1 = evalBSpline(sx+j-(1+mx)/2,dx);
    Bx2 = evalBSpline(sx-j+(1-mx)/2,dx);
    for i=1:ceil((dy+1)/2) % loop over BSpline degree in y dimension
        By1 = evalBSpline(sy+i-(1+my)/2,dy);
        By2 = evalBSpline(sy-i+(1-my)/2,dy);
        for ind = 1:numelXi % loop over evaluation points, computing interpolated value
            ziMat(ind) = ziMat(ind) + cMat(ky(ind)-i+my,kx(ind)-j+mx).*By1(ind).*Bx1(ind) + ...
                cMat(ky(ind)+i-1+my,kx(ind)-j+mx).*By2(ind).*Bx1(ind) + ...
                cMat(ky(ind)-i+my,kx(ind)+j-1+mx).*By1(ind).*Bx2(ind) + ...
                cMat(ky(ind)+i-1+my,kx(ind)+j-1+mx).*By2(ind).*Bx2(ind);
        end
    end
    if yflag % add a correction factor if BSpline in y direction is shifted and of odd degree
        By1 = evalBSpline(sy+i+1/2,dy);
        for ind = 1:numelXi
            ziMat(ind) = ziMat(ind) + By1(ind).*...
                (cMat(ky(ind)-i-1,kx(ind)-j+mx).*Bx1(ind) + cMat(ky(ind)-i-1,kx(ind)+j-1+mx).*Bx2(ind));
        end
    end
end
if xflag % add a correction factor if BSpline in x direction is shifted and of odd degree
    Bx1 = evalBSpline(sx+j+1/2,dx);
    for i=1:ceil((dy+1)/2)
        By1 = evalBSpline(sy+i-(1+my)/2,dy);
        By2 = evalBSpline(sy-i+(1-my)/2,dy);
        for ind = 1:numelXi
            ziMat(ind) = ziMat(ind) + Bx1(ind).*...
                (cMat(ky(ind)-i+my,kx(ind)-j-1).*By1(ind) + cMat(ky(ind)+i-1+my,kx(ind)-j-1).*By2(ind));
        end
    end
end

% perform extrapolation
outOfBounds = xi(:)<xDataCol(1) | xi(:)>xDataCol(nDatax) | yi(:)<yDataCol(1) | yi(:)>yDataCol(nDatay);
ziMat(p(outOfBounds)) = extrapval;

% reshape result to have same size as input xi
zi = reshape(ziMat,siz_zi);
varargout{1} = zi;

%% subfunction parseInputs
function [b,xi,yi,extrapval] = parseInputs(varargin)

nargs = length(varargin);
error(nargchk(3,4,nargs));

% Process B
b = varargin{1};
if ~isequal(b.tensorOrder,2)
    error([mfilename,'parseInputs:WrongOrder'], ...
        'bsarray/interp2 can only be used with bsarray objects having tensor order 2.');
end

% Process XI
xi = varargin{2};
if ~isreal(xi)
    error([mfilename,'parseInputs:ComplexInterpPts'], ...
        'The interpolation points XI should be real.')
end

% Process YI
yi = varargin{3};
if ~isreal(yi)
    error([mfilename,'parseInputs:ComplexInterpPts'], ...
        'The interpolation points YI should be real.')
end
if ~isequal(size(xi),size(yi))
    error([mfilename,'parseInputs:YIXINotSameSize'], ...
        'YI must be the same size as XI');
end

% Process EXTRAPVAL
if nargs > 3
    extrapval = varargin{4};
else
    extrapval = [];
end
if isempty(extrapval)
    extrapval = NaN;
end
if ~isscalar(extrapval)
    error([mfilename,':NonScalarExtrapValue'],...
        'EXTRAP option must be a scalar.')
end

author: Sable

repo: mcbench-benchmarks


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

function [Mdh]=DHmatrix(theta,d,a,alpha)
%%output: DH matrix
%%input: DH parameters:
%%theta: rotation along x(n) axis
%%d: traslation along z(n) axis
%%a: translation along x(n+1) axis
%%alpha: rotation along x(n+1) axis
%%all angles expressed in degrees

Mdh=[cosd(theta) -sind(theta)*cosd(alpha) sind(theta)*sind(alpha) a*cosd(theta);
     sind(theta) cosd(theta)*cosd(alpha)  -cosd(theta)*sind(alpha) a*sind(theta);
     0,sind(alpha),cosd(alpha),d;
     0,0,0,1];

end


text:
function [X,Y,Z]=oblate0(eta,theta,psi)
% [X,Y,Z]=oblate0(eta,theta,psi) defines
% oblate spheroidal coordinates
X=cosh(eta).*sin(theta).*cos(psi);
Y=cosh(eta).*sin(theta).*sin(psi);
Z=sinh(eta).*cos(theta);


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IR Tools
% April 2018
%
% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as
% part of the package.
%
% Copyright 2018 Silvia Gazzola, University of Bath, Per Christian Hansen,
% Technical University of Denmark and James G. Nagy, Emory University.
%
% Contact: s.gazzola@bath.ac.uk, pcha@dtu.dk, jnagy@emory.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Installating
%   IRtools_setup       - Set up search paths to IRtools
%
% Example scripts:
%   EXblur_cgls_hybrid  - Example script, speckle deblurring problem
%   EXdiffusion_rrgmres - Example script, inverse diffusion problem
%   EXinvinterp2_cgls   - Example script, inverse interpolation problem
%   EXnmr_cgls_mrnsd    - Example script, 2D NMR relaxometry
%   EXsparsity          - Example script, deblurring without and with sparsity
%
% Iterative methods:
%   IRart              - Algebraic Reconstruction Technique
%   IRcgls             - Conjugate Gradient algorithm for Least Squares problems
%   IRconstr_ls        - Least squares solver with box and energy constraints
%   IRell1             - Least squares solver with 1-norm penalization term
%   IRenrich           - Enriched CGLS method
%   IRfista            - FISTA algorithm for constrained least squares problems
%   IRget              - Get options for IR Tools functions
%   IRhtv              - Least squares solver with heuristic total variation penalization
%   IRhybrid_fgmres    - Hybrid version of FGMRES for enforcing 1-norm penalization
%   IRhybrid_gmres     - Hybrid version of GMRES algorithm for square systems
%   IRhybrid_lsqr      - Hybrid version of LSQR algorithm
%   IRirn              - Least squares solver with 1-norm penalization term
%   IRmrnsd            - Modified Residual Norm Steepest Descent method
%   IRnnfcgls          - Modified flexible CGLS for nonnegatively constrained LS problems
%   IRrestart          - Restarted Krylov subspace methods
%   IRrrgmres          - Range Restricted GMRES for square systems
%   IRset              - Set options for IR Tools functions
%   IRsirt             - Simuletaneous Iterative Reconstruction Technique
%
% Operators (functions) for some test problems:
%   OPdiffusion        - The forward computation and its adjoint for PRdiffusion
%   OPinvinterp2       - The forward computation and its adjoint for PRinvinterp2
%   OPnmr              - The forward computation and its adjoint for PRnmr
%
% Test problems and associated functions:
%   PRblur             - Generates data for use in image deblurring problems
%   PRblurdefocus      - Image deblurring problem with a defocus point spread function
%   PRblurgauss        - Image deblurring problem with a Gaussian point spread function
%   PRblurmotion       - Image deblurring problem with linear motion blur
%   PRblurrotation     - Image deblurring problem with rotational motion blur
%   PRblurshake        - Image deblurring problem with random camera motion (shaking)
%   PRblurspeckle      - Image deblurring problem with a speckle point spread function
%   PRdiffusion        - Generates data for use in an inverse diffusion problems
%   PRget              - Get options for IR Tools test problems
%   PRinvinterp2       - Generates data for use in 2D inverse interpolation problems
%   PRnmr              - Generates data for 2D NMR relaxometry problems
%   PRnoise            - Add noise of a specified distribution and level
%   PRseismic          - Generates data for seismic travel-time tomography problems
%   PRset              - Set options for IR Tools test problems
%   PRshowb            - Show the right-hand side b (the data) in IR Tools
%   PRshowx            - Show the solution x in IR Tools
%   PRspherical        - Generates data for spherical means tomography problems
%   PRtomo             - Generates data for X-ray tomographic reconstruction problems

repo: IRtools


text:
function [atomicWeights] = parse_Atomic_Weights_and_Isotopic_Compositions_for_All_Elements
% Parses NIST data on atomic weights
%
% USAGE:
%
%    [atomicWeights] = parse_Atomic_Weights_and_Isotopic_Compositions_for_All_Elements
%
% OUTPUT:
%    atomicWeights:    atomic weights of the isotopes
%
% The atomic weight J. S. Coursey, D. J. Schwab, and R. A. Dragoset
% NIST, Physics Laboratory, Office of Electronic Commerce in Scientific and
% Engineering Data.
% The atomic weights are available for elements 1 through 112, 114, & 116 and
% isotopic compositions or abundances are given when appropriate. The atomic
% weights data were published by T.B. Coplen in `Atomic Weights of the Elements
% 1999`, (and include changes reported from the 2001 review in `Chem. Int., 23, 179 (2001)`)
% and the isotopic compositions data were published by K.J.R. Rosman2 and P.D.P. Taylor3
% in `Isotopic Compositions of the Elements 1997`. The relative atomic masses of the
% isotopes data were published by G. Audi4 and A. H. Wapstra5 in `The 1995 Update To
% The Atomic Mass Evaluation`.
% http://physics.nist.gov/PhysRefData/Compositions/
%
% .. Author: - Ronan Fleming, 9 March 09

fid=fopen('Atomic_Weights_and_Isotopic_Compositions_for_All_Elements.txt','r');

nElements=2816/8;
p=1;
ele=1;
fieldNames={'AtomicNumber';
'AtomicSymbol';
'MassNumber';
'RelativeAtomicMass';
'IsotopicComposition';
'StandardAtomicWeight';
'Notes'};

for x=1:2816
    line= fgetl(fid);
    if ~isempty(line)
        if p==2 || p==7
            tmp=textscan(line,'%s%s\n',2815,'Delimiter','=(','TreatAsEmpty',' ');
            if ~isempty(tmp{2})
                atomicWeights.data(ele).(fieldNames{p})=tmp{2};
            end
            if p==2
                tmp2=tmp{2};
                if ~isempty(tmp2)
                    atomicWeights.AtomicSymbol{ele}=tmp2{1};
                end
            end
        else
            tmp=textscan(line,'%s%f\n',2815,'Delimiter','=(','TreatAsEmpty',' ');
            atomicWeights.data(ele).(fieldNames{p})=tmp{2};
        end
        p=p+1;
    else
        ele=ele+1;
        p=1;
    end
end

fclose(fid);

author: opencobra

repo: cobratoolbox


text:
function c = demean(c,varargin);

% DEMEAN removes the mean of each trace.
%
% C = DEMEAN(C) removes the DC offset from each trace. In most cases it is
% unnecessary to call this function directly. By default, all traces are
% demeaned and detrended when they are loaded into a correlation object.
% This is one of the assumptions of the correlation toolbox.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

c.W = demean(c.W);

author: geoscience-community-codes

repo: GISMO


text:
close all
clear all
rootDir = '\\Snarp\u1\data\reading_longitude\fmri';

cd(rootDir);
[fileNum dirName] = countDirs(pwd);

for count = 3:fileNum
    try
        currentDir = dirName{count};

        cd(fullfile(rootDir,currentDir));
        mrVista;
        vw = getSelectedInplane;
        motionCompPlotMI(vw);

    end

    close all;
    clear all;

    rootDir = '\\Snarp\u1\data\reading_longitude\fmri';

	cd(rootDir);
	[fileNum dirName] = countDirs(pwd);

    pause(0.001);
end


text:
function    g = lfmGradientUpsilon(gamma,sigma2,gradThetaGamma,t1, t2, mode)

% LFMGRADIENTUPSILON Gradient of the function \upsilon(z) with respect to
% one of the hyperparameters of the kernel: m_k, C_k, D_k, m_r, C_r or D_r.
% FORMAT
% DESC Computes the gradient of the function \upsilon(z) with respect to
% one of the parameters of the system (mass, spring or damper).
% ARG gamma : Gamma value of the system.
% ARG sigma2 : length scale of latent process.
% ARG gradThetaGamma : Vector with the gradient of gamma1 and gamma2 with
% respect to the desired parameter.
% ARG Tt1 : first time input (number of time points 1 x number of time points 2).
% ARG Tt2 : second time input (number of time points 1 x number of time points 2).
% ARG mode : indicates in which way the variables t1 and t2 must be
% transposed
% RETURN g : Gradient of the kernel with respect to the desired
% parameter.
%
% COPYRIGHT : David Luengo, Mauricio Alvarez 2008
%
% SEEALSO : lfmKernGradient, lfmXlfmKernGradient, lfmGradientH

% KERN

%Parameters of the function

sigma = sqrt(sigma2);

% Initialization of vectors and matrices

switch mode
    case 1
        % In this mode the Tt1 is actually the Tt1 and Tt2 is actually Tt2
        % Z1 is a matrix with order given by Tt1 and z2 is a vector with
        % order given by Tt2
        Tt11 = repmat(t1, 1, size(t2, 1));
        Tt22 = repmat(t2', size(t1, 1), 1);
        Z1 = (Tt11-Tt22)/sigma - sigma*gamma/2;
        z2 = t2/sigma + sigma*gamma/2;
        Z2 = repmat(z2.', size(t1, 1), 1);
        [preFactorZ1, gradPreFactorZ1] = quadrantWofzGrad(Z1, sigma,1);
        [preFactorz2, gradPreFactorz2] = quadrantWofzGrad(z2, sigma,1);
        Tt1minusTt2 = Tt11 - Tt22;
        preFactor1 = exp(sigma2*(gamma^2)/4 - gamma*(Tt1minusTt2)...
            + log(sigma2*gamma - 2*(Tt1minusTt2)) + log(gradThetaGamma));
        preFactor2 = exp(-(Tt1minusTt2).*(Tt1minusTt2)/sigma2 ...
            + log(gradThetaGamma) + log(gradPreFactorZ1));
        preFactor3 = exp(repmat((-t2.*t2/sigma2).',size(t1, 1), 1) ...
            - gamma*Tt11 + log(gradThetaGamma)...
            + log(Tt11.*repmat(preFactorz2.',size(t1, 1), 1) ...
            + repmat((gradPreFactorz2).',size(t1, 1), 1)));
    case 2

        % In this mode the Tt1 is actually the Tt2 and Tt2 is actually Tt1
        % Z1 is a matrix with the order given by Tt2 and z2 is a vector
        % with the order given by Tt1

        Tt11 = repmat(t2', size(t1, 1), 1);
        Tt22 = repmat(t1, 1, size(t2, 1));
        Z1 = (Tt11-Tt22)/sigma - sigma*gamma/2;
        z2 = t1/sigma + sigma*gamma/2;
        Z2 = repmat(z2,  1, size(t2, 1));
        [preFactorZ1, gradPreFactorZ1] = quadrantWofzGrad(Z1, sigma,1);
        [preFactorz2, gradPreFactorz2] = quadrantWofzGrad(z2, sigma,1);
        Tt1minusTt2 = Tt11 - Tt22;
        preFactor1 = exp(sigma2*(gamma^2)/4 - gamma*(Tt1minusTt2)...
            + log(sigma2*gamma - 2*(Tt1minusTt2)) + log(gradThetaGamma));
        preFactor2 = exp(-(Tt1minusTt2).*(Tt1minusTt2)/sigma2 ...
            + log(gradThetaGamma) + log(gradPreFactorZ1));
        preFactor3 = exp( repmat(-t1.*t1/sigma2,1, size(t2, 1))...
            - gamma*Tt11 + log(gradThetaGamma) ...
            + log(Tt11.*repmat(preFactorz2,1, size(t2, 1)) + repmat(gradPreFactorz2,1, size(t2, 1))));

    case {3,4}

        % In this mode Tt1 is actually Tt1 and Tt2 is zero.
        % Z1 is a vector and Z2 is a scalar
        z1 = t1/sigma - sigma*gamma/2;
        z2 = sigma*gamma/2;
        Z1 = repmat(z1, 1, size(t2,1));
        Z2 = repmat(z2, size(t1, 1), size(t2,1));
        [preFactorz1, gradPreFactorz1] = quadrantWofzGrad(z1, sigma,1);
        [preFactorz2, gradPreFactorz2] = quadrantWofzGrad(z2, sigma,1);
        preFactor1 = repmat(exp(sigma2*(gamma^2)/4 - gamma*(t1) ...
            + log(sigma2*gamma - 2*t1) + log(gradThetaGamma)), 1, size(t2, 1));
        preFactor2 = repmat(exp(-t1.*t1/sigma2 + log(gradThetaGamma)...
            + log(gradPreFactorz1)), 1, size(t2, 1));
        preFactor3 = repmat(exp(- gamma*t1 + log(gradThetaGamma) ...
            + log(t1*preFactorz2 + gradPreFactorz2)),1, size(t2, 1));

end

g = zeros(size(preFactor1));
% Evaluation of Upsilon when real(Z1)>=0 and real(Z2)>=0
ind = (real(Z1)>=0) & (real(Z2)>=0);
if any(any(ind))
    g(ind) = preFactor1(ind) - preFactor2(ind) + preFactor3(ind);
end;
% Evaluation of Upsilon when real(Z1)<0 and real(Z2)>=0
ind = (real(Z1)<0) & (real(Z2)>=0);
if any(any(ind))
    g(ind) = preFactor2(ind) + preFactor3(ind);
end
% Evaluation of Upsilon when real(Z1)>=0 and real(Z2)<0
ind = (real(Z1)>=0) & (real(Z2)<0);
if any(any(ind))
    g(ind) = - preFactor2(ind) - preFactor3(ind);
end;
% Evaluation of Upsilon when real(Z1)<0 and real(Z2)<0
ind = (real(Z1)<0) & (real(Z2)<0);
if any(any(ind))
    g(ind) = -preFactor1(ind) + preFactor2(ind) - preFactor3(ind);
end;

if mode==4
    g = g.';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Y,SetupStruc] = OnMVDR(S,Transfer,SetupStruc)
K = SetupStruc.K;
hop = SetupStruc.hop;
win = SetupStruc.win;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(S,2);
frame_N = size(S,3);
blockWin = repmat(win,[1 N frame_N]);
X = fft(S.*blockWin);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(K,frame_N,Num);
%%%%%%%%%%%%%%%%%%%%%%%%%% Initial the coefficients stored in 'SetupStruc'
if(~isfield(SetupStruc,'W'))
    SetupStruc.W = permute(Transfer,[2 3 1])/N;
    SetupStruc.V = zeros(N,Num,K_m);
    sign_Int = 1;
    SetupStruc.R = zeros(N,N,K_m);
    SetupStruc.WNG = [];
%     SetupStruc.miu = [];
%     SetupStruc.mag = [];
else
    sign_Int = 0;
end
WNG_set = -5;
WNG = zeros(1,K_m);
% miu_ = zeros(1,K_m);
% mag = zeros(1,K_m);
delta2 = 10^(WNG_set/10);
for i = 2:K_m
    X_f = permute(X(i,:,:),[2 3 1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% MVDR processing
    W_f = zeros(N,Num);
    Steer = permute(Transfer(i,:,:),[2 3 1]);
    if(sign_Int==1 || SetupStruc.alpha==1)
        R = X_f*X_f'/frame_N;
    else
        R = (1-SetupStruc.alpha)*SetupStruc.R(:,:,i)+SetupStruc.alpha*(X_f*X_f')/frame_N;
        SetupStruc.R(:,:,i) = R;
    end
    miu = 2/3/sum(real(diag(R)));
    for j = 1:Num
        h = Steer(:,j);
        mo2 = h'*h;
        Wc = h/mo2;
        Pc0 = eye(N)-Wc*h';
        V = SetupStruc.V(:,j,i);
        W = SetupStruc.W(:,j,i);
        V = Pc0*(V-miu*R*W);
        b2 = 1/delta2-1/mo2;
        if(V'*V<=b2)
            W = Wc+V;
        else
            W = Wc+sqrt(b2)*V/sqrt(V'*V);
        end
        W_f(:,j) = W;
        SetupStruc.V(:,j,i) = V;
    end
    SetupStruc.W(:,:,i) = W_f;
    WNG(i) = -10*log10(W_f(:,1)'*W_f(:,1));
%     mag(i) = abs(W_f(:,1)'*Steer(:,1));
%     miu_(i) = miu;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Y_ = W_f'*X_f;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
SetupStruc.WNG = [SetupStruc.WNG;WNG];
% SetupStruc.miu = [SetupStruc.miu;miu_];
% SetupStruc.mag = [SetupStruc.mag;mag];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end
return;

repo: Voice-Separation-and-Enhancement


text:
function ret = spm_ov_contour(varargin)
% Contour tool - plugin for spm_orthviews
%
% This routine is a plugin to spm_orthviews. For general help about
% spm_orthviews and plugins type
%             help spm_orthviews
% at the MATLAB prompt.
%__________________________________________________________________________
% Copyright (C) 2014-2016 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_ov_contour.m 6774 2016-04-20 12:25:08Z guillaume $

switch lower(varargin{1})
    % Context menu and callbacks
    case 'context_menu'
        item0 = uimenu(varargin{3}, ...
            'Label', 'Contour',...
            'UserData',struct('nblines',3,'style','r-','id',varargin{2}));
        item1 = uimenu(item0,...
            'Label','Display onto');
        item1_1 = uimenu(item1,...
            'Label','all',...
            'Callback',@(hObj,event) contour_menu(hObj,event,varargin{2},Inf));
        item1_2 = uimenu(item1,...
            'Label','all but current',...
            'Callback',@(hObj,event) contour_menu(hObj,event,varargin{2},NaN));
        item1_3 = uimenu(item1,...
            'Label','select...',...
            'Callback',@(hObj,event) contour_menu(hObj,event,varargin{2}));
        item2 = uimenu(item0,...
            'Label','Options');
        item2_1 = uimenu(item2,...
            'Label', 'Number of lines...',...
            'UserData',struct('str','Number of lines', 'field','nblines'),...
            'Callback',@(hObj,event) contour_options(hObj,event,item0));
        item2_2 = uimenu(item2,...
            'Label','Line style...',...
            'UserData',struct('str','Line style', 'field','style'),...
            'Callback',@(hObj,event) contour_options(hObj,event,item0));
        %'Callback',@(hObj,event) contour_display(hObj,event,varargin{2}));
        ret = item0;
    case 'display'
        if nargin == 1
            varargin{2} = spm_input('Select image', '!+1', 'e');
        end
        contour_display(varargin{2:end});
    case 'redraw'
        contour_redraw(varargin{2:end});
    otherwise
end

%==========================================================================
function contour_options(hObj,event,hM)
Finter = spm_figure('FindWin', 'Interactive');
spm_input('!DeleteInputObj',Finter);
UDm = get(hObj,'UserData');
UDc = get(hM,'UserData');
if isnumeric(UDc.(UDm.field))
    N = 1;
    T = 'e';
else
    N = Inf;
    T = 's';
end
UDc.(UDm.field) = spm_input(UDm.str, '!+1', T, UDc.(UDm.field), N);
set(hM,'UserData',UDc);
try, contour_redraw(UDc.id); end

%==========================================================================
function contour_menu(hObj,event,i,varargin)

global st

allOpts = get(get(hObj,'parent'),'children');
current = strcmp(get(hObj,'Checked'),'on');
if any(strcmp(get(allOpts,'Checked'),'on'))
    set(allOpts, 'Checked','off');
    contour_delete(i);
    st.vols{i} = rmfield(st.vols{i},'contour');
end
if ~current
    set(hObj, 'Checked','on');
    contour_display(i,varargin{:});
end

%==========================================================================
function contour_display(i,o)

global st

if nargin < 2
    o = spm_input('Select image(s)', '!+1', 'e', ...
        num2str(spm_orthviews('valid_handles')));
    o = intersect(spm_orthviews('valid_handles'),o);
elseif isinf(o)
    o = spm_orthviews('valid_handles');
elseif isnan(o)
    o = setxor(spm_orthviews('valid_handles'),i);
end

try
    hM = findobj(st.vols{i}.ax{1}.cm,'Label','Contour');
    UD = get(hM,'UserData');
    nblines = UD.nblines;
    linestyle = UD.style;
catch
    nblines = 3;
    linestyle = 'r-';
end
linewidth = 1;

contour_delete(i);

lh = {};
sw = warning('off','MATLAB:contour:ConstantData');
for d = 1:3
    CData = sqrt(sum(get(st.vols{i}.ax{d}.d,'CData').^2, 3));
    CData(isinf(CData)) = NaN;
    CData(isnan(CData)) = 0;
    for h = o(:)'
        set(st.vols{h}.ax{d}.ax,'NextPlot','add');
        [C,lh{end+1}] = ...
            contour(st.vols{h}.ax{d}.ax,CData,...
            nblines,linestyle,'LineWidth',linewidth);
    end
end
warning(sw);
set(cat(1,lh{:}),'HitTest','off');

st.vols{i}.contour.images = o;
st.vols{i}.contour.handles = lh;

%==========================================================================
function contour_redraw(i,varargin) %i, TM0, TD, CM0, CD, SM0, SD

global st

contour_delete(i);
contour_display(i,st.vols{i}.contour.images);

%==========================================================================
function contour_delete(i)

global st

try, delete(cat(1,st.vols{i}.contour.handles{:})); end

author: spm


text:
function L = computeL(nu)
    L_tmp = cell(nu,1);

    for i=1:nu
       L_tmp{i,1} = eye(2);
    end

    L = cell2mat(L_tmp);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
classdef TestFastFeatureDetector
    %TestFastFeatureDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            img = imread(TestFastFeatureDetector.im);
            obj = cv.FastFeatureDetector('Threshold',10, 'NonmaxSuppression',true);
            assert(obj.Threshold == 10);
            assert(obj.NonmaxSuppression == true);
            typename = obj.typeid();

            kpts = obj.detect(img);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestFastFeatureDetector.kfields, fieldnames(kpts))));
        end

        function test_detect_imgset
            img = imread(TestFastFeatureDetector.im);
            imgs = {img, img};
            obj = cv.FastFeatureDetector();

            kpts = obj.detect(imgs);
            validateattributes(kpts, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, {'vector'}), kpts);
            cellfun(@(kpt) assert(all(ismember(TestFastFeatureDetector.kfields, fieldnames(kpt)))), kpts);
        end

        function test_detect_mask
            img = imread(TestFastFeatureDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.FastFeatureDetector();
            kpts = obj.detect(img, 'Mask',mask);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestFastFeatureDetector.kfields, fieldnames(kpts))));
            xy = cat(1, kpts.pt);
            assert(all(xy(:,1) <= ceil(size(img,2)/2)));
        end

        function test_detect_types
            types = {'TYPE_5_8', 'TYPE_7_12', 'TYPE_9_16'};
            img = imread(TestFastFeatureDetector.im);
            obj = cv.FastFeatureDetector();
            for i=1:numel(types)
                obj.Type = types{i};
                kpts = obj.detect(img);
                validateattributes(kpts, {'struct'}, {'vector'});
            end
        end

        function test_error_1
            try
                cv.FastFeatureDetector('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function name = p17_name ( )

%*****************************************************************************80
%
%% P17_NAME returns the name of problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'L2(X-Z)^2';

  return
end


text:
function qual_obj_array = load_qual_files(filenames)
% This function loads quality measures from files into an array of
% Quality objects.
% Inputs:
%   filenames: A list of files containing quality measures.
% Output:
%   qual_obj_array: An array of Quality objects containing the same
%     number of elements as filenames.

assert(iscell(filenames))

num_systems = length(filenames);
qual_obj_array = Quality.empty(num_systems,0);

log_info('quality:\n');

for ii=1:num_systems
    qual_obj_array(ii) = Quality.read(filenames{ii});
    log_info(' %s\n',filenames{ii});
end

author: nesl


text:
function [ resized_patch ] = get_pixels(im, pos, sz, resize_target)

if isscalar(sz)  %square sub-window
    sz = [sz, sz];
end

%make sure the size is not to small
if sz(1) < 1
    sz(1) = 2;
end
if sz(2) < 1
    sz(2) = 2;
end

xs = floor(pos(2)) + (1:sz(2)) - floor(sz(2)/2);
ys = floor(pos(1)) + (1:sz(1)) - floor(sz(1)/2);

%check for out-of-bounds coordinates, and set them to the values at
%the borders
xs(xs < 1) = 1;
ys(ys < 1) = 1;
xs(xs > size(im,2)) = size(im,2);
ys(ys > size(im,1)) = size(im,1);

%extract image
im_patch = im(ys, xs, :);

if isempty(resize_target)
    resized_patch = im_patch;
else
    resized_patch = mexResize(im_patch,resize_target,'auto');
end
end

repo: AutoTrack


text:
function [varargout] = rdivide(varargin)

% RDIVIDE compile the missing mex file on the fly

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);
  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear all;
close all;
clc;
export = true;
configs = cell(1, 5);
configs{1} = create_config([1 1 1 1 1 1], 4, '1-4');
configs{2} = create_config([2 1 1 1 1 1], 4, '2.1-4');
configs{3} = create_config([4 2 1 1 1 1], 4, '42.1-4');
configs{4} = create_config([2 2 2 2 2 2], 4, '2-4');
configs{5} = create_config([2 2 2 2 2 2], 8, '2-8');

nc = numel(configs);

points = [150  275  490  885 1600 2885 5205 9400 16980];
nns = numel(points);
clustering_acc_perc_mat = zeros(nc, nns);
spr_perc_mat = zeros(nc, nns);
spr_error_perc_mat = zeros(nc, nns);

for c=1:nc
    config = configs{c};
    if isempty(config)
        continue;
    end
    solver_name = sprintf('ssc_mc_omp_%s', config.name);
    filepath = sprintf('bin/solver_%s.mat', solver_name);
    load(filepath);

    clustering_acc_perc_mat(c, :)  = result.clustering_acc_perc_arr;
    spr_perc_mat(c, :) = result.spr_perc_arr;
    spr_error_perc_mat(c, :) = 100 * result.spr_error_arr;
    fprintf('%s Points: \n', solver_name);
    fprintf('%4d ', result.num_points_arr);
    fprintf('\n');

    fprintf('Accuracy: \n');
    fprintf('%.2f ', result.clustering_acc_perc_arr);
    fprintf('\n');

    % plot_bench_results(result, solver_name);
end

legends = {'OMP', 'MC-OMP 2.1-4', 'MC-OMP 42.1-4', 'MC-OMP 2-4', 'MC-OMP 2-8'};

fprintf('\n\nAccuracy: \n');
disp(clustering_acc_perc_mat);

fprintf('\nSubspace preserving representations: \n');
disp(spr_perc_mat);

fprintf('\nSubspace representation error: \n');
disp(spr_error_perc_mat);

styles = spx.graphics.plot_styles();
mf = spx.graphics.Figures;
mf.new_figure('Clustering Accuracy');
for c=1:nc
    style = styles{c};
    plot(points, clustering_acc_perc_mat(c, :), style);
    hold on;
end
xlabel('Points');
ylabel('Clustering accuracy (%)');
ylim([40 100]);
grid on;
legend(legends, 'Location', 'SouthEast');
if export
set(gcf, 'units', 'inches', 'position', [.8 .8 4 3]);
set(gca,'FontSize',8);
set(findall(gcf,'type','text'),'FontSize',8);
export_fig plots/clustering_accuracy.png -r120 -nocrop;
export_fig plots/clustering_accuracy.pdf;
end

mf.new_figure('Subspace preserving representations');
for c=1:nc
    style = styles{c};
    plot(points, spr_perc_mat(c, :), style);
    hold on;
end
xlabel('Points');
ylabel('Subspace preserving representation percentage (%)');
ylim([0 100]);
grid on;
legend(legends, 'Location', 'SouthEast');
if export
set(gcf, 'units', 'inches', 'position', [.8 .8 4 3]);
set(gca,'FontSize',8);
set(findall(gcf,'type','text'),'FontSize',8);
export_fig plots/subspace_preserving_representation_perc.png -r120 -nocrop;
export_fig plots/subspace_preserving_representation_perc.pdf;
end

mf.new_figure('Subspace representation error');
for c=1:nc
    style = styles{c};
    plot(points, spr_error_perc_mat(c, :), style);
    hold on;
end
xlabel('Points');
ylabel('Subspace representation error  (%)');
ylim([0 50]);
grid on;
legend(legends);
if export
set(gcf, 'units', 'inches', 'position', [.8 .8 4 3]);
set(gca,'FontSize',8);
set(findall(gcf,'type','text'),'FontSize',8);
export_fig plots/subspace_representation_error_perc.png -r120 -nocrop;
export_fig plots/subspace_representation_error_perc.pdf;
end

author: indigits

repo: sparse-plex


text:
%%**********************************************************************
%% HSDNTpred: Compute (dX,dy,dZ) for NT direction.
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

 function [par,dX,dy,dZ,coeff,L,hRd] = ...
           HSDNTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);

    global schurfun schurfun_par
%%
%% compute NT scaling matrix
%%
    [par.W,par.G,par.sv,par.gamx,par.gamz,par.dd,par.ee,par.ff] = ...
     NTscaling(blk,X,Z,Zchol,invZchol);
%%
%% compute schur matrix
%%
    m = par.m;
    schur = sparse(m+2,m+2);
    UU = []; EE = [];
    dX = cell(size(blk,1),1); dy = []; dZ = cell(size(blk,1),1);
%%
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
       elseif strcmp(pblk{1},'q');
          [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.ee);
       elseif strcmp(pblk{1},'s')
          if isempty(schurfun{p})
             schur = schurmat_sblk(blk,At,par,schur,p,par.W);
          elseif isstr(schurfun{p})
             schurtmp = sparse(m,m);
             if ~isempty(par.permZ{p})
                Wp = par.W{p}(par.permZ{p},par.permZ{p});
             else
                Wp = par.W{p};
             end
             eval(['schurtmp = ',schurfun{p},'(Wp,Wp,schurfun_par(p,:));']);
             schur = schur + schurtmp;
          end
       end
    end
%%
%% compute rhs
%%
    [rhs,EinvRc,hRd] = HSDNTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
    par.addschur   = par.kap/par.tau;
    schur(m+1,m+1) = schur(m+1,m+1) + par.kap/par.tau;
    schur(m+2,m+2) = schur(m+2,m+2) + par.addschur;
    [xx,coeff,L] = HSDlinsysolve(par,schur,UU,EE,par.Umat,rhs);
%%
%% compute (dX,dZ)
%%
    [par,dX,dy,dZ] = HSDNTdirfun(blk,At,par,Rd,EinvRc,xx);
%%**********************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function b = r8to_to_r8ge ( n, a )

%*****************************************************************************80
%
%% R8TO_TO_R8GE copies a R8TO matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8TO storage format is used for a Toeplitz matrix, which is constant
%    along diagonals.  Thus, in an N by N Toeplitz matrix, there are at most
%    2*N-1 distinct entries.  The format stores the N elements of the first
%    row, followed by the N-1 elements of the first column (skipping the
%    entry in the first row).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(2*N-1), the R8TO matrix.
%
%    Output, real B(N,N), the R8GE matrix.
%
  for i = 1 : n
    b(i,1:i-1) = a(n+i-1:-1:n+1,1);
    b(i,i:n) = a(1:n-i+1,1);
  end

  return
end


text:
function trialmodels = single_trial_setup_models(EXPT, onsets, eventdesign, s, varargin)

    % trialmodels = single_trial_setup_models(EXPT, onsets, eventdesign, s, varargin)
    %Constructs a trial-level design matrix, estimates betas and writes them to
    %img files.
    %
    %USAGE: single_trial_analysis(EXPT,onsets,eventdesign,submat);
    %       e.g.,
    %       [onsets2 runlabels eventlabels rundesign eventdesign] = ...
    %       onsets2singletrial(subjects,EXPT.SPM.nscan,'TR',2,eventnames,'noadd');
    %
    %       single_trial_analysis(EXPT,onsets2,eventdesign,1,'gamma3',60);
    %
    %       single_trial_analysis(EXPT,onsets2,eventdesign,2,'gamma3',60, ...
    %       'ratings',ratings,'conditioncontrasts',[-3 -1 1 3], ...
    %       'conditionnames',{'warm' 'low' 'medium' 'hot'}, 'mvmt', ravols_movement_params);
    %
    %INPUTS:
    %   EXPT - EXPT object
    %       Required fields:
    %       subjects
    %       TR
    %       mask
    %       SPM.nscan
    %       FILES.im_files
    %
    %   onsets - 'cell' format onsets (see below)
    %       to convert from structure format, see onsets2singletrial
    %
    %   eventdesign - stick (delta) function indicating condition onsets;
    %                 returned as output of onsets2singletrial
    %
    %   submat - subject (or vector of subject) indices
    %   basistype (optional) - basis type (default: 'gamma3')
    %   hipass (optional) - high-pass filter (default: 60)
    %   startslice (optional) - starting slice
    %
    % Optional inputs:
    % conditioncontrasts: followed by contrasts across expt'l conditions
    % which are coded in eventdesign.  e.g., [-3 -1 1 3]';
    %
    % conditionnames: names for each condition in eventdesign: {'warm' 'low' 'medium' 'hot'};
    %
    % ratings : followed by behavioral data values
    % in same format as onsets
    % See onsets2singletrial for ways to do this.
    %
    %OUTPUTS:
    %   -saves a mat file with the suffix 'design' that contains the design mtx
    %
    %see trial_level_beta3.m for more information
    %
    %A NOTE ABOUT ONSETS:
    %The SCAN lab uses two different formats for its onsets. The
    %'structure' format contains a different field for each subject, and a
    %different sub-field for each run, with sub-fields for each condition that
    %contain the onsets. The 'cell' format is a cell vector, where each cell
    %contains 1 subject; within each subject-cell is another cell vector, where
    % each cell contains onsets for a particular run/condition. If, for example, you have 2
    %conditions and 2 runs, the cell array will be organized as follows:
    %The 1st cell will contain the onsets of all events from condition 1 in
    %run 1, the 2nd cell will contain the onsets of events from condition 2 in
    %run 1, the 3rd cell will contain onsets of events from condition 1 in run
    %2, etc. Specifically for single-trial analysis, the onsets must be in cell
    %format (and in TRs), with the added requirement that the onsets be adjusted
    %so that the first onset of each run is equal to the onset of the last scan
    %(*not event*) of the previous run. For example, if you have 100 scans in
    %each run, the first onset of the 1st run should be 0, the first onset of
    %the 2nd run should be 100, etc.
    %A useful function in this regard is 'onsets2singletrial',which converts
    %structure-format onsets to cell-format and adjusts the timing. The output of
    %onsets2singletrial can be fed directly into single_trial_analysis.
    %
    % Sam Gershman, 2006
    % Modified by Tor Wager, March 2007

    %-------- SOME PARAMETERS ---------------%
    TR = EXPT.TR;                   %TR
    nscans = EXPT.SPM.nscan;        %vector containing number of scans in each run
    nruns = length(nscans);         %number of runs
    %----------------------------------------%

    %----- SET DEFAULTS ---------------------%
    basistype = 'gamma3';
    hipass = 100;
    domvmt2 = 0;
    %----------------------------------------%

    disp(' '); fprintf('Building design matrix for %3.0f...\n',s);
    disp(['Name of this subject in EXPT.subjects : ' EXPT.subjects{s}]);

    ratings = []; conditioncontrasts = []; conditionnames = []; startslice = 1;

    for i = 1:length(varargin)
        if ischar(varargin{i})

            switch varargin{i}

                case 'ratings', ratings =  varargin{i+1};
                    disp('Found ratings. Will include in model.');

                case {'conditionnames'}, conditionnames = varargin{i+1};
                case {'conditioncontrasts'}, conditioncontrasts = varargin{i+1};

                case {'mvmt','movement'}, mvmt =  varargin{i+1}; disp('Found movement parameters.')

                case 'mvmt2', domvmt2 = 1; disp('Squared/lagged movement model is ON.');

                case {'basis','basisset','bases','basistype'} % name(s) of basis sets to use
                    disp('Found basistype.');
                    basistype = varargin{i+1}; varargin{i+1} = [];

                case 'whichbasis' % enter index in basis sets for each event
                    disp('Found whichbasis.');
                    whichbasis = varargin{i+1};

                case {'hipass','hplength','hp'}, hipass = varargin{i+1};
                    disp('Found high-pass filter length.')

                otherwise
                    disp('Unknown string input.')
            end
        end
    end

        % ---------------------------------------------------------------
        % set up images and onsets for this subject
        % ---------------------------------------------------------------

        ons = onsets{s};    %onsets for subj

        % Tor added this: we now need to enter ntrials to
        % whole_brain_filter
        ntrials = length(ons);

        % ---------------------------------------------------------------
        % CONSTRUCT TRIAL-LEVEL DESIGN MATRIX (models single trials)
        % ---------------------------------------------------------------

        disp('  Constructing design matrix...');

        if ~exist('whichbasis', 'var')
            fprintf('Using the same basis set for all events : %s\n', basistype);

            [trialX, trial_px, bf] = trial_level_beta3('onsets',ons,'output','design', ...
                'rows',sum(nscans),'basistype',basistype,'tr',TR);

            trialX = trialX(:, 1:end-1);    %remove 'experiment' intercept

        else
            fprintf('   Using basis sets indexed by whichbasis for each event\n');
            fprintf('   This will ONLY work if all basis sets have same # parameters.\n');

            trialX = [];

            for i = 1:length(ons)

                fprintf(1,'%3.0f ', i);

                if i == 1
                    [trialX_thistrial, trial_px, bf] = trial_level_beta3('onsets',ons(i),'output','design', ...
                    'rows',sum(nscans),'basistype',basistype{whichbasis(i)},'tr',TR);
                else
                    [trialX_thistrial] = trial_level_beta3('onsets',ons(i),'output','design', ...
                    'rows',sum(nscans),'basistype',basistype{whichbasis(i)},'tr',TR);
                end

                trialX_thistrial = trialX_thistrial(:,1:end-1); % get rid of intercept

                trialX = [trialX trialX_thistrial];

            end

            fprintf(1, '\n');
        end

        fprintf(1,'   Number of basis functions: %3.0f\n',size(bf,2))
        disp('NOTE: all basis sets for different event types MUST have same number of basis functions!')

        %CONSTRUCT HIGH-PASS FILTERING MATRIX
        % ---------------------------------------------------------------
        disp('  Adding high-pass filter to design matrix...');
        khstr = [];
        for r = 1:nruns

            [dummy,KL,KH{r}] = use_spm_filter(TR,nscans(r),'none','specify',hipass);
            if r==1; kh = ['KH{',num2str(r),'}']; else kh = [',KH{',num2str(r),'}']; end
            khstr = strcat(khstr,kh);

        end
        KH = eval(['blkdiag(',khstr,')']);

        % Add movement param stuff here
        % ---------------------------------------------------------------
        MV = [];
        if exist('mvmt', 'var')

            disp('  Adding movement parameters...');

            if size(mvmt, 1) ~= size(trialX, 1)
                error('mvmt params and trial design are different lengths.');
            end

            mvmt = scale(mvmt);

            if domvmt2
                disp('   Adding lagged and squared movement parameters...');
                m2 = [zeros(1, 6); mvmt(1:end-1,:)];
                %gradient(mvmt);
                MV = [mvmt mvmt .^ 2 m2 m2 .^ 2];

            else
                MV = mvmt;
            end

        end

        % Add intercepts and dummy covs for the first 2 trials
        % ---------------------------------------------------------------
        disp('  Adding intercepts and dummy covs for first 1 trials...');
        whremove = 1;  % could be 1:2, 1:3, etc.
        I = intercept_model(nscans, 1, whremove);

        % Put all regressors together
        % ---------------------------------------------------------------
        trialX_of_interest = trialX;

        trialX = [trialX KH MV I];
        trialpx = pinv(trialX);

        % save
        % ---------------------------------------------------------------
        disp(['  Saving to: ',EXPT.subjects{s},'_design.mat']);
        save([EXPT.subjects{s},'_design'], 'trialX', 'trialpx', 'trialX_of_interest');

        % ---------------------------------------------------------------
        % Get info for SUBJECT-LEVEL DESIGN MATRIX (models condition
        % effects across trials)
        % ---------------------------------------------------------------

        % Each row is a trial, each col. is a task condition
        % indicator matrix of 1/0, 1's in col. k where trials belong to condition k

        % Build subject-level design matrix, for analysis of condition
        % effects across trials
        % No intercept: interpret change from abs. baseline of zero

        subjlevelX = eventdesign{s};

        if size(subjlevelX, 1) ~= ntrials
            error('eventdesign is wrong size for ntrials.  Error in either onsets2singletrial or ons.');
        end

        subjlevelpx = pinv(subjlevelX);

        % contrasts should be in column vectors
        if isrow(conditioncontrasts), conditioncontrasts = conditioncontrasts'; end

        % stuff needed for whole_brain_filter trial-level models
        % ---------------------------------------------------------------

        trialmodels.trialX = trialX;
        trialmodels.bf = bf;
        trialmodels.ntrials = ntrials;
        trialmodels.trialpx = trialpx;
        trialmodels.subjlevelX = subjlevelX;
        trialmodels.subjlevelpx = subjlevelpx;
        trialmodels.conditioncontrasts = conditioncontrasts;
        trialmodels.conditionnames = conditionnames;

        if exist('ratings','var')
            disp('Found ratings and will include in trialmodels var, which is used in model fitting.');
            disp('So: will analyze by ratings in model fitting.');

            trialmodels.ratings = ratings{s};
        end

        % stuff to save even tho we don't need it in fitting.
        trialmodels.ons = ons;
        trialmodels.basistype = basistype;
        trialmodels.ntrials = ntrials;

        save([EXPT.subjects{s},'_design'],'-append','trialmodels');

        % tor added this to make image;
        figure; imagesc(trialX), colormap gray; drawnow

disp('Done setting up models.');

    end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%  generate a pair of nxn integer matrices that are inverse to each other
%
%  syntax  >> [X,Y] = IntegerInverse(n)
%          >> [X,Y] = IntegerInverse(n,m)  % entries in interval [-m,m]

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function V = hogDraw( H, w, fhog )
% Create visualization of hog descriptor.
%
% USAGE
%  V = hogDraw( H, [w], [fhog] )
%
% INPUTS
%  H          - [m n oBin*4] computed hog features
%  w          - [15] width for each glyph
%  fhog       - [0] if true draw features returned by fhog
%
% OUTPUTS
%  V          - [m*w n*w] visualization of hog features
%
% EXAMPLE
%
% See also hog, fhog
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.23
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% fold normalizations
if(nargin<3 || isempty(fhog)), fhog=0; end
m=size(H,3); if(fhog), m=(m-5)/3; H=H(:,:,1:m*3); m=3; else m=4; end
s=size(H); s(3)=s(3)/m; w0=H; H=zeros(s);
for o=0:m-1, H=H+w0(:,:,(1:s(3))+o*s(3)); end;

% construct a "glyph" for each orientaion
if(nargin<2 || isempty(w)), w=15; end
bar=zeros(w,w); bar(:,round(.45*w):round(.55*w))=1;
bars=zeros([size(bar) s(3)]);
for o=1:s(3), bars(:,:,o)=imrotate(bar,-(o-1)*180/s(3),'crop'); end

% make pictures of positive weights by adding up weighted glyphs
H(H<0)=0; V=zeros(w*s(1:2));
for r=1:s(1), rs=(1:w)+(r-1)*w;
  for c=1:s(2), cs=(1:w)+(c-1)*w;
    for o=1:s(3), V(rs,cs)=V(rs,cs)+bars(:,:,o)*H(r,c,o); end
  end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function y = UpDyadLo(x,qmf)
% UpDyadLo -- Lo-Pass Upsampling operator; periodized
%  Usage
%    u = UpDyadLo(d,f)
%  Inputs
%    d    1-d signal at coarser scale
%    f    filter
%  Outputs
%    u    1-d signal at finer scale
%
%  See Also
%    DownDyadLo, DownDyadHi, UpDyadHi, IWT_PO, iconv
%
	y =  iconv(qmf, UpSampleN(x) );

%  Revision History
%  10/1/05       AM      UpSample is changed to UpSampleN

%
% Copyright (c) 1993. Iain M. Johnstone
% Last modified on October 2005.
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner


text:
model = 'resnext50' ; % {'cafferef','resnext50','resnext101'}
input = 'rgb' ; % {'rgb','of'}
dataset = 'ucf101' ; % {'ucf101','hmdb51'}  hmdb51 requires more iterations to train (add more epochs to learning rate)
opts.train.batchSize = 128 ;
opts.train.numSubBatches = 32 ; % increase the number (16,32) if it does not fit into gpu mem
opts.epochFactor = 5 ;
opts.split = 1 ;

opts.train.gpus = 1 ;

run matconvnet/matlab/vl_setupnn.m ;
vl_contrib install mcnExtraLayers ; vl_contrib setup mcnExtraLayers ;
vl_contrib install autonn ; vl_contrib setup autonn ;

% addpath(fullfile('matconvnet','contrib','mcnExtraLayers','matlab')) ;

opts.expDir = ['exp/' model 'rgb-arpool-split' num2str(opts.split)] ;
if strcmp(input,'rgb')
  opts.DropOutRate = 0.5 ;
  trainfn = @cnn_dicnn_rgb ;
elseif strcmp(input,'of')
  opts.DropOutRate = 0.8 ;
  trainfn = @cnn_dicnn_of ;
end

if strcmp(model,'cafferef')

  opts.pool1Layer = 'conv1' ;
  % download from http://www.vlfeat.org/matconvnet/models/imagenet-caffe-ref.mat
  opts.modelPath = fullfile('models','imagenet-caffe-ref.mat') ;
  opts.networkFn = @cnn_init_cafferef ;

  if strcmp(input,'rgb')
    opts.train.learningRate = 1e-3 * [ones(1,2) 0.1*ones(1,2)] ;
  else
    opts.train.learningRate = 3e-3 * [ones(1,10) 0.1*ones(1,2)] ;
  end

  opts.train.numEpochs = numel(opts.train.learningRate) ;
elseif strcmp(model,'resnext50') || strcmp(model,'resnext101')
  % download from http://www.robots.ox.ac.uk/~albanie/models/pytorch-imports/resnext_50_32x4d-pt-mcn.mat
  % download from http://www.robots.ox.ac.uk/~albanie/models/pytorch-imports/resnext_101_32x4d-pt-mcn.mat
  if strcmp(model,'resnext50')
    opts.modelPath = fullfile('models','resnext_50_32x4d-pt-mcn.mat') ;
  else
    opts.modelPath = fullfile('models','resnext_101_32x4d-pt-mcn.mat') ;
  end
  opts.modelPath = fullfile('models','resnext_50_32x4d-pt-mcn.mat') ;
  opts.networkFn = @cnn_init_resnext ;
  if strcmp(input,'rgb')
    opts.train.learningRate = 1e-2 * [ones(1,2) 0.1*ones(1,8) ] ;
  else
    opts.train.learningRate = 1e-2 * [ones(1,2) 0.1*ones(1,2) ] ;
  end
end

addpath dicnn ;

[net, info] = trainfn(opts)

author: hbilen

repo: dynamic-image-nets


text:
classdef GemSeparator < GemPanel
    % GemSeparator GEM class for a separator line and label between panels
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        TitleLeftPadding = 5
        TitleHeight = 15
        SeparatorHeight = 20
        TitleWidth = 80
        TitleFontSize = 11
        Bold = true
    end

    properties (Access = protected)
        TitlePanel
    end

    methods
        function obj = GemSeparator(parent_handle, title_text)
            obj = obj@GemPanel(parent_handle);
            obj.TopBorder = true;
            obj.TitlePanel = GemText(obj, upper(title_text), title_text, title_text);
            obj.TitlePanel.FontSize = obj.TitleFontSize;
            obj.TitlePanel.Bold = obj.Bold;
            obj.AddChild(obj.TitlePanel);
        end

        function Resize(obj, position)
            Resize@GemPanel(obj, position);
            title_panel_position = [1 + obj.TitleLeftPadding, 1, position(3) - obj.TitleLeftPadding, obj.TitleHeight];
            obj.TitlePanel.Resize(title_panel_position);
        end

        function height = GetRequestedHeight(obj, width)
            height = obj.SeparatorHeight;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [y,f]=v_zoomfft(x,n,m,s,d)
%V_ZOOMFFT    DTFT evaluated over a linear frequency range Y=(X,N,M,S,D)
% Inputs:
%    x    vector (or matrix)
%    n    reciprocal of normalized frequency increment (can be non-integer).
%         The frequency increment is fs/n where fs is the sample frequency
%         [default n=size(x,d)]
%    m    mumber of output points is floor(m) [default m=n]
%    s    starting frequency index (can be non-integer).
%         The starting frequency is s*fs/n [default s=0]
%    d    dimension along which to do fft [default d=first non-singleton]
%
% Outputs:
%    y       Output dtft coefficients. y has the same dimensions as x except
%            that size(y,d)=floor(m).
%    f(1,m)  normalized frequencies (1 corresponds to fs)
%
% This routine allows the evaluation of the DFT over an arbitrary range of
% frequencies; as its name implies this lets you zoom into a narrow portion
% of the spectrum.
% The DTFT of X will be evaluated along dimension D at the M frequencies
% f=fs*(s+(0:m-1))/n where fs is the sample frequency. Note that N and S
% need not be integers although M will be rounded down to an integer.
% Thus v_zoomfft(x,n,n,0,d) is equivalent to fft(x,n,d) for n>=length(x).

% [1] L.R.Rabiner,  R.W.Schafer and C.M.Rader, "The chirp z-transform algorithm"
%     IEEE Trans. Audio Electroacoustics 17 (2), 86�92 (1969).

%      Copyright (C) Mike Brookes 2007
%      Version: $Id: v_zoomfft.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent n0 k0 s0 m0 b c h g
e=size(x);
p=prod(e);
if nargin<5
    d=find(e>1);
    if ~isempty(d)
        d=d(1);
    else
        d=1;
    end
end
k=e(d);
q=p/k;
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
if nargin<2 || isempty(n)
    n=k;
end
if nargin<3 || isempty(m)
    m=floor(n);
else
    m=floor(m);
end
if nargin<4 || isempty(s)
    s=0;
end
l=pow2(nextpow2(m+k-1));    % round up to next power of 2
if n==fix(n) && s==fix(s) && n<2*l && n>=k
    a=fft(z,n,1);           % quickest to do a normal fft
    y=a(1+mod(s:s+m-1,n),:);
else
    % can precaluclate all this for fixed n, k, s and m
    if isempty(b) || n~=n0 || k~=k0 || s~=s0 || m~=m0
        n0=n;
        k0=k;
        s0=s;
        m0=m;
        b=exp(1i*pi*mod((s+(1-k:m-1)').^2,2*n)/n);
        c=conj(b(k:k+m-1));
        h=fft(b,l,1);
        g=exp(-1i*pi*mod(((0:k-1)').^2,2*n)/n);
    end
    a=ifft(fft(z.*repmat(g,1,q),l,1).*repmat(h,1,q)); % calculate correlation
    y=a(k:k+m-1,:).*repmat(c,1,q);
end
if d==1
    e(d)=m;
    y=reshape(y,e);
else
    r(1)=m;
    y=shiftdim(reshape(y,r),length(e)+1-d);
end
f=(s+(0:m-1))/n;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% This class is for computing superpixel-based affinities described in the paper.
% This class requires the image graphs methods by Steve Eddins. Find it here:
% http://www.mathworks.com/matlabcentral/fileexchange/53614-image-graphs

classdef Superpixels
    properties
        labels
        spcount
        neigh
        centroids
    end
    methods
        function obj = Superpixels(im, spcnt)
            if ~exist('spcnt', 'var') || isempty(spcnt)
                spcnt = 2500;
            end
            [L, N] = superpixels(im, spcnt, 'Compactness', 1e-20);
            obj.labels = L;
            obj.spcount = N;
            % Find neighboring superpixels
            g = adjacentRegionsGraph(L);
            obj.neigh = g.Edges.Labels;
            % Find centroids
            s = regionprops(L, 'centroid');
            cent = cat(1, s.Centroid);
            obj.centroids = round(cent(:, 2:-1:1));
            [h, w, ~] = size(im);
            obj.centroids(:, 3) = sub2ind([h, w], obj.centroids(:, 1), obj.centroids(:, 2));
        end

        function regmeans = computeRegionMeans(obj, image)
            [h, w, c] = size(image);
            image = reshape(image, [h*w, c]);
            regmeans = zeros(obj.spcount, c);
            idx = label2idx(obj.labels);
            for i = 1 : length(idx)
                regmeans(i, :) = mean(image(idx{i}, :), 1);
            end
        end

        % This is for the semantic affinity, generates affinities in [-1, 1]
        function W = neighborAffinities(obj, features, erfSteepness, erfCenter)
            if ~exist('erfSteepness', 'var') || isempty(erfSteepness)
                erfSteepness = 20;
            end
            if ~exist('erfCenter', 'var') || isempty(erfCenter)
                erfCenter = 0.85;
            end
            [h, w, ~] = size(features);
            N = h * w;
            spMeans = obj.computeRegionMeans(features);
            affs = zeros(size(obj.neigh, 1), 1);
            inds1 = affs;
            inds2 = affs;
            for i = 1 : size(obj.neigh, 1)
                ind1 = obj.neigh(i, 1);
                ind2 = obj.neigh(i, 2);
                affs(i) = sigmoidAff(spMeans(ind1, :), spMeans(ind2, :), erfSteepness, erfCenter);
                inds1(i) = obj.centroids(ind1, 3);
                inds2(i) = obj.centroids(ind2, 3);
            end
            W = sparse(inds1, inds2, affs, N, N);
            W = W' + W;
        end

        % This is for the nonlocal color affinity, generates affinities in [0, 1]
        function W = nearbyAffinities(obj, image, erfSteepness, erfCenter, proxThresh)
            if ~exist('erfSteepness', 'var') || isempty(erfSteepness)
                erfSteepness = 50;
            end
            if ~exist('erfCenter', 'var') || isempty(erfCenter)
                erfCenter = 0.95;
            end
            if ~exist('proxThresh', 'var') || isempty(proxThresh)
                proxThresh = 0.2;
            end
            [h, w, ~] = size(image);
            N = h * w;
            spMeans = obj.computeRegionMeans(image);
            combinationCnt = obj.spcount;
            combinationCnt = combinationCnt * (combinationCnt - 1) / 2;
            affs = zeros(combinationCnt, 1);
            inds1 = affs;
            inds2 = affs;
            cnt = 1;
            cents = obj.centroids(:, 1:2);
            cents(:,1) = cents(:,1) / h;
            cents(:,2) = cents(:,2) / w;
            for i = 1 : obj.spcount
                for j = i + 1 : obj.spcount
                    centdist = cents(i, 1:2) - cents(j, 1:2);
                    centdist = sqrt(centdist * centdist');
                    if centdist > proxThresh
                        affs(cnt) = 0;
                    else
                        affs(cnt) = sigmoidAffPos(spMeans(i, :), spMeans(j, :), erfSteepness, erfCenter);
                    end
                    inds1(cnt) = obj.centroids(i, 3);
                    inds2(cnt) = obj.centroids(j, 3);
                    cnt = cnt + 1;
                end
            end
            W = sparse(inds1, inds2, affs, N, N);
            W = W' + W;
        end

        function vis = visualizeRegionMeans(obj, im)
            vis = label2rgb(obj.labels, obj.computeRegionMeans(im));
        end

    end
end

function aff = sigmoidAff(feat1, feat2, steepness, center)
    aff = abs(feat1 - feat2);
    aff = 1 - sqrt(aff * aff');
    aff = (erf(steepness * (aff - center)));
end

function aff = sigmoidAffPos(feat1, feat2, steepness, center)
    aff = abs(feat1 - feat2);
    aff = 1 - sqrt(aff * aff');
    aff = (erf(steepness * (aff - center)) + 1) / 2;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function y=design(x)
% DESIGN
% DESIGN(x) creates a matrix of 1 and 0 corresponding to x
% y: Rows(x) x Max(x)
  i=1:1:rows(x);
  k=ones(rows(x),1);
  s=sparse(i,x,k,rows(x),max(x));
  y=full(s);

author: Sable

repo: mcbench-benchmarks


text:
function [MHz] = rpm2MHz(rpm)
% Convert frequency from revolutions per minute to Megahertz.
% Chad A. Greene 2012
MHz = rpm/60/1000000;


text:
function dh=constDhDataOf14R820()
%% This function returns the DH parameters
%     for the KUKA iiwa 7R800

% The dimnstions are taken from Kuka manuals:
% [1] Medien-Flansch
%     Für Produktfamilie LBR iiwa
%     Montage- und Betriebsanleitung

%% DH PARAMETERS FOR THE ROBOT
dh.alfa={0,-pi/2,pi/2,pi/2,-pi/2,-pi/2,pi/2};
% following are "d" parameters for iiwa 14 R 820
dh.d={0.36,0.0,0.42,0.0,0.4,0.0,0.126};
dh.a={0,0,0,0,0,0,0};
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [MPa] = Torr2MPa(Torr)
% Convert pressure from torr (same as mmHg) to megapascals
% Chad Greene 2012
MPa = Torr*0.000133322;