text:
function [ n_data, h, a, t ] = owen_values ( n_data )

%*****************************************************************************80
%
%% OWEN_VALUES returns some values of Owen's T function.
%
%  Discussion:
%
%    Owen's T function is useful for computation of the bivariate normal
%    distribution and the distribution of a skewed normal distribution.
%
%    Although it was originally formulated in terms of the bivariate
%    normal function, the function can be defined more directly as
%
%      T(H,A) = 1 / ( 2 * pi ) *
%        Integral ( 0 <= X <= A ) e^(H^2*(1+X^2)/2) / (1+X^2) dX
%
%    In Mathematica, the function can be evaluated by:
%
%      fx = 1/(2*Pi) * Integrate [ E^(-h^2*(1+x^2)/2)/(1+x^2), {x,0,a} ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real H, a parameter.
%
%    Output, real A, the upper limit of the integral.
%
%    Output, real T, the value of the function.
%
  n_max = 22;

  a_vec = [ ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.1000000000000000E+02, ...
    0.1000000000000000E+03 ];

  h_vec = [ ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02 ];

  t_vec = [ ...
    0.4306469112078537E-01, ...
    0.6674188216570097E-01, ...
    0.7846818699308410E-01, ...
    0.7929950474887259E-01, ...
    0.6448860284750376E-01, ...
    0.1066710629614485E+00, ...
    0.1415806036539784E+00, ...
    0.1510840430760184E+00, ...
    0.7134663382271778E-01, ...
    0.1201285306350883E+00, ...
    0.1666128410939293E+00, ...
    0.1847501847929859E+00, ...
    0.7317273327500385E-01, ...
    0.1237630544953746E+00, ...
    0.1737438887583106E+00, ...
    0.1951190307092811E+00, ...
    0.7378938035365546E-01, ...
    0.1249951430754052E+00, ...
    0.1761984774738108E+00, ...
    0.1987772386442824E+00, ...
    0.2340886964802671E+00, ...
    0.2479460829231492E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    h = 0.0;
    a = 0.0;
    t = 0.0;
  else
    h = h_vec(n_data);
    a = a_vec(n_data);
    t = t_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function parsave(filename, structX)
    vl_xmkdir(fileparts(filename));
    save(filename,'-struct','structX');
end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
%==================================================================
%                        General Data File
% Title: Default_title
% Units: SI
% Dimensions: 2D
% Type of problem: Plane_Stress
% Type of Phisics: Stokes
% Micro/Macro: MACRO
%
%==================================================================

%% Data
order={'QUADRATIC','LINEAR'};
Data_prb = {
'TRIANGLE';
'SI';
'2D';
'Plane_Stress';
'Stokes';
'MACRO'
};
state = 'Transient';
dtime = 0.01;
finalTime = 1;

%% Coordinates
% Node                X                Y                Z

coord = [
1            0            0            0
2            0         0.25            0
3         0.25            0            0
4         0.25         0.25            0
5            0          0.5            0
6          0.5            0            0
7         0.25          0.5            0
8          0.5         0.25            0
9          0.5          0.5            0
10         0.75            0            0
11            0         0.75            0
12         0.75         0.25            0
13         0.25         0.75            0
14         0.75          0.5            0
15          0.5         0.75            0
16            1            0            0
17            0            1            0
18            1         0.25            0
19         0.25            1            0
20         0.75         0.75            0
21            1          0.5            0
22          0.5            1            0
23         0.75            1            0
24            1         0.75            0
25            1            1            0
];

%% Conectivities
% Element        Node(1)                Node(2)                Node(3)                Material

connec = [
1 6 8 3 0
2 8 9 4 0
3 3 4 1 0
4 8 4 3 0
5 9 7 4 0
6 7 5 2 0
7 4 2 1 0
8 7 2 4 0
9 16 18 10 0
10 18 21 12 0
11 10 12 6 0
12 18 12 10 0
13 21 14 12 0
14 14 9 8 0
15 12 8 6 0
16 14 8 12 0
17 9 15 7 0
18 15 22 13 0
19 7 13 5 0
20 15 13 7 0
21 22 19 13 0
22 19 17 11 0
23 13 11 5 0
24 19 11 13 0
25 21 24 14 0
26 24 25 20 0
27 14 20 9 0
28 24 20 14 0
29 25 23 20 0
30 23 22 15 0
31 20 15 9 0
32 23 15 20 0
];

%% Variable Prescribed
% Node            Dimension                Value

velocityBC.domain = 'Border';
velocityBC.value  = 0;

velocity = [
1 1 0
1 2 0
5 1 0
5 2 0
6 1 0
6 2 0
16 1 0
16 2 0
17 1 0
17 2 0
21 1 0
21 2 0
22 1 0
22 2 0
25 1 0
25 2 0
];

pressure = [
25 1 0
];

nu=1;
Vol_force = @(x,y){nu*4*(x^3*((6 - 12*y)) + x^4*((-3 + 6*y)) + y*(1 - 3*y + 2*y^2)-6*x*y*(1 - 3*y + 2*y^2)...
+ 3*x^2*((-1 + 4*y - 6*y^2 + 4*y^3)))+(2*x-1)*(y-1);
-4*nu*(-3*((-1 + y)^2)*y^2 - 3*x^2*((1 - 6*y + 6*y^2))+2*x^3*((1 - 6*y + 6*y^2)) +...
x*(1 - 6*y + 12*y^2 - 12*y^3 + 6*y^4))+x*(x-1)};

% nu=1;
% Vol_force = @(x,y){nu*4*(x^3*((6 - 12*y)) + x^4*((-3 + 6*y)) + y*(1 - 3*y + 2*y^2)-6*x*y*(1 - 3*y + 2*y^2)...
% + 3*x^2*((-1 + 4*y - 6*y^2 + 4*y^3)))+nu*y*(1 - y)*(1 - 2*x);
% -4*nu*(-3*((-1 + y)^2)*y^2 - 3*x^2*((1 - 6*y + 6*y^2))+2*x^3*((1 - 6*y + 6*y^2)) +...
% x*(1 - 6*y + 12*y^2 - 12*y^3 + 6*y^4)+x*(1 - x)*(1 - 2*y))};

%% Group Elements
% Element        Group_num

Group = [
];

%% Initial Holes
% Elements that are considered holes initially
% Element

Initial_holes = [
];

%% Boundary Elements
% Elements that can not be removed
% Element

Boundary_elements = [
];

%% Micro gauss post
%
% Element

Micro_gauss_post = [
];

%% Micro Slave-Master
% Nodes that are Slaves
% Nodes             Value (1-Slave,0-Master)

Micro_slave = [
];

%% Nodes solid
% Nodes that must remain
% Nodes

% nodesolid = unique(pointload_complete(:,1));

%% External border Elements
% Detect the elements that define the edge of the domain
% Element               Node(1)           Node(2)

External_border_elements = [
1 3 6
3 1 3
6 5 2
7 2 1
9 16 18
9 10 16
10 18 21
11 6 10
21 22 19
22 17 11
22 19 17
23 11 5
25 21 24
26 24 25
29 25 23
30 23 22
];

%% External border Nodes
% Detect the nodes that define the edge of the domain
% Node

External_border_nodes = [
];

%% Materials
% Materials that have been used
% Material_Num              Mat_density        Young_Modulus        Poisson

Materials = [
];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%Compile pCT radiography

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Stefanie Kaser, Benjamin Kirchmayer
%--------------------------------------------------------------------------
clear all;

%% Precompile actions to write intercepts vector lengths
prompt = 'Enter pixel side length in mm (or hit enter): ';
input0 = input(prompt);

if input0 > 0
    pix = single(input0/0.25);

    prompt = 'Enter SL intercepts before phantom (or hit enter for default values): ';
    input1 = input(prompt);
    prompt = 'Enter CS intercepts within hull (or hit enter for default values): ';
    input2 = input(prompt);
    prompt = 'Enter SL intercepts behind phantom (or hit enter for default values): ';
    input3 = input(prompt);

    try
        if input1 > 0
            interceps_in = uint16(input1);
        else
            interceps_in = uint16(10/pix);
            disp(strcat('Setting intercepts vector size (SL, in) to default: ', int2str(uint16(10/pix))));
        end
    catch
        interceps_in = uint16(10/pix);
        disp(strcat('Setting intercepts vector size (SL, in) to default: ', int2str(uint16(10/pix))));
    end

    try
        if input2 > 0
            interceps_cs = uint16(input2);
        else
            interceps_cs = uint16(220/pix);
            disp(strcat('Setting intercepts vector size (CS within hull) to default: ', int2str(uint16(220/pix))));
        end
    catch
        interceps_cs = uint16(220/pix);
        disp(strcat('Setting intercepts vector size (CS within hull) to default: ', int2str(uint16(220/pix))));
    end

    try
        if input3 > 0
            interceps_out = uint16(input3);
        else
            interceps_out = 100;
            disp(strcat('Setting intercepts vector size (SL, out) to default: ', int2str(uint16(100/pix))));
        end
    catch
        interceps_out = 100;
        disp(strcat('Setting intercepts vector size (SL, out) to default: ', int2str(uint16(100/pix))));
    end
else
    warning('No valid pixel size was entered. Intercepts vector size can be set manually or default values can be used. Be aware that default values are optimal for a pixel size of 0.25x0.25mm^2.')
    prompt = 'Enter SL intercepts before phantom (or hit enter for default values): ';
    input1 = input(prompt);
    prompt = 'Enter CS intercepts within hull (or hit enter for default values): ';
    input2 = input(prompt);
    prompt = 'Enter SL intercepts behind phantom (or hit enter for default values): ';
    input3 = input(prompt);

    try
        if input1 > 0
            interceps_in = uint16(input1);
        else
            interceps_in = 10;
            disp('Setting intercepts vector size (SL, in) to default (10)');
        end
    catch
        interceps_in = 10;
        disp('Setting intercepts vector size (SL, in) to default (10)');
    end

    try
        if input2 > 0
            interceps_cs = uint16(input2);
        else
            interceps_cs = 220;
            disp('Setting intercepts vector size (CS within hull) to default (220).');
        end
    catch
        interceps_cs = 220;
        disp('Setting intercepts vector size (CS within hull) to default (220).');
    end

    try
        if input3 > 0
            interceps_out = uint16(input3);
        else
            interceps_out = 100;
            disp('Setting intercepts vector size (SL, out) to default (100).');
        end
    catch
        interceps_out = 100;
        disp('Setting intercepts vector size (SL, out) to default (100).');
    end
end

% SL (before phantom)
try
    textCell = readlines('../Common/CUDA/improvedForwardProjections.hpp');
catch
   error('Assure to start the compile script from the correct directory')
end
searchMask = cell2mat(cellfun(@(x) contains(x, '#define vecSizeIn'), textCell','UniformOutput', false))';
textCell{searchMask} = ['#define vecSizeIn ' num2str(interceps_in)];

writeID = fopen('../Common/CUDA/improvedForwardProjections.hpp', 'w');

for i=1:numel(textCell)
    fprintf(writeID, '%s\n', textCell{i});
end
fclose(writeID);

% CS (within hull)
textCell = readlines('../Common/CUDA/improvedForwardProjections.hpp');
searchMask = cell2mat(cellfun(@(x) contains(x, '#define vecSizeCS'), textCell','UniformOutput', false))';
textCell{searchMask} = ['#define vecSizeCS ' num2str(interceps_cs)];

writeID = fopen('../Common/CUDA/improvedForwardProjections.hpp', 'w');

for i=1:numel(textCell)
    fprintf(writeID, '%s\n', textCell{i});
end
fclose(writeID);

% SL (behind phantom)
textCell = readlines('../Common/CUDA/improvedForwardProjections.hpp');
searchMask = cell2mat(cellfun(@(x) contains(x, '#define vecSizeOut'), textCell','UniformOutput', false))';
textCell{searchMask} = ['#define vecSizeOut ' num2str(interceps_out)];

writeID = fopen('../Common/CUDA/improvedForwardProjections.hpp', 'w');

for i=1:numel(textCell)
    fprintf(writeID, '%s\n', textCell{i});
end
fclose(writeID);

%% Compile
clear all;
addpath('pCTMexFiles');
addpath('./Utilities/Setup');
mex -setup
[cudapath, cuda_ver]=locate_cuda();
if isempty(cudapath)
    error(sprintf('CUDA Path not found. \nAdd the path by writting in MATLAB:\nsetenv(''CUDA_PATH'',''your path'')\nWhere "your path" is C:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v11.2, for example, \nor /usr/local/cuda on linux')) ;
end
if ispc
   setenv('CUDA_PATH',cudapath);
end
set_cuda_cc_flags(cuda_ver);

rmpath('./Utilities/Setup');
% if isempty(getenv('CUDA_PATH'))
%    setenv('CUDA_PATH','usr/local/cuda') ;
% end

fprintf("Compiling pCT source...\n");
if isunix
    if ~isempty(strfind(computer('arch'),'64'))
        mex -largeArrayDims ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/linux64
    else
        mex ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/linux32
    end
elseif ispc
    if ~isempty(strfind(computer('arch'),'64'))
        mex -largeArrayDims ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/win64
    else
        mex ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/win32
    end
elseif ismac
    if ~isempty(strfind(computer('arch'),'64'))
        mex -largeArrayDims ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/mac64
    else
        mex ./Utilities/cuda_interface/pCTCubicSpline_mex.cpp ../Common/CUDA/improvedForwardProjections.cu ../Common/CUDA/improvedForwardProjections_cone.cu -outdir ./pCTMexFiles/mac32
    end
end

if ispc
    if ~isempty(strfind(computer('arch'),'64'))
        addpath('./pCTMexFiles/win64');
    else
        addpath('./pCTMexFiles/win32');
    end
elseif ismac
    if ~isempty(strfind(computer('arch'),'64'))
        addpath('./pCTMexFiles/mac64');
    else
        addpath('./pCTMexFiles/mac32');
    end
else
    if ~isempty(strfind(computer('arch'),'64'))
        addpath('./pCTMexFiles/linux64');
    else
        addpath('./pCTMexFiles/linux32');
    end
end

fprintf('Compilation successful!\n');

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function outpoints = mni2tal(inpoints)
% Converts coordinates from MNI brain to best guess
% for equivalent Talairach coordinates
% FORMAT outpoints = mni2tal(inpoints)
% Where inpoints is N by 3 or 3 by N matrix of coordinates
%  (N being the number of points)
% outpoints is the coordinate matrix with Talairach points
% Matthew Brett 10/8/99

% ensure that SPM is available, needed for spm_matrix
hasspm = ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

dimdim = find(size(inpoints) == 3);
if isempty(dimdim)
  ft_error('input must be a N by 3 or 3 by N matrix')
end
if dimdim == 2
  inpoints = inpoints';
end

% Transformation matrices, different zooms above/below AC
upT   = spm_matrix([0 0 0 0.05 0 0 0.99 0.97 0.92]);
downT = spm_matrix([0 0 0 0.05 0 0 0.99 0.97 0.84]);

tmp = inpoints(3,:)<0;  % 1 if below AC
inpoints = [inpoints; ones(1, size(inpoints, 2))];
inpoints(:, tmp) = downT * inpoints(:, tmp);
inpoints(:, ~tmp) = upT * inpoints(:, ~tmp);
outpoints = inpoints(1:3, :);
if dimdim == 2
  outpoints = outpoints';
end

author: spm


text:
function x = hatmap(u)
%HATMAP Summary of this function goes here
%   return skew symmetric matrix from a vector
x = [0, -u(3), u(2);
    u(3), 0, -u(1);
    -u(2), u(1), 0];
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function [signature] = getAertsSign_DICOMconvention(vol,mask,spatialRef,scanType)

% PARAMETERS
if strcmp(scanType,'CT')
    res = [1,1,3]; % In (Aerts,2014), all features were calculated from 1 X 1 X 3 mm^3 raw images
    levels = 0:25:3000; % Division of quantized bins in (Aerts,2014) for texture analysis
elseif strcmp(scanType,'PET') % Novelty tested in this work
    levels = 0:1:100; % Division of quantized bins for texture analysis, similarly to (Aerts,2014)
    res = [4,4,4]; % Resolution is here set up to the approximate resolution of PET scans
end

% INITIALIZATION
signature = zeros(1,4);
nLevels = numel(levels) - 1;

% STEP 0: RESAMPLING OF VOLUME AND MASK
pixelW = spatialRef.PixelExtentInWorldX;
sliceS = spatialRef.PixelExtentInWorldZ;
a = pixelW/res(1); b = pixelW/res(2); c = sliceS/res(3);
maskBox = imresize3D(mask,[],[round(double(size(mask,1))*a),round(double(size(mask,2))*b),round(double(size(mask,3))*c)],'nearest','fill');
ROIbox = imresize3D(vol,[],[round(double(size(vol,1))*a),round(double(size(vol,2))*b),round(double(size(vol,3))*c)],'cubic','fill');
ROIonly = ROIbox; ROIonly(~maskBox) = NaN;

% STEP 1: CALCULATE FEATURE 1
energy = ROIonly .^2; energy = sum(energy(~isnan(energy)));
signature(1) = energy;

% STEP 2: CALCULATE FEATURE 2
maskArea = imfill(maskBox,'holes'); maskArea = padarray(maskArea,[1,1,1],0); % Not using inner surfaces
volume = sum(maskBox(:)) * prod(res); % In mm^3
[p,q,r] = meshgrid(res(1).*(1:1:size(maskArea,2)),res(2).*(1:1:size(maskArea,1)),res(3).*(1:1:size(maskArea,3)));
[faces,vertices] = isosurface(p,q,r,maskArea,.5);
a = vertices(faces(:,2),:) - vertices(faces(:,1),:);
b = vertices(faces(:,3),:) - vertices(faces(:,1),:);
c = cross(a,b,2);
surface = 1/2 * sum(sqrt(sum(c.^2,2)));
compactness = volume/(sqrt(pi)*surface^(2/3));
signature(2) = compactness;

% STEP 3: CALCULATE FEATURE 3
ROIquant = ROIonly;
ROIquant(ROIonly < levels(1)) = 1;
for n = 1:nLevels
    ROIquant(ROIonly>=levels(n) & ROIonly<levels(n+1)) = n;
end
ROIquant(ROIonly >= levels(end)) = nLevels;
newLevels = 1:max(ROIquant(~isnan(ROIquant)));
[GLN] = getGLN_Aerts(ROIquant,newLevels);
signature(3) = GLN;

% STEP 4: CALCULATE FEATURE 4 (HLH --> 1) H: DICOM = left-right, Axial image = left-right; 2) L: DICOM = anterior posterior, Axial image = top-bottom; 3) H: DICOM = superior-inferior, Axial image = 3rd dimension)
% - DICOM XYZ coordinates are assumed in HLH feature definition of (Aerts et al, 2014)
% - From the swt2 function, we have in DICOM XYZ coordinates: [LL,LH,HL,HH] = swt2(image,N,'wavName')
sizeV = size(ROIbox);
if mod(sizeV(1),2)
    sizeV(1) = sizeV(1) + 1;
end
if mod(sizeV(2),2)
    sizeV(2) = sizeV(2) + 1;
end
if mod(sizeV(3),2)
    sizeV(3) = sizeV(3) + 1;
end
HLH = zeros(sizeV);
maskBox = imresize3D(maskBox,[],sizeV,'nearest','fill'); % Making sure the dimensions are even for decomposition at level 1
ROIbox = imresize3D(ROIbox,[],sizeV,'cubic','fill');
for k = 1:sizeV(3) % In this image plane (axial), we need the HL sub-band
    image = ROIbox(:,:,k);
    [~,~,HL,~] = swt2(image,1,'coif1');
    HLH(:,:,k) = HLH(:,:,k) + HL;
end
for i = 1:sizeV(1) % In this image plane (coronal), we need the HH sub-band
    image = squeeze(ROIbox(i,:,:))';
    [~,~,~,HH] = swt2(image,1,'coif1');
    for k = 1:sizeV(3)
        HLH(i,:,k) = HLH(i,:,k) + HH(k,:);
    end
end
for j = 1:sizeV(2) % In this image plane (sagittal), we need the HL sub-band
    image = squeeze(ROIbox(:,j,:));
    [~,~,HL,~] = swt2(image,1,'coif1');
    for k = 1:sizeV(3)
        HLH(:,j,k) = HLH(:,j,k) + HL(:,k);
    end
end
HLH = HLH/3;

ROIquant = HLH; ROIquant(~maskBox) = NaN;
Ng = numel(newLevels); [ROIquant,levels] = uniformQuantization(ROIquant,Ng);
[GLN] = getGLN_Aerts(ROIquant,levels);
signature(4) = GLN;

end

author: mvallieres

repo: radiomics


text:
function hra = computeHourAngle(ut, long, bodyInfo)
%     if(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.TwoBody)
%         chain = bodyInfo.getOrbitElemsChain();
%         rVectBodyToSun = -1.0 * getPositOfBodyWRTSun_alg(ut, chain{:});
%
% 	elseif(bodyInfo.propTypeEnum == BodyPropagationTypeEnum.Numerical)
%
%     else
%
%     end
    [rVectBodyToSun, ~] = getPositOfBodyWRTSun(ut, bodyInfo, bodyInfo.celBodyData);
    rVectBodyToSun = -1.0 * rVectBodyToSun;

    rVectBodyToSun = reshape(rVectBodyToSun, 3,numel(ut));
    rVectBodyToSunNorm = sqrt(sum(rVectBodyToSun.^2,1));

    inputs = bodyInfo.getFixedFrameFromInertialFrameInputsCache();
    rVectSunECEF = getFixedFrameVectFromInertialVect_alg(ut, rVectBodyToSun, inputs{:});

    rECEF = getrVectEcefFromLatLongAlt(zeros(size(long)), long, zeros(size(long)), bodyInfo);

    planarRvectEcef = [rECEF(1,:); rECEF(2,:); zeros(1,size(rECEF,2))];
    planarRvectSunEcef = [rVectSunECEF(1,:); rVectSunECEF(2,:); zeros(1,size(rVectSunECEF,2))];

    hra = angleNegPiToPi(dang(planarRvectEcef,planarRvectSunEcef));
    hra(rVectBodyToSunNorm == 0) = 1.0;
    hra = reshape(hra,size(ut));
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% Color Mapping
%
%%
% A central issue when interpreting plots is to have a consistent color
% coding among all plots. In MTEX this can be achieved in two ways. If the
% the minimum and maximum value are known then one can specify the color
% range directly using the options *colorrange* or *contourf*, or the
% command <setColorRange.html setColorRange> is used which allows to set
% the color range afterwards.
%
%% A sample ODFs and Simulated Pole Figure Data
%
% Let us first define some model <SO3Fun.SO3Fun.html ODFs> to be plotted later
% on.

cs = crystalSymmetry('-3m');
odf = fibreODF(Miller(1,1,0,cs),zvector)
pf = calcPoleFigure(odf,[Miller(1,0,0,cs),Miller(1,1,1,cs)],...
  equispacedS2Grid('points',500,'antipodal'));

%% Tight Colorcoding
%
% When <PoleFigure.plot.html plot> is called without any colorcoding option,
% the plots are constructed using the  *tight* option to the range of the data
% independently from the other plots. This means that different pole
% figures may have different color coding and in principle cannot be
% compared to each other.

close all
plot(pf)
mtexColorbar

%% Equal Colorcoding
%
% The *tight* colorcoding can make the reading and comparison of two pole figures
% a bit hard. If you want to have one colorcoding for all plots within one figure use the
% option *colorrange* to *equal*.

plot(pf,'colorRange','equal')
mtexColorbar

%% Setting an Explicite Colorrange
%
% If you want to have a unified colorcoding for several figures you can
% set the colorrange directly in the <SO3Fun.plotPDF.html plot command>

close all
plotPDF(odf,[Miller(1,0,0,cs),Miller(1,1,1,cs)],...
  'colorrange',[0 4],'antipodal');
mtexColorbar

figure
plotPDF(.5*odf+.5*uniformODF(cs),[Miller(1,0,0,cs),Miller(1,1,1,cs)],...
  'colorrange',[0 4],'antipodal');
mtexColorbar

%% Setting the Contour Levels
%
% In the case of contour plots you can also specify the *contour levels*
% directly

close all
plotPDF(odf,[Miller(1,0,0,cs),Miller(1,1,1,cs)],...
  'contourf',0:1:5,'antipodal')
mtexColorbar

%% Modifying the Colorrange After Plotting
%
% The color range of the figures can also be adjusted afterwards using the
% command <mtexFigure.CLim.html CLim>

CLim(gcm,[0.38,3.9])

%% Logarithmic Plots
%
% Sometimes logarithmic scaled plots are of interest. For this case all
% plots in MTEX understand the option *logarithmic*, e.g.

close all;
plotPDF(odf,[Miller(1,0,0,cs),Miller(1,1,1,cs)],'antipodal','logarithmic')
CLim(gcm,[0.01 12]);
mtexColorbar

%% Changing the Colormap
%
% The colormap can be changed by the command mtexColorMap, e.g., in order
% to set a white to black colormap one has the commands

plotPDF(odf,[Miller(1,0,0,cs),Miller(1,1,1,cs)],'antipodal')
mtexColorMap white2black
mtexColorbar

%% Multiple Colormaps
%
% One can even use different colormaps within one figure

% initialize an MTEXFigure
mtexFig = newMtexFigure;

% for three different colormaps
for cm = {'hot', 'cool', 'parula'}

  % generate a new axis
  nextAxis

  % plot some random data in different axis
  plot(vector3d.rand(100),'smooth','grid','grid_res',90*degree,'upper');

  % and apply an individual colormap
  mtexColorMap(mtexFig.gca,char(cm))

  % set the title to be the name of the colormap
  mtexTitle(char(cm))
end

% plot a colorbar for each plot
mtexColorbar

author: mtex-toolbox

repo: mtex

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


text:
function [flag] = certify(mesh)
%CERTIFY error checking for JIGSAW mesh objects.

%-----------------------------------------------------------
%   Darren Engwirda
%   github.com/dengwirda/jigsaw-matlab
%   07-Aug-2019
%   darren.engwirda@columbia.edu
%-----------------------------------------------------------
%

    np = +0 ; flag = -1;

    if (~isstruct(mesh))
        error('certify:incorrectInputClass', ...
            'Incorrect input class.') ;
    end

    if (inspect(mesh,'point'))
        if (~isempty (mesh.point.coord))
        if (isnumeric(mesh.point.coord))
%----------------------------------------- check MESH coords
        np = size(mesh.point.coord,1) ;

        if (ndims(mesh.point.coord) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid POINT.COORD dimensions.') ;
        end
        if ( size(mesh.point.coord,2)< 3)
        error('certify:incorrectDimensions', ...
            'Invalid POINT.COORD dimensions.') ;
        end

        if (any(isinf(mesh.point.coord)))
        error('certify:invalidMeshPosition', ...
            'Invalid POINT.COORD values.') ;
        end
        if (any(isnan(mesh.point.coord)))
        error('certify:invalidMeshPosition', ...
            'Invalid POINT.COORD values.') ;
        end

        if (isfield(mesh,'mshID'))
        if (strcmpi(mesh.mshID,'euclidean-grid'))
        error('certify:incompatiblemshID', ...
            'Incompatible msh-ID flag.') ;
        end
        if (strcmpi(mesh.mshID,'ellipsoid-grid'))
        error('certify:incompatiblemshID', ...
            'Incompatible msh-ID flag.') ;
        end
        end

        elseif(iscell(mesh.point.coord))
%----------------------------------------- check GRID coords
        if (~isvector(mesh.point.coord))
        error('certify:incorrectDimensions', ...
            'Invalid POINT.COORD dimensions.') ;
        end

        for ii = +1:length(mesh.point.coord)

        if (~isvector(mesh.point.coord{ii}))
        error('certify:incorrectDimensions', ...
            'Invalid POINT.COORD dimensions.') ;
        end
        if (any(isinf(mesh.point.coord{ii})))
        error('certify:invalidMeshPosition', ...
            'Invalid POINT.COORD values.') ;
        end
        if (any(isnan(mesh.point.coord{ii})))
        error('certify:invalidMeshPosition', ...
            'Invalid POINT.COORD values.') ;
        end

        end

        if (isfield(mesh,'mshID'))
        if (strcmpi(mesh.mshID,'euclidean-mesh'))
        error('certify:incompatiblemshID', ...
            'Incompatible msh-ID flag.') ;
        end
        if (strcmpi(mesh.mshID,'ellipsoid-mesh'))
        error('certify:incompatiblemshID', ...
            'Incompatible msh-ID flag.') ;
        end
        end

        else
%----------------------------------------- wrong POINT class
        error('certify:incorrectInputClass', ...
            'Invalid POINT.COORD type.') ;

        end
        end
    end

    if (inspect(mesh,'radii'))
%----------------------------------------- check RADII value
        if (~isempty  (mesh.radii))
        if (~isnumeric(mesh.radii))
        error('certify:incorrectInputClass', ...
            'Invalid RADII class.') ;
        end
        if (~isvector(mesh.radii))
        error('certify:incorrectDimensions', ...
            'Invalid RADII dimensions.') ;
        end
        if (length(mesh.radii) ~= +1 && ...
            length(mesh.radii) ~= +3 )
        error('certify:incorrectDimensions', ...
            'Invalid RADII dimensions.') ;
        end
        if (any(isinf(mesh.radii)))
        error('certify:invalidRadiiEntries', ...
            'Invalid RADII entries.') ;
        end
        if (any(isnan(mesh.radii)))
        error('certify:invalidRadiiEntries', ...
            'Invalid RADII entries.') ;
        end
        end
    end

    if (inspect(mesh,'value'))
%----------------------------------------- check VALUE value
        if (~isempty  (mesh.value))
        if (isnumeric(mesh.point.coord))
%----------------------------------------- for MESH obj kind

        if (ndims(mesh.value) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end
        if ( size(mesh.value,1) ~= np)
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end

        if (any(isinf(mesh.value)))
        error('certify:invalidValueEntries', ...
            'Invalid VALUE entries.') ;
        end
        if (any(isnan(mesh.value)))
        error('certify:invalidValueEntries', ...
            'Invalid VALUE entries.') ;
        end

        elseif(iscell(mesh.point.coord))
%----------------------------------------- for GRID obj kind

        if (length(mesh.point.coord) ~= ...
            ndims (mesh.value))
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end

        if (length(mesh.point.coord) == 2)

        if (isvector(mesh.value))
        if (length(mesh.point.coord{2}) ...
           *length(mesh.point.coord{1}) ...
            ~= numel(mesh.value))
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end
        else
        if (length(mesh.point.coord{2}) ...
            ~= size(mesh.value,1) || ...
            length(mesh.point.coord{1}) ...
            ~= size(mesh.value,2) )
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end
        end

        end

        if (length(mesh.point.coord) == 3)

        if (isvector(mesh.value))
        if (length(mesh.point.coord{2}) ...
           *length(mesh.point.coord{1}) ...
           *length(mesh.point.coord{3}) ...
            ~= numel(mesh.value))
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end
        else
        if (length(mesh.point.coord{2}) ...
            ~= size(mesh.value,1) || ...
            length(mesh.point.coord{1}) ...
            ~= size(mesh.value,2) || ...
            length(mesh.point.coord{3}) ...
            ~= size(mesh.value,3) )
        error('certify:incorrectDimensions', ...
            'Invalid VALUE dimensions.') ;
        end
        end

        end

        if (any(isinf(mesh.value)))
        error('certify:invalidValueEntries', ...
            'Invalid VALUE entries.') ;
        end
        if (any(isnan(mesh.value)))
        error('certify:invalidValueEntries', ...
            'Invalid VALUE entries.') ;
        end

        else
%----------------------------------------- wrong VALUE class
        error('certify:incorrectInputClass', ...
            'Invalid VALUE class.') ;

        end
        end
    end

    if (inspect(mesh,'slope'))
%----------------------------------------- check SLOPE value
        if (~isempty  (mesh.slope))
        if (isnumeric(mesh.point.coord))
%----------------------------------------- for MESH obj kind

        if (ndims(mesh.slope) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end
        if ( size(mesh.slope,1) ~= np)
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end

        if (any(isinf(mesh.slope)))
        error('certify:invalidValueEntries', ...
            'Invalid SLOPE entries.') ;
        end
        if (any(isnan(mesh.slope)))
        error('certify:invalidValueEntries', ...
            'Invalid SLOPE entries.') ;
        end

        elseif(iscell(mesh.point.coord))
%----------------------------------------- for GRID obj kind

        if (length(mesh.point.coord) ~= ...
            ndims (mesh.slope))
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end

        if (length(mesh.point.coord) == 2)

        if (isvector(mesh.slope))
        if (length(mesh.point.coord{2}) ...
           *length(mesh.point.coord{1}) ...
            ~= numel(mesh.slope))
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end
        else
        if (length(mesh.point.coord{2}) ...
            ~= size(mesh.slope,1) || ...
            length(mesh.point.coord{1}) ...
            ~= size(mesh.slope,2) )
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end
        end

        end

        if (length(mesh.point.coord) == 3)

        if (isvector(mesh.slope))
        if (length(mesh.point.coord{2}) ...
           *length(mesh.point.coord{1}) ...
           *length(mesh.point.coord{3}) ...
            ~= numel(mesh.slope))
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end
        else
        if (length(mesh.point.coord{2}) ...
            ~= size(mesh.slope,1) || ...
            length(mesh.point.coord{1}) ...
            ~= size(mesh.slope,2) || ...
            length(mesh.point.coord{3}) ...
            ~= size(mesh.slope,3) )
        error('certify:incorrectDimensions', ...
            'Invalid SLOPE dimensions.') ;
        end
        end

        end

        if (any(isinf(mesh.slope)))
        error('certify:invalidValueEntries', ...
            'Invalid SLOPE entries.') ;
        end
        if (any(isnan(mesh.slope)))
        error('certify:invalidValueEntries', ...
            'Invalid SLOPE entries.') ;
        end

        else
%----------------------------------------- wrong VALUE class
        error('certify:incorrectInputClass', ...
            'Invalid SLOPE class.') ;

        end
        end
    end

    if (inspect(mesh,'edge2'))
%----------------------------------------- check EDGE2 index
        if (~isempty  (mesh.edge2.index))
        if (~isnumeric(mesh.edge2.index))
        error('certify:incorrectInputClass', ...
            'Invalid EDGE2.INDEX type.') ;
        end
        if (ndims(mesh.edge2.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid EDGE2.INDEX dimensions.') ;
        end
        if (size(mesh.edge2.index,2)~= 3)
        error('certify:incorrectDimensions', ...
            'Invalid EDGE2.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.edge2.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid EDGE2.INDEX indexing.') ;
        end
        if (any(isnan(mesh.edge2.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid EDGE2.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.edge2.index(:,1:2))) < +1 || ...
            max(max( ...
        mesh.edge2.index(:,1:2))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid EDGE2.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'tria3'))
%----------------------------------------- check TRIA3 index
        if (~isempty  (mesh.tria3.index))
        if (~isnumeric(mesh.tria3.index))
        error('certify:incorrectInputClass', ...
            'Invalid TRIA3.INDEX type.') ;
        end
        if (ndims(mesh.tria3.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid TRIA3.INDEX dimensions.') ;
        end
        if (size(mesh.tria3.index,2)~= 4)
        error('certify:incorrectDimensions', ...
            'Invalid TRIA3.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.tria3.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA3.INDEX indexing.') ;
        end
        if (any(isnan(mesh.tria3.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA3.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.tria3.index(:,1:3))) < +1 || ...
            max(max( ...
        mesh.tria3.index(:,1:3))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA3.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'quad4'))
%----------------------------------------- check QUAD4 index
        if (~isempty  (mesh.quad4.index))
        if (~isnumeric(mesh.quad4.index))
        error('certify:incorrectInputClass', ...
            'Invalid QUAD4.INDEX type.') ;
        end
        if (ndims(mesh.quad4.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid QUAD4.INDEX dimensions.') ;
        end
        if (size(mesh.quad4.index,2)~= 5)
        error('certify:incorrectDimensions', ...
            'Invalid QUAD4.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.quad4.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid QUAD4.INDEX indexing.') ;
        end
        if (any(isnan(mesh.quad4.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid QUAD4.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.quad4.index(:,1:4))) < +1 || ...
            max(max( ...
        mesh.quad4.index(:,1:4))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid QUAD4.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'tria4'))
%----------------------------------------- check TRIA4 index
        if (~isempty  (mesh.tria4.index))
        if (~isnumeric(mesh.tria4.index))
        error('certify:incorrectInputClass', ...
            'Invalid TRIA4.INDEX type.') ;
        end
        if (ndims(mesh.tria4.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid TRIA4.INDEX dimensions.') ;
        end
        if (size(mesh.tria4.index,2)~= 5)
        error('certify:incorrectDimensions', ...
            'Invalid TRIA4.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.tria4.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA4.INDEX indexing.') ;
        end
        if (any(isnan(mesh.tria4.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA4.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.tria4.index(:,1:4))) < +1 || ...
            max(max( ...
        mesh.tria4.index(:,1:4))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid TRIA4.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'hexa8'))
%----------------------------------------- check HEXA8 index
        if (~isempty  (mesh.hexa8.index))
        if (~isnumeric(mesh.hexa8.index))
        error('certify:incorrectInputClass', ...
            'Invalid HEXA8.INDEX type.') ;
        end
        if (ndims(mesh.hexa8.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid HEXA8.INDEX dimensions.') ;
        end
        if (size(mesh.hexa8.index,2)~= 9)
        error('certify:incorrectDimensions', ...
            'Invalid HEXA8.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.hexa8.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid HEXA8.INDEX indexing.') ;
        end
        if (any(isnan(mesh.hexa8.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid HEXA8.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.hexa8.index(:,1:8))) < +1 || ...
            max(max( ...
        mesh.hexa8.index(:,1:8))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid HEXA8.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'wedg6'))
%----------------------------------------- check WEDG6 index
        if (~isempty  (mesh.wedg6.index))
        if (~isnumeric(mesh.wedg6.index))
        error('certify:incorrectInputClass', ...
            'Invalid WEDG6.INDEX type.') ;
        end
        if (ndims(mesh.wedg6.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid WEDG6.INDEX dimensions.') ;
        end
        if (size(mesh.wedg6.index,2)~= 7)
        error('certify:incorrectDimensions', ...
            'Invalid WEDG6.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.wedg6.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid WEDG6.INDEX indexing.') ;
        end
        if (any(isnan(mesh.wedg6.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid WEDG6.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.wedg6.index(:,1:6))) < +1 || ...
            max(max( ...
        mesh.wedg6.index(:,1:6))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid WEDG6.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'pyra5'))
%----------------------------------------- check PYRA5 index
        if (~isempty  (mesh.pyra5.index))
        if (~isnumeric(mesh.pyra5.index))
        error('certify:incorrectInputClass', ...
            'Invalid PYRA5.INDEX type.') ;
        end
        if (ndims(mesh.pyra5.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid PYRA5.INDEX dimensions.') ;
        end
        if (size(mesh.pyra5.index,2)~= 6)
        error('certify:incorrectDimensions', ...
            'Invalid PYRA5.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.pyra5.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid PYRA5.INDEX indexing.') ;
        end
        if (any(isnan(mesh.pyra5.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid PYRA5.INDEX indexing.') ;
        end
        if (min(min( ...
        mesh.pyra5.index(:,1:5))) < +1 || ...
            max(max( ...
        mesh.pyra5.index(:,1:5))) > np  )
        error('certify:invalidMeshIndexing', ...
            'Invalid PYRA5.INDEX indexing.') ;
        end
        end
    end

    if (inspect(mesh,'bound'))
%----------------------------------------- check BOUND index
        if (~isempty  (mesh.bound.index))
        if (~isnumeric(mesh.bound.index))
        error('certify:incorrectInputClass', ...
            'Invalid BOUND.INDEX type.') ;
        end
        if (ndims(mesh.bound.index) ~= 2)
        error('certify:incorrectDimensions', ...
            'Invalid BOUND.INDEX dimensions.') ;
        end
        if (size(mesh.bound.index,2)~= 3)
        error('certify:incorrectDimensions', ...
            'Invalid BOUND.INDEX dimensions.') ;
        end
        if (any(isinf(mesh.bound.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid BOUND.INDEX indexing.') ;
        end
        if (any(isnan(mesh.bound.index)))
        error('certify:invalidMeshIndexing', ...
            'Invalid BOUND.INDEX indexing.') ;
        end
        if (min(min( ...
            mesh.bound.index(:,2:2))) < +1 )
        error('certify:invalidMeshIndexing', ...
            'Invalid BOUND.INDEX indexing.') ;
        end
        end
    end

%----------------------------------------- ok if we get here
    flag = +1 ;

end

author: dengwirda


text:
function h = p30_h ( n, x )

%*****************************************************************************80
%
%% P30_H evaluates the Hessian for problem 30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  a = 1.0;
  d = 6.0;
  e = 10.0;

  b = 5.1 / ( 4.0 * pi^2 );
  c = 5.0 / pi;
  ff = 1.0 / ( 8.0 * pi );

  h(1,1) = 2.0 * a * ( - 2.0 * b * x(1) + c ) ...
    * ( - 2.0 * b * x(1) + c ) ...
    - 4.0 * a * b * ( x(2) - b * x(1)^2 + c * x(1) - d ) ...
    - e * ( 1.0 - ff ) * cos ( x(1) );

  h(1,2) = 2.0 * a * ( - 2.0 * b * x(1) + c );

  h(2,1) = 2.0 * a * ( - 2.0 * b * x(1) + c );
  h(2,2) = 2.0 * a;

  return
end


text:
function [DataTest DataTrain CTest CTrain Loc_test] = samplesdivide(indian_pines_corrected,indian_pines_gt,train,randpp);

CTrain = [];
CTest = [];
DataTest  = [];
DataTrain = [];

[m n p] = size(indian_pines_corrected);
indian_pines_map = uint8(zeros(m,n));
data_col = reshape(indian_pines_corrected,m*n,p);
[mm nn] = ind2sub([m n],1:m*n);
data_col = [mm' nn' data_col];

for i = 1:max(indian_pines_gt(:))
    ci = length(find(indian_pines_gt==i));
    [v]=find(indian_pines_gt==i);
    datai = data_col(find(indian_pines_gt==i),:);
    if train>1
        cTrain = round(train);
    else
        cTrain  = round(train*ci);
    end
    cTest  = ci-cTrain;
    CTrain = [CTrain cTrain];
    CTest = [CTest cTest];
    index = randpp{i};
    DataTest = [DataTest; datai(index(1:cTest),:)];
    DataTrain = [DataTrain; datai(index(cTest+1:cTest+cTrain),:)];
end

Normalize = max(max(DataTrain(:,3:end)));
DataTrain(:,3:end) = DataTrain(:,3:end)./Normalize;
DataTest(:,3:end)  = DataTest(:,3:end)./Normalize;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function Q = cumsummat(N)
%CUMSUMMAT   Trigonometric Fourier integration matrix.
%   Q = CUMSUMMAT(N) is the matrix that maps function values at N equi-spaced
%   points to values of the integral of the interpolating trigonometric
%   polynomial at those points.

% [TODO]: Add support.
error('CHEBFUN:TRIGCOLLOC:cumsummat:notSupported', ...
    ['Indefinite integration is currently not supported for ' ...
    'TRIGCOLLOC discretization.\nPlease consider using ' ....
    'CHEBCOLLOC2 or ULTRAS discretization.']);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [uOutput] = scecdcimp(nFunction, sFilename)

% Filter function switchyard
if nFunction == FilterOp.getDescription
    uOutput = 'USGS PDE Data Center Compressed Format (stringconvert)';
elseif nFunction == FilterOp.importCatalog
    % Read formated data
    mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
    % Create empty catalog
    uOutput = zeros(length(mData), 9);
    % Loop thru all lines of catalog and convert them
    for i = 1:length(mData)
        if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
        try

            uOutput(i,1) = str2num(mData{i}(35:42));  % Longitude
            uOutput(i,2) = str2num(mData{i}(28:34));  % Latitude
            uOutput(i,3) = str2num(mData{i}(8:11));    % Year
            uOutput(i,4) = str2num(mData{i}(13:14));    % Month
            uOutput(i,5) = str2num(mData{i}(15:16));   % Day
            uOutput(i,6) = str2num(mData{i}(66:68));  % Magnitude
            uOutput(i,7) = str2num(mData{i}(43:47));  % Depth
            uOutput(i,8) = str2num(mData{i}(17:18));  % Hour
            uOutput(i,9) = str2num(mData{i}(19:20));  % Minute

             %Create decimal year
            % uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9)]);
         catch
            msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
            uOutput(i,:)=nan;
        end
    end
    l = isnan(uOutput(:,1));
    uOutput(l,:) = [];
end

author: CelsoReyes


text:
function p = mvnormpdf(varargin)
%MVNORMPDF    Multivariate normal probability density function.
% MVNORMPDF(x) returns a row vector giving the density at each column of x
%   under a standard multivariate normal.
% MVNORMPDF(x,m) subtracts m from x first.
%   If cols(m) == 1, subtracts m from each column of x.
%   If cols(m) == cols(x), subtracts corresponding columns.
%   If cols(x) == 1, x is repeated to match cols(m).
% MVNORMPDF(x,m,S) specifies the standard deviation, or more generally
%   an upper triangular Cholesky factor of the covariance matrix.
%   In the univariate case, multiple standard deviations can be specified.
%   If m is empty, no subtraction is done (zero mean).
% MVNORMPDF(x,m,[],V) specifies the variance or covariance matrix.
% MVNORMPDF(x,m,'inv',iV) specifies the inverse of the covariance matrix, i.e.
%   the precision matrix.
% MVNORMPDF(x,m,iS,'inv') specifies the reciprocal of the standard deviation,
%   or more generally the upper triangular Cholesky factor of the
%   inverse covariance matrix.
%   This is the most efficient option.
% See test_normpdf for a timing test.

% this may look strange, but computing normpdf directly is no faster or
% more stable than exp(mvnormpdfln).
p = exp(mvnormpdfln(varargin{:}));

author: andrewssobral

repo: lrslibrary


text:
% classify function
% TO RUN : 1) MatConvNet must be installed.
%          2) Modify variable 'setup_path'
% Input : 101 x 101 x 3 feature image
% Output : scores [0, 1] of each class (1~8)

function scores = classifyFeature(img, images, net)
    % 1) Input Normalization : change to 50x100x3 and subtract data_mean
    nor_im = single(img(52:101,1:100,:))- images.data_mean;

    % 2) Classify via vl_simplenn
    net.layers{end}.type = 'softmax';
    res = vl_simplenn(net, nor_im) ;

    % show the classification result
    scores = squeeze(gather(res(end).x)) ;

end


text:
function t = subsasgn(t,s,rhs)
%SUBSASGN Subscripted assignment for sparse tensor.
%
%   We can assign elements to a sptensor in three ways.
%
%   Case 1: X(R1,R2,...,RN) = Y, in which case we replace the
%   rectangular subtensor (or single element) specified by the ranges
%   R1,...,RN with Y. The right-hand-side can be a scalar or an
%   sptensor.
%
%   Case 2: X(S) = V, where S is a p x n array of subscripts and V is
%   a scalar value or a vector containing p values.
%
%   Linear indexing is not supported for sparse tensors.
%
%   Examples
%   X = sptensor([30 40 20]) %<-- Create an emtpy 30 x 40 x 20 sptensor
%   X(30,40,20) = 7 %<-- Assign a single element to be 7
%   X([1,1,1;2,2,2]) = 1 %<-- Assign a list of elements to the same value
%   X(11:20,11:20,11:20) = sptenrand([10,10,10],10) %<-- subtensor!
%   X(31,41,21) = 7 %<-- grows the size of the tensor
%   X(111:120,111:120,111:120) = sptenrand([10,10,10],10) %<-- grows
%   X(1,1,1,1) = 4 %<-- increases the number of dimensions from 3 to 4
%
%   X = sptensor([30]) %<-- empty one-dimensional tensor
%   X([4:6]) = 1 %<-- set subtensor to ones (does not increase dimension)
%   X([10;12;14]) = (4:6)'  %<-- set three elements
%   X(31) = 7 %<-- grow the first dimension
%   X(1,1) = 0 %<-- add a dimension, but no nonzeros
%
%   Note regarding singleton dimensions: It is not possible to do, for
%   instance, X(1,1:10,1:10) = sptenrand([1 10 10],5). However, it is okay
%   to do X(1,1:10,1:10) = squeeze(sptenrand([1 10 10],5)).
%
%   See also SPTENSOR, TENSOR/SUBSASGN.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

switch s.type

    case '.'
        error(['Cannot change field ', s.subs, ' directly.']);

    case '()'

        % Do nothing if both subscripts and RHS are empty
        if isempty(s.subs{1}) && isempty(t.vals)
            return;
        end

        % Figure out if we are doing a subtensor or a list of subscripts...
        type = tt_assignment_type(t,s.subs,rhs);

        %% Case I: Replace a sub-tensor
        if isequal(type,'subtensor')

            %% Case I(a): RHS is another sparse tensor
            if isa(rhs,'sptensor')

                %% First, Resize the tensor and check the size match with
                %% the tensor that's being inserted.
                m = 1;
                for n = 1:numel(s.subs)
                    if ischar(s.subs{n}) && (s.subs{n} == ':')
                        if ndims(t) < n
                            newsz(1,n) = rhs.size(m);
                        else
                            newsz(1,n) = max([t.size(n), rhs.size(m)]);
                        end
                        m = m + 1;
                    elseif numel(s.subs{n}) == 1
                        if ndims(t) < n
                            newsz(1,n) = s.subs{n};
                        else
                            newsz(1,n) = max([t.size(n) s.subs{n}]);
                        end
                    else
                        if numel(s.subs{n}) ~= rhs.size(m)
                            error('RHS does not match range size');
                        end
                        if ndims(t) < n
                            newsz(1,n) = max(s.subs{n});
                        else
                            newsz(1,n) = max([t.size(n) s.subs{n}]);
                        end
                        m = m + 1;
                    end
                end
                t.size = newsz;

                % Expand subs array if there are new modes, i.e., if the order
                % has increased.
                if ~isempty(t.subs) && (size(t.size,2) > size(t.subs,2))
                    t.subs(:,end+1:size(t.size,2)) = 1;
                end

                % Delete what currently occupies the specified range
                rmloc = subdims(s.subs,t);
                kploc = setdiff(1:nnz(t),rmloc);
                newsubs = t.subs(kploc,:);
                newvals = t.vals(kploc);

                % Renumber the subscripts
                addsubs = irenumber(rhs, t.size, s.subs);
                t.subs = [newsubs; addsubs];
                t.vals = [newvals; rhs.vals];

                return;
            end

            %% Case I(b): RHS is zero or scalar

            % First, Resize the tensor.
            % Determine new size of existing modes
            for n = 1:ndims(t)
                if ischar(s.subs{n}) && (s.subs{n} == ':')
                    newsz(1,n) = t.size(n);
                else
                    newsz(1,n) = max([t.size(n) s.subs{n}]);
                end
            end
            % Determine size of new modes, if any
            for n = ndims(t)+1:numel(s.subs)
                newsz(1,n) = max(s.subs{n});
            end
            t.size = newsz;

            % Expand subs array if there are new modes, i.e., if the order
            % has increased.
            if ~isempty(t.subs) && (size(t.size,2) > size(t.subs,2))
                t.subs(:,end+1:size(t.size,2)) = 1;
            end

            % Case I(b)i: Zero right-hand-side
            if numel(rhs) == 1 && rhs == 0

                % Delete what currently occupies the specified range
                rmloc = subdims(s.subs,t);
                kploc = setdiff(1:nnz(t),rmloc);
                t.subs = t.subs(kploc,:);
                t.vals = t.vals(kploc);
                return;
            end

            % Case I(b)ii: Scalar right-hand-side
            if numel(rhs) == 1

                % Determine number of dimensions (may be larger than
                % current number)
                N = numel(s.subs);

                % Figure out how many indices are in each dimension
                nssubs = zeros(N,1);
                for n = 1:N
                    if ischar(s.subs{n}) && s.subs{n} == ':'
                        s.subs{n} = 1:size(t,n);
                    end
                    nssubs(n) = numel(s.subs{n});
                end

                % Preallocate (discover any memory issues here!)
                addsubs = zeros(prod(nssubs),N);

                % Generate appropriately sized ones vectors.
                o = cell(N,1);
                for n = 1:N
                    o{n} = ones(nssubs(n),1);
                end

                % Generate each column of the subscripts in turn
                for n = 1:N
                    i = o;
                    i{n} = s.subs{n}';
                    addsubs(:,n) = khatrirao(i);
                end

                if ~isempty(t.subs)
                    % replace existing values
                    [junk,loc] = intersect(t.subs,addsubs,'rows');
                    t.vals(loc) = rhs;
                    % pare down list of subscripts to add
                    addsubs = setdiff(addsubs,t.subs,'rows');
                end
                t.subs = [t.subs; addsubs];
                t.vals = [t.vals; rhs*ones(size(addsubs,1),1)];
                return;
            end

            error('Invalid RHS')

        end

        % Case II: Subscripts
        if isequal(type,'subscripts')

            % Case II: Replacing values at specified indices

            newsubs = [s.subs{1}];
            tt_subscheck(newsubs);

            % Error check on subscripts
            if size(newsubs,2) < ndims(t)
                error('Invalid subscripts');
            end

            % Check for expanding the order
            if size(newsubs,2) > ndims(t)
                t.size(end+1:size(newsubs,2)) = 1;
                if ~isempty(t.subs)
                    t.subs(:,end+1:size(newsubs,2)) = 1;
                end
            end

            % Copy rhs to newvals
            newvals = rhs;

            % Error check the RHS is a column vector. We do not bother to
            % handle any other type of RHS with the sparse tensor.
            tt_valscheck(newvals);

            % Determine number of nonzeros being inserted. (This is
            % determined by the number of subscripts. Later we will check
            % to see that it matches the size of the RHS.)
            newnnz = size(newsubs,1);

            % Error check on size of newvals
            if numel(newvals) == 1

                % Special case where newvals is a single element to be
                % assigned to multiple RHS. Fix to be correct size.
                newvals = newvals * ones(newnnz,1);

            elseif size(newvals,1) ~= newnnz

                % Sizes don't match!
                error('Number of subscripts and number of values do not match!');

            end

            % Remove duplicates & print warning if any duplicates were
            % removed.
            [newsubs,idx] = unique(newsubs,'rows');
            if size(newsubs,1) ~= newnnz
                warning('Duplicate assignments discarded.');
            end
            newvals = newvals(idx);

            % Find which subscripts already exist and their locations
            [tf,loc] = ismember(newsubs,t.subs,'rows');

            % Split into three groups for processing:
            %
            % Group A: Elements that already exist and need to be changed
            % Group B: Elements that already exist and need to be removed
            % Group C: Elements that do not exist and need to be added
            %
            % Note that we are ignoring any new zero elements, because
            % those obviously do not need to be added. Also, it's
            % important to process Group A before Group B because the
            % processing of Group B may change the locations of the
            % remaining elements.

            idxa = find((tf .* newvals) ~= 0);
            idxb = find((tf .* ~abs(newvals)) ~= 0);
            idxc = find((~tf .* newvals) ~= 0);

            % Process Group A: Changing values
            if ~isempty(idxa)
                t.vals(loc(idxa)) = newvals(idxa);
            end

            % Process Group B: Removing values
            if ~isempty(idxb)
                removesubs = loc(idxb);
                keepsubs = setdiff(1:nnz(t),removesubs);
                t.subs = t.subs(keepsubs,:);
                t.vals = t.vals(keepsubs);
            end

            % Process Group C: Adding new, nonzero values
            if ~isempty(idxc)
                t.subs = [t.subs; newsubs(idxc,:)];
                t.vals = [t.vals; newvals(idxc)];
            end

            % Resize the tensor!
            for n = 1:length(t.size)
                smax = max(newsubs(:,n));
                t.size(n) = max(t.size(n), smax);
            end

            return;

        end

        error('Invalid call to sptensor/subsasgn');

    case '{}'
        error('Subscript cell reference not supported for sptensor.');

    otherwise
        error('Incorrect indexing into sptensor.')

end

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function cpp = subsetconnectivity ( t, p, nds )

%*****************************************************************************80
%
%% SUBSETCONNECTIVITY: connectivity of boundary node subset in a triangulation.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    29 April 2014
%
%  Author:
%
%    Marcus Garvie
%
%  Parameters:
%
%    Input, real P(NP,2), the coordinates of a set of nodes.
%
%    Input, integer T(NT,3), a list of the nodes which make up each
%    triangle of a triangulation of the nodes in P.
%
%    Input, integer NDS(NN,1), a list of nodes comprising a subset of the
%    boundary nodes.
%
%    Output, CPP(NE,2), a list of edges, each edge defined by 2 nodes from
%    the list NDS.
%
  cpp = [];
%
% Work out connectivity for the whole boundary.
%
  edges = boundedges (p,t);
%
%  NO_EDGES = Number of edges in whole boundary.
%
  [no_edges,~] = size(edges);
%
%  NN = Number of nodes in the subset of boundary nodes.
%
  NN = length(nds);
%
%  Find the edges in the subset of the boundary NDS.
%
  count = 0;

  for i = 1:no_edges

    node1 = edges(i,1);
    node2 = edges(i,2);
    bin1 = ismember(node1,nds);
    bin2 = ismember(node2,nds);

    if (bin1==1) & (bin2==1)
      count = count + 1;
      cpp(count,:) = edges(i,:);
    end

  end

  return
end


text:
function [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
%LTFATDIDITFAIL  Did a test fail
%
%  [test_fail,fail]=LTFATDIDITFAIL(res,test_fail) updates test_fail if
%  res is above threshhold and outputs the word FAIL in the variable
%  fail. Use only in testing scripts.
global LTFAT_TEST_TYPE;
if nargin<3
  tolerance=1e-10;
  if strcmpi(LTFAT_TEST_TYPE,'single')
     tolerance=2e-4;
  end
end;

fail='';
if (abs(res)>tolerance) || isnan(res)
  fail='FAILED';
  test_failed=test_failed+1;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function applyAllPVC_HN(pathWORK,patientNames)
% -------------------------------------------------------------------------
% function applyAllPVC_HN(pathWORK,patientNames)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function applies partial-volume corrections for all PET volume of
% the head and neck cohort according to the methodology developed in ref.
% [1].
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Boussion, N. et al. (2009). Incorporation of wavelet-based denoising
%     in iterative deconvolution for partial volume correction in
%     whole-body PET imaging. Eur J Nucl Med Mol Imaging, 36(7), 1064-1075.
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the HN WORKSPACE directory.
% - patientNames: Cell of strings specifying the patient names to read.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd([pathWORK,'/DATA'])
nPatient = length(patientNames);

% PVE CORRECTIONS
for i = 1:nPatient
    fprintf('\n***** CORRECTING %s *****\n',patientNames{i})
    load(patientNames{i}) % Variable 'sData' now in MATLAB Workspace
    volumePVC = PVEcorrect(sData{2}.scan.volume.data,2); % 2 iterations have proved to be sufficient
    sData{2}.scan.volume.data = volumePVC;
    sData{2}.scan.info = 'PVE corrected';
    save(patientNames{i},'sData')
end

cd(startpath)
end


text:
function demo_tcca()

% Set the various paths
rootpath = filesepStandard(fileparts(which('Homer3.m')));
currpath = pwd;
subjFolder = 'UnitTests/Example6_GrpTap';
subjFolder = [rootpath, subjFolder];
procStreamConfigFile = './processOpt_tcca.cfg';   % proc stream config file path relative to subject folder
cd(subjFolder);
resetGroupFolder('');

fprintf('Loading subject folder:  %s\n', subjFolder);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change processing stream param values to newval
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dataTree = LoadDataTree(subjFolder, [], procStreamConfigFile);
iG = dataTree.GetCurrElemIndexID();
dataTree.groups(iG).Calc();
dataTree.groups(iG).Save();

cd(currpath)

author: BUNPC


text:
function DCM = spm_dcm_erp_update(DCM,oldDCM,fields)
% Set priors over DCM model parameters for Bayesian updating
% FORMAT DCM = spm_dcm_erp_update(DCM,oldDCM,fields)
%
% DCM    - DCM structure to be updated
% oldDCM - inverted DCM with posterior moments
% fields - character array of fields to be updated: e.g.,{'A','B'}
%__________________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_erp_update.m 5801 2013-12-10 18:45:20Z guillaume $

% get prior structure for the sort of model
%==========================================================================
try, model = DCM.options.model; catch, model = 'CMC'; end

% get the posterior (expectation and covariance) moments of parameters
%--------------------------------------------------------------------------
Ep      = oldDCM.Ep;
Cp      = spm_unvec(diag(oldDCM.Cp),Ep);

% prior moments of parameters
%--------------------------------------------------------------------------
[pE,pC] = spm_dcm_neural_priors(DCM.A,DCM.B,DCM.C,model);

% fill-in prior moments
%--------------------------------------------------------------------------
for i = 1:length(fields)

    pE.(fields{i}) = Ep.(fields{i});
    pC.(fields{i}) = Cp.(fields{i});

end

% they surprising model structure
%==========================================================================
DCM.M.pE = pE;
DCM.M.pC = pC;


text:
%  M file that starts a Mosaic Browser displaying
%  the HTML version of the ZMAP Users Guide
%
%  Stefan Wiemer   6/96

%    !netscape /Seis/A/stefan/zmapwww/title.htm  &

report_this_filefun(mfilename('fullpath'));
disp('Attempting to open browser - please be patient...');
disp('If the browser does not open, please open the browser first and try again, or open the file ./zmapwww/title.hmt manually');
do = [ 'web ' hodi '/zmapwww/title.htm ;' ];
err=['errordlg('' Error while opening, please open the browser first and try again or open the file ./zmapwww/title.hmt manually'');'];
eval(do,err)


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

function test_k_nearest(testCase)
    m = [1.0000    0.6000    0.2000
    0.6000    1.0000    0.3000
    0.2000    0.3000    1.0000];
    m = spx.cluster.similarity.filter_k_nearest_neighbors(m, 2);
    verifyEqual(testCase, m, [ 1.0000    0.6000         0
    0.6000    1.0000    0.3000
         0    0.3000    1.0000]);
    % We consider a graph with [[2,3,4], [1], [1, 4], [1, 3]]
    m = [
    0 1 1 1
    1 0 0 0
    1 0 0 1
    1 0 1 0];
    m2 = spx.cluster.similarity.filter_k_nearest_neighbors(m, 2);
    verifyEqual(testCase, m2, m);
    m2 = spx.cluster.similarity.filter_k_nearest_neighbors(m, 1);
    verifyEqual(testCase, m2, [0 1 1 1
        1 0 0 0
        1 0 0 0
        1 0 0 0]);
end

function test_sim_1(testCase)
    x = [1 2 3];
end

author: indigits

repo: sparse-plex


text:
 function [opt, extra] = vararg_pair(opt, varargs, varargin)
%function [opt, extra] = vararg_pair(opt, varargs, [options])
%|
%| Process name / value pairs, replacing the "default" field values
%| of the opt structure with the user-specified values.
%| This allows flexible argument order and "named arguments" somewhat like IDL.
%| This newer version allows option names to differ from the field names.
%| If two output arguments are given, then any "extra" name value pairs
%| that are not associated with the fields in structure opt will be
%| returned as name/value pairs in the "extra" output.
%|
%| in
%|	opt	struct	opt.name1 opt.name2 etc., containing default values
%|	varargs	{'name1', 'value1', ... } as a cell array
%|
%| options
%|	'subs'	Nx2 cell	{'arg_name1', 'struct_name1'; ...}
%|		short argument name followed by longer structure name
%|	'allow_new'	0|1	allow new structure elements? (default: 0)
%|
%| out
%|	opt		opt.(name) = value
%|	extra		{'name_i', 'value_i', ...} all "unused" pairs
%|
%| Copyright 2005-6-18, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(opt, 'test'), vararg_pair_test, clear, return, end

if isempty(varargs) && numel(varargin) == 0
	extra = {};
	return
end

base = [caller_name ' arg: ']; % for printing

local.subs = {};
local.allow_new = 0;
local = vararg_pair_do(local, varargin, {}, 0, 0);

allow_extra = nargout > 1;
[opt extra] = vararg_pair_do(opt, varargs, local.subs, local.allow_new, ...
	allow_extra, base);

% Structure elements of opt with value "NaN" are checked as being required.
% NO, this just adds overhead.  Removed 2007-2-9.
% vararg_pair_check(opt);

%
% vararg_pair_do()
%
function [opt, extra] = ...
	vararg_pair_do(opt, varargs, subs, allow_new, allow_extra, base)

if allow_extra && allow_new
	error 'only one of "new" or "extra" is sensible for vararg_pair'
end

npair = floor(length(varargs) / 2);
if 2*npair ~= length(varargs), fail('need names and values in pairs'), end
args = varargs(1:2:end);
vals = varargs(2:2:end);

if ~isempty(subs)
	if size(subs,2) ~= 2, fail('subs must be Nx2'), end
%	subs1 = strvcat(subs(:,1));
	subs1 = char(subs(:,1)); % strmatch doc says fastest for char array!
	subs2 = subs(:,2);
	% todo: apply 'subs' all at once here?
end

extra = {};
for ii=1:npair
	arg = args{ii};
	val = vals{ii};

	if ~ischar(arg)
		fail('unknown option of class %s', class(arg))
	end

	if ~isempty(subs)
		arg = vararg_pair_sub(arg, subs1, subs2);
	end

	[opt ex] = vararg_pair_put(opt, arg, val, allow_new, allow_extra);
	if allow_extra && ~isempty(ex)
		extra = [extra(:); ex(:)]'; % faster than {extra{:}, ex{:}};
	end

	if isfield(opt, 'chat') && opt.chat && isempty(ex)
		show_pair(arg, val, base)
	end
end

%
% vararg_pair_sub()
% substitute option name for field name, if a Nx2 cell array is provided.
%
function arg = vararg_pair_sub(arg, subs1, subs2)
ii = strmatch(arg, subs1, 'exact');
if isempty(ii)
	return
elseif length(ii) == 1
	arg = subs2{ii};
else
	printm([subs1(ii) ' ' subs2{ii}])
	error 'double match?  only one substitution allowed'
end

%
% vararg_pair_put()
% opt.(arg) = value
% with extra effort to handle opt.(arg1.arg2) = value
% because matlab does not handle that case by itself, except via subsagn
%
function [opt, extra] = vararg_pair_put(opt, arg, value, allow_new, allow_extra)

extra = {};
if isfield(opt, arg) % simple case of name / value pair
	opt.(arg) = value;
	return
end

idot = strfind(arg, '.');
if isempty(idot)
	if allow_new
		opt.(arg) = value;
		return
	elseif allow_extra
		extra = {arg, value};
		return
	else
		fail('unknown option name "%s"', arg)
	end
end

% tricky case of a.b.c
if length(idot) > 1, error 'a.b.c.d not done', end

arg1 = arg([1:(idot-1)]);
arg2 = arg([(idot+1):end]);
if ~isfield(opt, arg1)
	fail('unknown option1 name "%s"', arg1)
end
if ~isfield(opt.(arg1), arg2)
	fail('unknown option2 name %s', arg2)
end
s = struct('type', {'.', '.'}, 'subs', ...
	{ arg([1:(idot-1)]), arg([(idot+1):end]) });
try
	opt = subsasgn(opt, s, value);
catch
	printf(lasterr)
	fail('subsasgn? unknown option name %s', arg)
end

%
% vararg_pair_check()
% check that the user supplied the required arguments,
% which are signified by 'nan' values.
%
function vararg_pair_check(opt)
names = fieldnames(opt);
ok = 1;
for ii=1:length(names)
	t = opt.(names{ii});
	if isnumeric(t) && any(isnan(t(:)))
		printf('Required option "%s" missing', names{ii});
		ok = 0;
	end
end
if ~ok, error 'missing required "option(s)"', end

%
% show_pair()
%
function show_pair(name, value, base)

pri = @(varargin) printf([base varargin{1}], varargin{2:end});

if isnumeric(value) || islogical(value)
	if max(size(value)) == 1	% only print scalars
		if ~streq(name, 'chat')
			pri('%s = %g', name, value)
		end
	else	% otherwise print sum (e.g., for mask)
		pri('sum(%s) = %g', name, sum(value(:)))
	end

elseif ischar(value)
	pri('%s = %s', name, value)

elseif isa(value, 'cell')
	pri('%s {cell}', name)

elseif isa(value, 'function_handle')
	pri('%s = @%s', name, func2str(value))

elseif isstruct(value)	% show first structure element, if string
	fields = fieldnames(value);
%	name = getfield(arg, fields{1});
	name1 = value.(fields{1});
	if ischar(name1)
		pri('%s.%s = %s', name, fields{1}, name1)
	end

else
	warn('display type "%s" not implemented', class(value))
end

%
% test routine
%
function vararg_pair_test
args = {'a', 1, 'b', 2, 'c.d', 3, 'g', 9, 'l', 5, 'll', 6};%, 'c.f', 4};
opt.a = 0;
opt.b = 0;
opt.f = 0;
opt.long1 = 0;
opt.long2 = 0;
%opt.g = nan;
opt.c.d = 0;
opt.c.e = 0;
sub = {'g', 'f'; 'l', 'long1'; 'll', 'long2'};
opt = vararg_pair(opt, args, 'subs', sub);
jf_equal(opt.a, 1)
jf_equal(opt.b, 2)
jf_equal(opt.c.d, 3)
jf_equal(opt.c.e, 0)
jf_equal(opt.f, 9)
jf_equal(opt.long1, 5)
jf_equal(opt.long2, 6)

opt = vararg_pair(opt, {'z', 'new'}, 'allow_new', 1);
if ~isequal(opt.z, 'new'), fail('bug'), end
try % typo in name should yield error
	opt = vararg_pair(opt, {'verify_bad_test!', 'new'}, 'allow_new', 0)
catch
	return
end
error 'should not get here'

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [Q,Z,decomp,Den] = findsos(P,flag,options)
% FINDMATRIXSOS --- Find a sum of squares decomposition of a given matrix polynomial.
%
% [Q,Z,decomp,Den] = findsos(P,flag,options)
%
% P is a symmetric polynomial matrix.
%
% FLAG is an optional argument which, if set to 'rational', returns a
% sum of squares decomposition with rational coefficients.
%
% OPTIONS is an optional argument which specifies the solver and the
% solver-specific parameters on its fields as
%   options.solver
%   options.params
% the default value for options.solver is 'SeDuMi'. The solver parameters
% are fields of options.params as, for instance, options.params.tol = 1e-9.
%
% A positive semidefinite Q and a symbolic monomial vector Z will be
% computed such that
%
%    (Ir kron Z)' * Q * (Ir kron Z) = P(x)
%
% If P is not a sum of squares, the function will return empty Q and Z.
%
% If P is a polynomial with integer coefficients and is represented as a
% symbolic object, then [Q,Z] = findsos(P,'rational') will compute a rational
% matrix Q such that Z'*Q*Z = P.

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%

% 12/27/01 - SP
% 03/27/02 - SP

if nargin == 2
    if ~strcmp(flag,'rational')
        options=flag;
        flag = 'abcdefgh';
    else
        %options.solver='sedumi';
        options.solver='sdpt3';
        flag = 'abcdefgh';
    end
elseif nargin==1
    %options.solver='sedumi';
    options.solver='sdpt3';
end

if isa(P,'sym')

    dimp = size(P);
    if dimp(1)~=dimp(2) ;
        disp(['The polynomial matrix is not square, it cannot be a sum of squares']);
        Q = [];
        Z = [];
        return;
    end;

    P = expand(P);
    vars = findsym(P);
    vars = sym(['[',vars,']']);

    nvars = size(vars,2) ;

    for var = 1:nvars;
        for j = 1:dimp(1)
            if rem(double(feval(symengine,'degree',P(j,j))),2) ;
                disp(['Degree in ' char(vars(var)) ' is not even for the diagonal elements. The polynomial matrix cannot be a sum of squares']);
                Q = [];
                Z = [];
                return;
            end;
        end;
    end;

    if isequal(P,P.')~=1;
        disp(['The polynomial matrix is not symmetric, it can not be a sum of squares']);
        Q = [];
        Z = [];
        return;
    end;

    prog = sosprogram(vars);
    prog = sosineq(prog,P);
    [prog,info] = sossolve(prog,options);

    if strcmp(options.solver,'SDPNAL')
        disp('findsos function currently not supported for SDPNAL solver');
        Q = [];
        Z = [];
        if nargout >= 3
            decomp = [];
            Den = [];
        end
        return;
    elseif (info.dinf==1)||(info.pinf==1)
        disp('No sum of squares decomposition is found.');
        Q = [];
        Z = [];
        if nargout >= 3
            decomp = [];
            Den = [];
        end
        return;
    else
        Q = reshape(prog.solinfo.RRx,sqrt(length(prog.solinfo.RRx)),sqrt(length(prog.solinfo.RRx)));
        if isa(P,'sym');
            Z = mysympower(vars,prog.extravar.Z{1});
        else
            pvar Z
            coefftemp = speye(size(prog.extravar.Z{1},1));
            Z = polynomial(coefftemp,prog.extravar.Z{1},prog.vartable,[size(prog.extravar.Z{1},1),1]);
        end;
    end;

    if nargin > 1 & strcmp(flag,'rational')

        A = (full(prog.expr.At{1})') ;
        B = (full(prog.expr.b{1})) ;
        % Round x0
        N = 1 ; % This is the tolerance in the rounding
        xx = prog.solinfo.x;

        kmax = 10 ;
        pd = 0 ;
        while kmax ~= 0 ;
            kmax = kmax - 1 ;
            x0 = round(N*xx);
            [Qflat,NN] = proj3(x0,A,B,N);
            n = sqrt(length(Qflat));
            Qr = reshape(Qflat,n,n);
            % Qr should be PSD (should really check symbolically)
            if min(eig(Qr/NN))>-1e-14 ; kmax=0 ; pd = 1 ; end
            % Increase N, and try again
            N = 2*N ;
        end

        % Experimental, no good error checking yet, so we check that
        % expand(NN*P - Z.'*Qr*Z) is zero!

        if (expand(NN*P-Z.'*Qr*Z) ~= 0) | (pd == 0);
            Qr=[];Z=[];NN=[];
            disp('Could not compute a rational SOS!');
        end

        if nargout == 4
            Q = Qr;
            Den = NN;
        else
            if isa(P,'sym')
                Q = sym(Qr/NN);
            else
                Q = Qr/NN;
                disp('To obtain an exact rational solution, run the function with three output arguments.');
            end;
        end;

    end;

    L = sqrtm(double(Q));
    decomp   = L*(kron(eye(dimp(1)),Z));

else

    % PJS -- Handle polynomial variables
    % Most of this code simply mimics the symbolic case and hence the
    % overlap can be reduced.

    dimp = size(P);
    if dimp(1)~=dimp(2) ;
        disp(['The polynomial matrix is not square, it cannot be a sum of squares']);
        Q = [];
        Z = [];
        return;
    end;

    nvars = P.nvar;
    vars = polynomial(zeros(nvars,1));
    for j=1:nvars
        vars(j) = pvar(P.varname{j});
    end

    for var = 1:nvars;
        for j = 1:dimp(1)
            Pjj = P(j,j);
            maxdeg = max(Pjj.degmat(:,var));
            if rem(maxdeg,2)
                disp(['Degree in ' char(vars(var)) ' is not even for the diagonal elements. The polynomial matrix cannot be a sum of squares']);
                Q = [];
                Z = [];
                return;
            end;
        end;
    end;

    tmp = isequal(P,P.');
    if ~all(tmp(:))
        disp(['The polynomial matrix is not symmetric, it can not be a sum of squares']);
        Q = [];
        Z = [];
        return;
    end;

    prog = sosprogram(vars);
    prog = sosineq(prog,P);
    [prog,info] = sossolve(prog);

    if 'solver' == 'SDPNAL'
        disp('findsos function currently not supported for SDPNAL solver');
        Q = [];
        Z = [];
        if nargout >= 3
            decomp = [];
            Den = [];
        end
        return;
    elseif (info.dinf==1)||(info.pinf==1)
        disp('No sum of squares decomposition is found.');
        Q = [];
        Z = [];
        if nargout >= 3
            decomp = [];
            Den = [];
        end
        return;
    else
        Q = reshape(prog.solinfo.RRx,sqrt(length(prog.solinfo.RRx)),sqrt(length(prog.solinfo.RRx)));
        if isa(P,'sym');
            Z = mysympower(vars,prog.extravar.Z{1});
        else
            pvar Z
            coefftemp = speye(size(prog.extravar.Z{1},1));
            Z = polynomial(coefftemp,prog.extravar.Z{1},prog.vartable,[size(prog.extravar.Z{1},1),1]);
        end;
    end;

    if nargin > 1 & strcmp(flag,'rational')

        A = (full(prog.expr.At{1})') ;
        B = (full(prog.expr.b{1})) ;
        % Round x0
        N = 1 ; % This is the tolerance in the rounding
        xx = prog.solinfo.x;

        kmax = 10 ;
        pd = 0 ;
        while kmax ~= 0 ;
            kmax = kmax - 1 ;
            x0 = round(N*xx);
            [Qflat,NN] = proj3(x0,A,B,N);
            n = sqrt(length(Qflat));
            Qr = reshape(Qflat,n,n);
            % Qr should be PSD (should really check symbolically)
            if min(eig(Qr/NN))>-1e-14 ; kmax=0 ; pd = 1 ; end
            % Increase N, and try again
            N = 2*N ;
        end

        % Experimental, no good error checking yet, so we check that
        % expand(NN*P - Z.'*Qr*Z) is zero!

        if (expand(NN*P-Z.'*Qr*Z) ~= 0) | (pd == 0);
            Qr=[];Z=[];NN=[];
            disp('Could not compute a rational SOS!');
        end

        if nargout == 4
            Q = Qr;
            Den = NN;
        else
            if isa(P,'sym')
                Q = sym(Qr/NN);
            else
                Q = Qr/NN;
                disp('To obtain an exact rational solution, run the function with three output arguments.');
            end;
        end;

    end;

    L = sqrtm(double(Q));
    decomp   = L*(kron(eye(dimp(1)),Z));

    % PJS -- Comment out original code in this section
%     varname = P.var;
%     vars = [];
%     for i = 1:size(varname,1)
%         pvar(varname{i});
%         vars = [vars eval(varname{i})];
%     end;
%
%     deg = mod(max(P.degmat,[],1),2);
%     if sum(deg)~=0
%         i = find(deg~=0);
%         disp(['Degree in ' varname{i(1)} ' is not even. The polynomial cannot be a sum of squares']);
%         Q = [];
%         Z = [];
%         return;
%     end;

end;

author: yu-jiang

repo: radpbook

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


text:
%%%%%%%%%%%%%%%%%%%%%%%%% UPDATE PROCEDURES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% UPDATE THE GLOBAL STATE VECTORS FROM ENTITY.VIRTUAL
function [METAObjUpdate]	= OMAS_updateGlobalStates(SIM,objectID,globalVelocity_k,quaternion_k,idleStatus_k)
% This function reallocates the global properties of each entity to the
% META.OBJECT series to allow faster reference and increased independance
% of the main cycle from the object cycles.

% CONFIRM OBJECT VARIABLES
assert(isa(objectID,'uint16')   == 1,'The objectID must be a uint16 type');
assert(islogical(idleStatus_k)  == 1,'The idle agent status must be reported as a logical');
assert(size(globalVelocity_k,1) == 3 && isnumeric(globalVelocity_k),'Object velocity update must be given as a column vector [3x1].');
assert(size(quaternion_k,1)     == 4 && isnumeric(quaternion_k),'Object quaternion update must be given as a column vector [4x1].');
assert(any(isnan(quaternion_k)) == 0,sprintf('Quaternion for objectID %d is invalid; q = [%f %f %f %f]',int8(objectID),quaternion_k(1),quaternion_k(2),quaternion_k(3),quaternion_k(4)));

% IDENTIFY THE ASSOCIATED META OBJECT
logicalIndices = SIM.globalIDvector == objectID;                           % The logical position of the ID
index = inf;
for i = 1:numel(logicalIndices)
   if logicalIndices(i)
      index = i;
      break
   end
end

% EXTRACT THE META STRUCTURE TO BE UPDATED
METAObjUpdate = SIM.OBJECTS(1,index);                                      % Get the current META object associated with 'entity'

% //////////////////////// UPDATE META.OBJECT /////////////////////////////
% UPDATE THE META ROTATION DEFINING FIXED LOCAL >> ROTATED IN THE GLOBAL
METAObjUpdate.R = OMAS_geometry.quaternionToRotationMatrix(quaternion_k);
% UPDATE THE GLOBAL POSITION
globalPosition_k = METAObjUpdate.globalState(1:3) + globalVelocity_k*SIM.TIME.dt; % Calculate the new global position
% REBUILD GLOBAL STATES (ENTITY & META)
METAObjUpdate.globalState = [globalPosition_k;globalVelocity_k;quaternion_k];
% DETERMINE IF ENTITY HAS INDICATED THAT TASK IS COMPLETE
METAObjUpdate.idleStatus = idleStatus_k;
end

repo: OpenMAS


text:
function publish_html(filename, outputDir, stylesheet)

% publish_html - publish a file to HTML format
%
%    publish_html(filename, outputDir, stylesheet);
%
%   Copyright (c) 2008 Gabriel Peyre

if nargin<1
    filename = 'content';
end
if nargin<2
    outputDir = 'html';
end
if nargin<3
    stylesheet = [outputDir '/gpeyre.xsl'];
end

opts.outputDir = outputDir;
if not(isempty(stylesheet))
    opts.stylesheet = stylesheet;
end
file = publish(filename,opts);

web(file);

author: alecjacobson

repo: gptoolbox


text:
function [g,L,info] = gabpars_from_window(g,a,M,L,callfun)
%GABPARS_FROM_WINDOW  Compute g and L from window
%   Usage: [g,g.info,L] = gabpars_from_window(f,g,a,M);
%
%   Use this function if you know a window and a lattice
%   for the DGT. The function will calculate a transform length L and
%   evaluate the window g into numerical form.
%
%   If the transform length is unknown (as it usually is unless explicitly
%   specified by the user), set L to be [] in the input to this function.

if nargin<5
  stacknames=dbstack;
  callfun=stacknames(2).name;
end;

assert_squarelat(a,M,1,callfun,0);

if ~isempty(L)
  if (prod(size(L))~=1 || ~isnumeric(L))
    error('%s: L must be a scalar',callfun);
  end;

  if rem(L,1)~=0
    error('%s: L must be an integer',callfun);
  end;
end;

if isnumeric(g)
  Lwindow=length(g);
else
  Lwindow=0;
end;

if isempty(L)
  % Smallest length transform.
  Lsmallest=lcm(a,M);

  % Choose a transform length larger than both the length of the
  % signal and the window.
  L=ceil(Lwindow/Lsmallest)*Lsmallest;
else

  if rem(L,M)~=0
    error('%s: The length of the transform must be divisable by M = %i',...
          callfun,M);
  end;

  if rem(L,a)~=0
    error('%s: The length of the transform must be divisable by a = %i',...
          callfun,a);
  end;

  if L<Lwindow
    error('%s: Window is too long.',callfun);
  end;

end;

b=L/M;
N=L/a;

[g,info]=gabwin(g,a,M,L,[0 1],'callfun',callfun);


text:
function msm_to_st_test01 ( )

%*****************************************************************************80
%
%% MSM_TO_ST_TEST01 tests MSM_TO_ST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_ST_TEST01\n' );
  fprintf ( 1, '  We will create a simple 5 by 5 matrix,\n' );
  fprintf ( 1, '  convert it to MATLAB sparse format,\n' );
  fprintf ( 1, '  and then have MSM_TO_ST write the matrix to an ST file.\n' );
%
%  Create a simple matrix.
%
  a = [ 11,  0,  0, 14,  0;
         0, 22,  0,  0,  0;
        31, 32, 33, 34, 35;
         0,  0,  0, 44, 45;
        51, 52,  0,  0, 55 ];
%
%  Make a sparse version of the matrix.
%
  b = sparse ( a );
%
%  Have MSM_TO_ST write the sparse matrix to an ST file.
%
  msm_to_st ( b, 'msm_to_st_test01.st' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [u, v] = HS(im1, im2, alpha, ite, uInitial, vInitial, displayFlow, displayImg)
% Horn-Schunck optical flow method
% Horn, B.K.P., and Schunck, B.G., Determining Optical Flow, AI(17), No.
% 1-3, August 1981, pp. 185-203 http://dspace.mit.edu/handle/1721.1/6337
%
% Usage:
% [u, v] = HS(im1, im2, alpha, ite, uInitial, vInitial, displayFlow)
% For an example, run this file from the menu Debug->Run or press (F5)
%
% -im1,im2 : two subsequent frames or images.
% -alpha : a parameter that reflects the influence of the smoothness term.
% -ite : number of iterations.
% -uInitial, vInitial : initial values for the flow. If available, the
% flow would converge faster and hence would need less iterations ; default is zero.
% -displayFlow : 1 for display, 0 for no display ; default is 1.
% -displayImg : specify the image on which the flow would appear ( use an
% empty matrix "[]" for no image. )
%
% Author: Mohd Kharbat at Cranfield Defence and Security
% mkharbat(at)ieee(dot)org , http://mohd.kharbat.com
% Published under a Creative Commons Attribution-Non-Commercial-Share Alike
% 3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
% October 2008
% Rev: Jan 2009

%% Default parameters
if nargin<1 || nargin<2
    im1=imread('yos9.tif');
    im2=imread('yos10.tif');
end
if nargin<3
    alpha=1;
end
if nargin<4
    ite=100;
end
if nargin<5 || nargin<6
    uInitial = zeros(size(im1(:,:,1)));
    vInitial = zeros(size(im2(:,:,1)));
elseif size(uInitial,1) ==0 || size(vInitial,1)==0
    uInitial = zeros(size(im1(:,:,1)));
    vInitial = zeros(size(im2(:,:,1)));
end
if nargin<7
    displayFlow=1;
end
if nargin<8
    displayImg=im1;
end

%% Convert images to grayscale
if size(size(im1),2)==3
    im1=rgb2gray(im1);
end
if size(size(im2),2)==3
    im2=rgb2gray(im2);
end
im1=double(im1);
im2=double(im2);

im1=smoothImg(im1,1);
im2=smoothImg(im2,1);

tic;

%%
% Set initial value for the flow vectors
u = uInitial;
v = vInitial;

% Estimate spatiotemporal derivatives
[fx, fy, ft] = computeDerivatives(im1, im2);

% Averaging kernel
kernel_1=[1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12];

% Iterations
for i=1:ite
    % Compute local averages of the flow vectors
    uAvg=conv2(u,kernel_1,'same');
    vAvg=conv2(v,kernel_1,'same');
    % Compute flow vectors constrained by its local average and the optical flow constraints
    u= uAvg - ( fx .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2);
    v= vAvg - ( fy .* ( ( fx .* uAvg ) + ( fy .* vAvg ) + ft ) ) ./ ( alpha^2 + fx.^2 + fy.^2);
end

u(isnan(u))=0;
v(isnan(v))=0;

%% Plotting
if displayFlow==1
    plotFlow(u, v, displayImg, 5, 5);
end

author: Sable

repo: mcbench-benchmarks


text:
function [w,cxe,w_pen,optimizerState,converged] = ...
    train_binary_classifier(classifier,classf,w0,objective_function,prior,...
                            penalizer,lambda,maxiters,maxCG,optimizerState,...
                            quiet,cstepHessian)
%
%   Supervised training of a regularized fusion.
%
%
% Inputs:
%
%   classifier: MV2DF function handle that maps parameters to llr-scores.
%               Note: The training data is already wrapped in this handle.
%
%   classf: 1-by-N row of class labels:
%                -1 for non_target,
%                +1 for target,
%                 0 for ignore
%
%   w0: initial parameters. This is NOT optional.
%
%   objective_function: A function handle to an Mv2DF function that
%                       maps the output (llr-scores) of classifier, to
%                       the to-be-minimized objective (called cxe).
%                       optional, use [] to invoke 'cllr_obj'.
%
%  prior: a prior probability for target to set the 'operating point'
%         of the objective function.
%         optional: use [] to invoke default of 0.5
%
%  penalizer: MV2DF function handle that maps parameters to a positive
%              regularization penalty.
%
%  lambda: a weighting for the penalizer
%
%  maxiters: the maximum number of Newton Trust Region optimization
%            iterations to perform. Note, the user can make maxiters
%            small, examine the solution and then continue training:
%            -- see w0 and optimizerState.
%
%
%
%  optimizerState: In this implementation, it is the trust region radius.
%                  optional:
%                    omit or use []
%                    If not supplied when resuming iteration,
%                    this may cost some extra iterations.
%                  Resume further iteration thus:
%   [w1,...,optimizerState] = train_binary_classifier(...);
%   ... examine solution w1  ...
%   [w2,...,optimizerState] = train_binary_classifier(...,w1,...,optimizerState);
%
%
%  quiet: if false, outputs more info during training
%
%
%  Outputs:
%    w: the solution.
%    cxe: normalized multiclass cross-entropy of the solution.
%         The range is 0 (good) to 1(useless).
%
%    optimizerState: see above, can be used to resume iteration.
%

if nargin==0
    test_this();
    return;
end

if ~exist('maxCG','var') || isempty(maxCG)
    maxCG = 100;
end

if ~exist('optimizerState','var')
    optimizerState=[];
end

if ~exist('prior','var') || isempty(prior)
    prior = 0.5;
end

if ~exist('objective_function','var') || isempty(objective_function)
    objective_function = @(w,T,weights,logit_prior) cllr_obj(w,T,weights,logit_prior);
end

%prior_entropy = -prior*log(prior)-(1-prior)*log(1-prior);
prior_entropy = objective_function([0;0],[1,-1],[prior,1-prior],logit(prior));

classf = classf(:)';

ntar = length(find(classf>0));
nnon = length(find(classf<0));
N = nnon+ntar;

weights = zeros(size(classf));
weights(classf>0) = prior/(ntar*prior_entropy);
weights(classf<0) = (1-prior)/(nnon*prior_entropy);
%weights remain 0, where classf==0

w=[];

if exist('penalizer','var') && ~isempty(penalizer)
    obj1 = objective_function(classifier,classf,weights,logit(prior));
    obj2 = penalizer(w);
    obj = sum_of_functions(w,[1,lambda],obj1,obj2);
else
    obj = objective_function(classifier,classf,weights,logit(prior));
end

w0 = w0(:);

if exist('cstepHessian','var') &&~ isempty(cstepHessian)
    obj = replace_hessian([],obj,cstepHessian);
end

[w,y,optimizerState,converged] = trustregion_newton_cg(obj,w0,maxiters,maxCG,optimizerState,[],1,quiet);

if exist('penalizer','var') && ~isempty(penalizer)
    w_pen = lambda*obj2(w);
else
    w_pen = 0;
end

cxe = y-w_pen;
if ~quiet
    fprintf('cxe = %g, pen = %g\n',cxe,w_pen);
end

function test_this()

%invoke test for linear_fuser, which calls train_binary_classifier
linear_fuser();

author: nesl