text:
function  noisyOut = findNoisyChannels(signal, noisyIn)
% Identify bad channels in EEG using a two-stage approach
%
% reference = findNoisyChannels(signal)
% reference = findNoisyChannels(signal, reference)
%
% First remove bad channels by amplitude, noise level, and correlation
% Apply ransac after these channels have been removed.
%
% Input parameters:
%     signal - structure with srate, chanlocs, chaninfo, and data fields
%     noisyIn - structure with input parameters
%
%  Notes: the signal is assumed to be high-passed. Removing line noise
%  is a good idea too.
%
%  noisyIn: (fields are filled in on input if not present and propagated to output)
%     name - name of the input file
%     srate - sample rate in HZ
%     samples - number of samples in the data
%     evaluationChannels - a vector of channels to use
%     channelLocations - a structure of EEG channel locations
%     chaninfo - standard EEGLAB chaninfo (nose direction is relevant)
%     chanlocs - standard EEGLAB chanlocs structure
%     robustDeviationThreshold - z score cutoff for robust channel deviation
%     highFrequencyNoiseThreshold -  z score cutoff for SNR (signal above 50 Hz)
%     correlationWindowSeconds - correlation window size in seconds (default = 1 sec)
%     correlationThreshold - correlation below which window is bad (default = 0.4)
%     badTimeThreshold - cutoff fraction of bad corr windows (default = 0.01)
%     ransacSampleSize - samples for computing ransac (default = 50)
%     ransacChannelFraction - fraction of channels for robust reconstruction (default = 0.25)
%     ransacCorrelationThreshold - cutoff correlation for abnormal wrt neighbors(default = 0.75)
%     ransacUnbrokenTime - cutoff fraction of time channel can have poor ransac predictability (default = 0.4)
%     ransacWindowSeconds - correlation window for ransac (default = 5 sec)
%
% Output parameters (c channels, w windows):
%    ransacPerformed - true if there were enough good channels to do ransac
%    noisyChannels - list of identified bad channel numbers
%    badChannelsFromCorrelation  - list of bad channels identified by correlation
%    badChannelsFromDeviation   - list of bad channels identified by amplitude
%    badChannelsFromHFNoise - list of bad channels identified by SNR
%    badChannelsFromRansac - list of channels identified by ransac
%    fractionBadCorrelationWindows - c x 1 vector with fraction of bad correlation windows
%    robustChannelDeviation - c x 1 vector with robust measure of average channel deviation
%    zscoreHFNoise - c x 1 vector with measure of channel noise level
%    maximumCorrelations - w x c array with max window correlation
%    ransacCorrelations = c x wr array with ransac correlations
%
% This function uses 4 methods for detecting bad channels after removing
% from consideration channels that have NaN data or channels that are
% identically constant.
%
% Method 1: too low or high amplitude. If the z score of robust
%           channel deviation falls below robustDeviationThreshold, the channel is
%           considered to be bad.
% Method 2: too low an SNR. If the z score of estimate of signal above
%           50 Hz to that below 50 Hz above highFrequencyNoiseThreshold, the channel
%           is considered to be bad.
%
% Method 3: low correlation with other channels. Here correlationWindowSize is the window
%           size over which the correlation is computed. If the maximum
%           correlation of the channel to the other channels falls below
%           correlationThreshold, the channel is considered bad in that window.
%           If the fraction of bad correlation windows for a channel
%           exceeds badTimeThreshold, the channel is marked as bad.
%
% After the channels from methods 2 and 3 are removed, method 4 is
% computed on the remaining signals
%
% Method 4: each channel is predicted using ransac interpolation based
%           on a ransac fraction of the channels. If the correlation of
%           the prediction to the actual behavior is too low for too
%           long, the channel is marked as bad.
%
% Assumptions:
%  - The signal is a structure of continuous data with data, srate, chanlocs,
%    and chaninfo fields.
%  - The signal.data has been high pass filtered.
%  - No segments of the EEG data have been removed

% Methods 1 and 4 are adapted from code by Christian Kothe and Methods 2
% and 3 are adapted from code by Nima Bigdely-Shamlo
%
%% Check the incoming parameters
if nargin < 1
    error('findNoisyChannels:NotEnoughArguments', 'requires at least 1 argument');
elseif isstruct(signal) && ~isfield(signal, 'data')
    error('findNoisyChannels:NoDataField', 'requires a structure data field');
elseif size(signal.data, 3) ~= 1
    error('findNoisyChannels:DataNotContinuous', 'data must be a 2D array');
elseif nargin < 2 || ~exist('noisyIn', 'var') || isempty(noisyIn)
    noisyIn = struct();
end

%% Set the defaults and initialize as needed
noisyOut = getNoisyStructure();
defaults = getPrepDefaults(signal, 'reference');
[noisyOut, errors] = checkPrepDefaults(noisyIn, noisyOut, defaults);
if ~isempty(errors)
    error('findNoisyChannels:BadParameters', ['|' sprintf('%s|', errors{:})]);
end
%% Fix the channel locations
channelLocations = noisyOut.channelLocations;
evaluationChannels = sort(noisyOut.evaluationChannels); % Make sure channels are sorted
evaluationChannels = evaluationChannels(:)';            % Make sure row vector
noisyOut.evaluationChannels = evaluationChannels;
originalChannels = 1:size(signal.data, 1);

%% Extract the data required
data = signal.data;
originalNumberChannels = size(data, 1);          % Save the original channels
data = double(data(evaluationChannels, :))';      % Remove the unneeded channels
signalSize = size(data, 1);
correlationFrames = noisyOut.correlationWindowSeconds * signal.srate;
correlationWindow = 0:(correlationFrames - 1);
correlationOffsets = 1:correlationFrames:(signalSize-correlationFrames);
WCorrelation = length(correlationOffsets);
ransacFrames = noisyOut.ransacWindowSeconds*noisyOut.srate;
ransacWindow = 0:(ransacFrames - 1);
ransacOffsets = 1:ransacFrames:(signalSize-ransacFrames);
WRansac = length(ransacOffsets);
noisyOut.zscoreHFNoise = zeros(originalNumberChannels, 1);
noisyOut.noiseLevels = zeros(originalNumberChannels, WCorrelation);
noisyOut.maximumCorrelations = ones(originalNumberChannels, WCorrelation);
noisyOut.dropOuts = zeros(originalNumberChannels, WCorrelation);
noisyOut.correlationOffsets = correlationOffsets;
noisyOut.channelDeviations = zeros(originalNumberChannels, WCorrelation);
noisyOut.robustChannelDeviation = zeros(originalNumberChannels, 1);
noisyOut.ransacCorrelations = ones(originalNumberChannels, WRansac);
noisyOut.ransacOffsets = ransacOffsets;

%% Detect constant or NaN channels and remove from consideration
nanChannelMask = sum(isnan(data), 1) > 0;
noSignalChannelMask = mad(data, 1, 1) < 10e-10 | std(data, 1, 1) < 10e-10;
noisyOut.noisyChannels.badChannelsFromNaNs = evaluationChannels(nanChannelMask);
noisyOut.noisyChannels.badChannelsFromNoData = evaluationChannels(noSignalChannelMask);
evaluationChannels = setdiff(evaluationChannels, ...
    union(noisyOut.noisyChannels.badChannelsFromNaNs, ...
          noisyOut.noisyChannels.badChannelsFromNoData));
data = signal.data;
data = double(data(evaluationChannels, :))';
[signalSize, numberChannels] = size(data);

%% Method 1: Unusually high or low amplitude (using robust std)
channelDeviation = 0.7413 *iqr(data); % Robust estimate of SD
channelDeviationSD =  0.7413 * iqr(channelDeviation);
channelDeviationMedian = nanmedian(channelDeviation);
noisyOut.robustChannelDeviation(evaluationChannels) = ...
    (channelDeviation - channelDeviationMedian) / channelDeviationSD;

% Find channels with unusually high deviation
badChannelsFromDeviation = ...
    abs(noisyOut.robustChannelDeviation) > ...
             noisyOut.robustDeviationThreshold | ...
             isnan(noisyOut.robustChannelDeviation);
badChannelsFromDeviation = originalChannels(badChannelsFromDeviation);
noisyOut.noisyChannels.badChannelsFromDeviation = badChannelsFromDeviation(:)';
noisyOut.channelDeviationMedian = channelDeviationMedian;
noisyOut.channelDeviationSD = channelDeviationSD;

%% Method 2: Compute the SNR (based on Christian Kothe's clean_channels)
% Note: RANSAC uses the filtered values X of the data
if noisyOut.srate > 100
    % Remove signal content above 50Hz and below 1 Hz
    B = design_fir(100,[2*[0 45 50]/noisyOut.srate 1],[1 1 0 0]);
    X = zeros(signalSize, numberChannels);
    parfor k = 1:numberChannels  % Could be changed to parfor
        X(:,k) = filtfilt_fast(B, 1, data(:, k)); end
    % Determine z-scored level of EM noise-to-signal ratio for each channel
    noisiness = mad(data- X, 1)./mad(X, 1);
    noisinessMedian = nanmedian(noisiness);
    noisinessSD = mad(noisiness, 1)*1.4826;
    zscoreHFNoiseTemp = (noisiness - noisinessMedian) ./ noisinessSD;
    noiseMask = (zscoreHFNoiseTemp > noisyOut.highFrequencyNoiseThreshold) | ...
                isnan(zscoreHFNoiseTemp);
    % Remap channels to original numbering
    badChannelsFromHFNoise  = evaluationChannels(noiseMask);
    noisyOut.noisyChannels.badChannelsFromHFNoise = badChannelsFromHFNoise(:)';
else
    X = data;
    noisinessMedian = 0;
    noisinessSD = 1;
    zscoreHFNoiseTemp = zeros(numberChannels, 1);
    noisyOut.noisyChannels.badChannelsFromHFNoise = [];
end

% Remap the channels to original numbering for the zscoreHFNoise
noisyOut.zscoreHFNoise(evaluationChannels) = zscoreHFNoiseTemp;
noisyOut.noisinessMedian = noisinessMedian;
noisyOut.noisinessSD = noisinessSD;

%% Method 3: Global correlation criteria (from Nima Bigdely-Shamlo)
channelCorrelations = ones(WCorrelation, numberChannels);
noiseLevels = zeros(WCorrelation, numberChannels);
channelDeviations = zeros(WCorrelation, numberChannels);
n = length(correlationWindow);
xWin = reshape(X(1:n*WCorrelation, :)', numberChannels, n, WCorrelation);
dataWin = reshape(data(1:n*WCorrelation, :)', numberChannels, n, WCorrelation);
parfor k = 1:WCorrelation
    eegPortion = squeeze(xWin(:, :, k))';
    dataPortion = squeeze(dataWin(:, :, k))';
    windowCorrelation = corrcoef(eegPortion);
    abs_corr = abs(windowCorrelation - diag(diag(windowCorrelation)));
    channelCorrelations(k, :)  = quantile(abs_corr, 0.98);
    noiseLevels(k, :) = mad(dataPortion - eegPortion, 1)./mad(eegPortion, 1);
    channelDeviations(k, :) =  0.7413 *iqr(dataPortion);
end
dropOuts = isnan(channelCorrelations) | isnan(noiseLevels);
channelCorrelations(dropOuts) = 0.0;
noiseLevels(dropOuts) = 0.0;
clear xWin;
clear dataWin;
noisyOut.maximumCorrelations(evaluationChannels, :) = channelCorrelations';
noisyOut.noiseLevels(evaluationChannels, :) = noiseLevels';
noisyOut.channelDeviations(evaluationChannels, :) = channelDeviations';
noisyOut.dropOuts(evaluationChannels, :) = dropOuts';
thresholdedCorrelations = ...
    noisyOut.maximumCorrelations < noisyOut.correlationThreshold;
fractionBadCorrelationWindows = mean(thresholdedCorrelations, 2);
fractionBadDropOutWindows = mean(noisyOut.dropOuts, 2);

% Remap channels to their original numbers
badChannelsFromCorrelation = find(fractionBadCorrelationWindows > noisyOut.badTimeThreshold);
noisyOut.noisyChannels.badChannelsFromCorrelation = badChannelsFromCorrelation(:)';
badChannelsFromDropOuts = find(fractionBadDropOutWindows > noisyOut.badTimeThreshold);
noisyOut.noisyChannels.badChannelsFromDropOuts = badChannelsFromDropOuts(:)';
noisyOut.medianMaxCorrelation =  median(noisyOut.maximumCorrelations, 2);

%% Bad so far by amplitude and correlation (take these out before doing ransac)
noisyChannels = union(noisyOut.noisyChannels.badChannelsFromDeviation, ...
    union(noisyOut.noisyChannels.badChannelsFromCorrelation, ...
          noisyOut.noisyChannels.badChannelsFromDropOuts));

%% Method 4: Ransac corelation (may not be performed)
% Setup for ransac (if a 2-stage algorithm, remove other bad channels first)
if noisyOut.ransacOff
    noisyOut.ransacBadWindowFraction = 0;
    noisyOut.ransacPerformed = false;
elseif isempty(channelLocations)
    warning('findNoisyChannels:noChannelLocation', ...
        'ransac could not be computed because there were no channel locations');
    noisyOut.ransacBadWindowFraction = 0;
    noisyOut.ransacPerformed = false;
else % Set up parameters and make sure enough good channels to proceed
    [ransacChannels, idiff] = setdiff(evaluationChannels, noisyChannels);
    X = X(:, idiff);

    % Calculate the parameters for ransac
    ransacSubset = round(noisyOut.ransacChannelFraction*size(data, 2));
    if noisyOut.ransacUnbrokenTime < 0
        error('find_noisyChannels:BadUnbrokenParameter', ...
            'ransacUnbrokenTime must be greater than 0');
    elseif noisyOut.ransacUnbrokenTime < 1
        ransacUnbrokenFrames = signalSize*noisyOut.ransacUnbrokenTime;
    else
        ransacUnbrokenFrames = srate*noisyOut.ransacUnbrokenTime;
    end

    nchanlocs = channelLocations(ransacChannels);
    if length(nchanlocs) ~= size(nchanlocs, 2)
        nchanlocs = nchanlocs';
    end
    if length(nchanlocs) < ransacSubset + 1 || length(nchanlocs) < 3 || ...
            ransacSubset < 2
        warning('find_noisyChannels:NotEnoughGoodChannels', ...
            'Too many channels have failed quality tests to perform ransac');
        noisyOut.ransacBadWindowFraction = 0;
        noisyOut.ransacPerformed = false;
    end
end

if noisyOut.ransacPerformed
    try
    % Calculate all-channel reconstruction matrices from random channel subsets
       locs = [cell2mat({nchanlocs.X}); cell2mat({nchanlocs.Y});cell2mat({nchanlocs.Z})];
    catch err
       error('findNoisyChannels:NoXYZChannelLocations', ...
             'Must provide valid channel locations');
    end
    if isempty(locs) || size(locs, 2) ~= length(ransacChannels) ...
            || any(isnan(locs(:)))
          error('find_noisyChannels:EmptyChannelLocations', ...
            'The signal chanlocs must have valid X, Y, and Z components');
    end
    P = hlp_microcache('cleanchans', @calc_projector, locs, ...
        noisyOut.ransacSampleSize, ransacSubset);
    ransacCorrelationsT = zeros(length(locs), WRansac);

    % Calculate each channel's correlation to its RANSAC reconstruction for each window
    n = length(ransacWindow);
    m = length(ransacChannels);
    p = noisyOut.ransacSampleSize;
    Xwin = reshape(X(1:n*WRansac, :)', m, n, WRansac);
    parfor k = 1:WRansac
        ransacCorrelationsT(:, k) = ...
            calculateRansacWindow(squeeze(Xwin(:, :, k))', P, n, m, p);
    end
    clear Xwin;
    noisyOut.ransacCorrelations(ransacChannels, :) = ransacCorrelationsT;
    flagged = noisyOut.ransacCorrelations < noisyOut.ransacCorrelationThreshold;
    badChannelsFromRansac = find(sum(flagged, 2)*ransacFrames > ransacUnbrokenFrames)';
    noisyOut.noisyChannels.badChannelsFromRansac = badChannelsFromRansac(:)';
    noisyOut.ransacBadWindowFraction = sum(flagged, 2)/size(flagged, 2);
end

% Combine bad channels detected from all methods
noisy = noisyOut.noisyChannels;
noisyOut.noisyChannels.badChannelsFromLowSNR = ...
    intersect(noisy.badChannelsFromHFNoise, noisy.badChannelsFromCorrelation);
noisyChannels = union(noisyChannels, ...
    union(union(noisy.badChannelsFromRansac, ...
          noisy.badChannelsFromHFNoise), ...
    union(noisy.badChannelsFromNaNs, ...
          noisy.badChannelsFromNoData)));
noisyOut.noisyChannels.all = noisyChannels(:)';
noisyOut.medianMaxCorrelation =  median(noisyOut.maximumCorrelations, 2);

%% Helper functions for findNoisyChannels
function P = calc_projector(locs, numberSamples, subsetSize)
% Calculate a bag of reconstruction matrices from random channel subsets

[permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples);
randomSamples = cell(1, numberSamples);
parfor k = 1:numberSamples
    tmp = zeros(size(locs, 2));
    slice = subsets(k, :);
    tmp(slice, :) = real(spherical_interpolate(permutedLocations(:, :, k), locs))';
    randomSamples{k} = tmp;
end
P = horzcat(randomSamples{:});

function [permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples)
 stream = RandStream('mt19937ar', 'Seed', 435656);
 numberChannels = size(locs, 2);
 permutedLocations = zeros(3, subsetSize, numberSamples);
 subsets = zeros(numberSamples, subsetSize);
 for k = 1:numberSamples
     subset = randsample(1:numberChannels, subsetSize, stream);
     subsets(k, :) = subset;
     permutedLocations(:, :,  k) = locs(:, subset);
 end

function Y = randsample(X, num, stream)
Y = zeros(1, num);
for k = 1:num
    pick = round(1 + (length(X)-1).*rand(stream));
    Y(k) = X(pick);
    X(pick) = [];
end

function rX = calculateRansacWindow(XX, P, n, m, p)
    YY = sort(reshape(XX*P, n, m, p),3);
    YY = YY(:, :, round(end/2));
    rX = sum(XX.*YY)./(sqrt(sum(XX.^2)).*sqrt(sum(YY.^2)));

function noisyOut = getNoisyStructure()
    noisyOut = struct('srate', [], ...
        'samples', [], ...
        'evaluationChannels', [], ...
        'channelLocations', [], ...
        'robustDeviationThreshold', [], ...
        'highFrequencyNoiseThreshold', [], ...
        'correlationWindowSeconds', [], ...
        'correlationThreshold', [], ...
        'badTimeThreshold', [], ...
        'ransacSampleSize', [], ...
        'ransacChannelFraction', [], ...
        'ransacCorrelationThreshold', [], ...
        'ransacUnbrokenTime', [], ...
        'ransacWindowSeconds', [], ...
        'noisyChannels', getBadChannelStructure(), ...
        'ransacPerformed', true, ...
        'channelDeviationMedian', [], ...
        'channelDeviationSD', [], ...
        'channelDeviations', [], ...
        'robustChannelDeviation', [], ...
        'noisinessMedian', [], ...
        'noisinessSD', [], ...
        'zscoreHFNoise', [], ...
        'noiseLevels', [], ...
        'maximumCorrelations', [], ...
        'dropOuts', [], ...
        'medianMaxCorrelation', [], ...
        'correlationOffsets', [], ...
        'ransacCorrelations', [], ...
        'ransacOffsets', [], ...
        'ransacBadWindowFraction', []);

author: VisLab

repo: EEG-Clean-Tools


text:
% OP_GRADSYMU_GRADSYMV_MP: assemble the stiffness matrix A = [a(i,j)], a(i,j) = (epsilon gradsym u_j, gradsym v_i), in a multipatch domain.
%
%   mat = op_gradsymu_gradsymv_mp (spu, spv, msh, [epsilon], [patches]);
%
% INPUT:
%
%   spu:     object representing the space of trial functions (see sp_multipatch)
%   spv:     object representing the space of test functions (see sp_multipatch)
%   msh:     object defining the domain partition and the quadrature rule (see msh_multipatch)
%   epsilon: function handle to compute the diffusion coefficient. Equal to one if left empty.
%   patches: list of patches where the integrals have to be computed. By default, all patches are selected.
%
% OUTPUT:
%
%   mat:    assembled stiffness matrix
%
% Copyright (C) 2015, 2016, 2020 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function A = op_gradsymu_gradsymv_mp (spu, spv, msh, coeff, patch_list)

  if (nargin < 5)
    patch_list = 1:msh.npatch;
  end

  if ((spu.npatch ~= spv.npatch) || (spu.npatch ~= msh.npatch))
    error ('op_gradsymu_gradsymv_mp: the number of patches does not coincide')
  end

  ncounter = 0;
  for iptc = patch_list
    if (nargin < 4 || isempty (coeff))
      [rs, cs, vs] = op_gradsymu_gradsymv_tp (spu.sp_patch{iptc}, spv.sp_patch{iptc}, msh.msh_patch{iptc});
    else
      [rs, cs, vs] = op_gradsymu_gradsymv_tp (spu.sp_patch{iptc}, spv.sp_patch{iptc}, msh.msh_patch{iptc}, coeff);
    end
    rows(ncounter+(1:numel (rs))) = spv.gnum{iptc}(rs);
    cols(ncounter+(1:numel (rs))) = spu.gnum{iptc}(cs);

    if (~isempty (spv.dofs_ornt))
      vs = spv.dofs_ornt{iptc}(rs)' .* vs;
    end
    if (~isempty (spu.dofs_ornt))
      vs = vs .* spu.dofs_ornt{iptc}(cs)';
    end

    vals(ncounter+(1:numel (rs))) = vs;
    ncounter = ncounter + numel (rs);
  end

  A = sparse (rows, cols, vals, spv.ndof, spu.ndof);
  clear rows cols vals rs cs vs

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
% icp_demo.m
% run
% >> icp_demo
% in MATLAB

% size of model/data
n_model=2000;
n_data=1300;

% model points
model=4*randn(3,n_model)-2*ones(3,n_model);
bol=(model(1,:).^2+model(2,:).^2)<2^2;
while any(not(bol))
    model(:,not(bol))=4*randn(3,sum(not(bol)))-2*ones(3,sum(not(bol)));
    bol=(model(1,:).^2+model(2,:).^2)<2^2;
end
model(:,not(bol))=4*randn(3,sum(not(bol)))-2*ones(3,sum(not(bol)));
model(3,:)=0.5*(model(1,:).^2-model(2,:).^2);
bol=(model(1,:).^2+model(2,:).^2)<0.9^2;
model(3,bol)=3.5-3.5*(1/0.9)*sqrt(model(1,bol).^2+model(2,bol).^2);

% data points
data=model(:,1:n_data);

% weights
weights=ones(1,n_data);weights=rand(1,n_data);

% Transform points in data (move away from start position).
deg=15;
v1=2*(rand-0.5)*deg*(pi/180);
v2=2*(rand-0.5)*deg*(pi/180);
v3=2*(rand-0.5)*deg*(pi/180);
TR0=[cos(v1) sin(v1) 0;-sin(v1) cos(v1) 0;0 0 1];
TR0=TR0*[cos(v2) 0 sin(v2);0 1 0;-sin(v2) 0 cos(v2)];
TR0=TR0*[1 0 0;0 cos(v3) sin(v3);0 -sin(v3) cos(v3)];
TT0=3.0*(rand(3,1)-0.5);
data=TR0*data+repmat(TT0,1,n_data);

% randvec
rndvec=uint32(randperm(n_data)-1);

% sizerand, number of point-matchings in each iteration.
sizernd=ceil(1.45*n_data);

% Number of iterations.
iter=40;

% Compile the c++ files (not nessesacily if it is already done).
try
    make
end

% Create the kd-tree, TreeRoot is the pointer to the kd-tree
[tmp, tmp, TreeRoot] = kdtree( model', []);

% Run the ICP algorithm.
[R,T]=icpCpp(model,data,weights,rndvec,sizernd,TreeRoot,iter);

% Free allocated memory for the kd-tree.
kdtree([],[],TreeRoot);

figure(1),plot3(model(1,:),model(2,:),model(3,:),'r.',data(1,:),data(2,:),data(3,:),'c.'),hold off;

% Transform the data points.
data=R*data+repmat(T,1,n_data);

figure(2),plot3(model(1,:),model(2,:),model(3,:),'r.',data(1,:),data(2,:),data(3,:),'b.'),hold off;

clear functions

author: Sable

repo: mcbench-benchmarks


text:
classdef nme_branch_acc < mp.nme_branch_ac & mp.form_acc

%   MATPOWER
%   Copyright (c) 2018-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Baljinnyam Sereeter, Delft University of Technology
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    methods
        function [h, dh] = ang_diff_fcn(obj, xx, Aang, lang, uang)
            if nargout > 1
                [va, dva] = obj.va_fcn(xx, [], 0);
                Adva = Aang * dva;
                dh = [ -Adva; Adva ];
            else
                va = obj.va_fcn(xx, [], 0);
            end
            Ax = Aang * va;
            h = [ lang - Ax;
                  Ax - uang ];
        end

        function d2H = ang_diff_hess(obj, xx, lambda, Aang)
            %% unpack data
            [vr, vi] = deal(xx{:});
            nn = length(vr);

            %% evaluate Hessian of voltage magnitude^2 function
            n = length(lambda) / 2;
            if n
                lam = lambda(n+1:2*n) - lambda(1:n);    %% lam_ub - lam_lb
            else
                lam = zeros(0,1);
            end

            d2H = obj.va_hess(xx, Aang' * lam, []);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ozf] = kN2ozf(kN)
% Convert force from kilonewtons to ounces-force.
% Chad A. Greene 2012
ozf = kN*3596.9431019;


text:
% anova1_cell() - compute F-values in cell array using ANOVA.
%
% Usage:
%    >> [F df] = anova1_cell( data );
%
% Inputs:
%   data       = data consisting of PAIRED arrays to be compared. The last
%                dimension of the data array is used to compute ANOVA.
% Outputs:
%   F   - F-value
%   df  - degree of freedom (array)
%
% Note: the advantage over the ANOVA1 function of Matlab statistical
%       toolbox is that this function works on arrays (see examples). Note
%       also that you still need the statistical toolbox to assess
%       significance using the fcdf() function. The other advantage is that
%       this function will work with complex numbers.
%
% Example:
%   a = { rand(1,10) rand(1,10) rand(1,10) }
%   [F df] = anova1_cell(a)
%   signif = 1-fcdf(F, df(1), df(2))
%
%   % for comparison
%   anova1( [ a{1,1}' a{1,2}' a{1,3}' ]) % look in the graph for the F value
%
%   b = { [ a{1,1}; a{1,1} ] [ a{1,2}; a{1,2} ] [ a{1,3}; a{1,3} ] }
%   [F df] = anova1_cell(b)
%
%   c{1,1} = reshape(repmat(b{1,1}, [2 1]),2,2,10);
%   c{1,2} = reshape(repmat(b{1,2}, [2 1]),2,2,10);
%   c{1,3} = reshape(repmat(b{1,3}, [2 1]),2,2,10);
%   [F df] = anova1_cell(c)
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005
%
% Reference:
%   Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill.

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

function [F, df] = anova1_cell(data)

    % This function does not return
    % correct values (see bug 336)
    % It should be fixed with Schaum's outlines p363
    % but requires some work. It now calls
    % anova2_cell which returns correct values

    warning off;
    [ F tmp tmp2 df] =  anova2_cell(data);
    warning on;
    return;

    % compute all means and all std
    % -----------------------------
    nd = myndims( data{1} );
    if nd == 1

        for i = 1:length(data)
            n( i) = length(data{i});
            m( i) = mymean(  data{i});
            sd(i) = mystd(   data{i});
        end;
        nt = sum(n);
        n   = n';
        m   = m';
        sd  = sd';

    elseif nd == 2

        for i = 1:length(data)
            n( :,i) = ones(size(data{i},1),1) * size(data{i},2, 'single');
            m( :,i) = mymean(  data{i},2);
            sd(:,i) = mystd(   data{i},[],2);
        end;
        nt = sum(n(1,:));

    elseif nd == 3

        for i = 1:length(data)
            n( :,:,i) = ones(size(data{i},1),size(data{i},2)) * size(data{i},3, 'single');
            m( :,:,i) = mymean(  data{i},3);
            sd(:,:,i) = mystd(   data{i},[],3);
        end;
        nt = sum(n(1,1,:));

    else

        for i = 1:length(data)
            n( :,:,:,i) = ones(size(data{i},1),size(data{i},2), size(data{i},3)) * size(data{i},4, 'single');
            m( :,:,:,i) = mymean(  data{i},4);
            sd(:,:,:,i) = mystd(   data{i},[],4);
        end;
        nt = sum(n(1,1,1,:));

    end;

    mt = mean(m,nd);
    ng = length(data); % number of conditions

    VinterG  = ( sum( n.*(m.^2), nd ) - nt*mt.^2 )/(ng-1);
    VwithinG = sum( (n-1).*(sd.^2), nd )/(nt-ng);
    F  = VinterG./VwithinG;
    df = [ ng-1 ng*(size(data{1},nd)-1) ];

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

function res = mymean( data, varargin) % deal with complex numbers
    res = mean( data, varargin{:});
    if ~isreal(data)
        res = abs( res );
    end;

function res = mystd( data, varargin) % deal with complex numbers
    res = std( abs(data), varargin{:});

author: PatternRecognition

repo: OpenBMI


text:
function pascal_to_i4_test ( )

%*****************************************************************************80
%
%% PASCAL_TO_I4_TEST tests PASCAL_TO_I4.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/pascal_to_i4.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PASCAL_TO_I4_TEST\n' );
  fprintf ( 1, '  PASCAL_TO_I4 converts Pascal triangle indices to a\n' );
  fprintf ( 1, '  linear index.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     J =>    K\n' );
  fprintf ( 1, '\n' );

  for d = 0 : 4
    for i = d : -1 : 0
      j = d - i;
      k = pascal_to_i4 ( i, j );
      fprintf ( 1, '  %4d  %4d    %4d\n', i, j, k );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_suite = test_regression_additive2

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION_ADDITIVE2
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testNeuralNetworkCFPrediction(testCase)
  verifyVarsEqual(testCase, getName(), {'Eft_map'}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [isBehavior] = bz_isBehavior(behavior)
% USAGE
% [isBehavior] = bz_isBehavior(behavior)
%
% INPUT
%       behavior   - struct with the following fields
%                   .timestamps
%                   .position
%                       .detectorname
%                       .detectionparms
%                       .detectioninterval
%                       .detectiondate
%
% OUTPUT
%      logical true if struct meets behavior criteria, false if otherwise
%
% written by david tingley, 2017

if isfield(behavior,'behaviorinfo') && isfield(behavior,'timestamps') && isfield(behavior,'samplingRate')% check that fields exist
     if isstruct(behavior.behaviorinfo) && isvector(behavior.timestamps) && isnumeric(behavior.samplingRate)
         % check that sub fields exist
         if isfield(behavior.behaviorinfo,'description') && isfield(behavior.behaviorinfo,'acquisitionsystem') && isfield(behavior.behaviorinfo,'substructnames')
            isBehavior = true;
         else
            warning('behavior struct seems to be missing something in .behaviorinfo')
            isBehavior = false;
         end

     else
         warning('one of the required fields for an behavior type is not formatted correctly')
         isBehavior = false;
     end
else
    warning('one of the required fields for an behavior type does not exist')
    isBehavior = false;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [lp,dlp] = priorT(mu,s2,nu,x)

% Univariate Student's t hyperparameter prior distribution.
% Compute log-likelihood and its derivative or draw a random sample.
% The prior distribution is parameterized as:
%
%   p(x) = ( 1 + (x-mu)^2/(s2*(nu-2)) )^(-(nu+1)/2) / Z, where
%      Z = gamma(nu/2) * sqrt(pi*(nu-2)*s2) / gamma((nu+1)/2),
%
% mu(1x1) is the mean parameter, s2(1x1) is the variance parameter,
% nu(1x1) > 2 is the degrees of freedom parameter and x(1xN) contains query
% hyperparameters for prior evaluation.
%
% For more help on design of priors, try "help priorDistributions".
%
% Copyright (c) by Roman Garnett and Hannes Nickisch, 2014-10-16.
%
% See also PRIORDISTRIBUTIONS.M.

if nargin<3, error('mu, s2 and nu parameters need to be provided'), end
if ~(isscalar(mu)&&isscalar(s2)&&isscalar(nu))
  error('mu, s2 and nu parameters need to be scalars')
end
if nargin<4, lp = mu + sqrt(s2*(nu-2)/priorGamma(nu/2,2))*randn; return, end

lZ = gammaln((nu+1)/2)-gammaln(nu/2)-log(pi*(nu-2)*s2)/2;
lp = -(nu+1)/2*log(1+(x-mu).^2/(s2*(nu-2))) + lZ;
dlp = (nu+1)*(mu-x)./(mu^2-2*mu*x+s2*(nu-2)+x.^2);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function C=multistep(x,N,ld)
x=im2double(x);
[m n]=size(x);
C=x;
for k=0:1:N-1
    R=getcorner(C,m/2^k,n/2^k);
    S=wave_transform(R,ld);
    C=putcorner(C,S);
end
end
function r=getcorner(C,m,n)
    r=C(1:m,1:n);
end
function x=putcorner(c,a)
    [m n]=size(a);
    c(1:m,1:n)=a;
    x=c;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Script chap::3::script
% Density Variance Gamma model with Cox-Ingersoll-Ross stochastic clock
%
%
%
T = 1;              % maturity
f0 = 100;           % spot value
r=0;                % risk less rate
d=0;                % dividend yield
% parameters for applying fourier transform to compute the density
a = 600;                            % parameter to compute density via fft
N = 1024;                           % number of grid points
x = ( (0:N-1) - N/2 ) / a;          % range

C = 2;
G = 2;
M = 2;
Y = .75;

legend = 'Base';
title_plot = 'CGMY Density';

func = @(x) cf_cgmy(x,T,0,r,d, C, G, M, Y);
y = fftdensity(func,a,N);
%% Changing C
C_low = 1;
C_high = 3;

func_low = @(x) cf_cgmy(x,T,0,r,d, C_low,G,M,Y);
y_low = fftdensity(func_low,a,N);
func_high = @(x) cf_cgmy(x,T,0,r,d, C_high,G,M,Y);
y_high = fftdensity(func_high,a,N);
legend_low = 'Changing C low';
legend_high = 'Changing C high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);

%% Changing G
G_low = 1.5;
G_high = 2.5;

func_low = @(x) cf_cgmy(x,T,0,r,d, C,G_low,M,Y);
y_low = fftdensity(func_low,a,N);
func_high = @(x) cf_cgmy(x,T,0,r,d, C,G_high,M,Y);
y_high = fftdensity(func_high,a,N);
legend_low = 'Changing G low';
legend_high = 'Changing G high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);

%% Changing M
M_low = 1.5;
M_high = 2.5;

func_low = @(x) cf_cgmy(x,T,0,r,d, C,G,M_low,Y);
y_low = fftdensity(func_low,a,N);
func_high = @(x) cf_cgmy(x,T,0,r,d, C,G,M_high,Y);
y_high = fftdensity(func_high,a,N);
legend_low = 'Changing M low';
legend_high = 'Changing M high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);

%% Changing Y
Y_low = -0.5;
Y_high = 1.2;

func_low = @(x) cf_cgmy(x,T,0,r,d, C,G,M,Y_low);
y_low = fftdensity(func_low,a,N);
func_high = @(x) cf_cgmy(x,T,0,r,d, C,G,M,Y_high);
y_high = fftdensity(func_high,a,N);
legend_low = 'Changing Y low';
legend_high = 'Changing Y high';

createfigure_density(x,y,y_low,y_high,title_plot,legend,legend_low,legend_high);


text:
function cvt_test02 ( )

%*****************************************************************************80
%
%% CVT_TEST02 repeats test 1, but uses twice as many iterations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  CVT computes a Centroidal Voronoi Tessellation.\n' );
  fprintf ( 1, '  Repeat test 1, but with twice the number of iterations.\n' );

  dim_num = 2;
  n = 10;
  batch = 1000;
  init = 0;
  init_string = 'uniform';
  it_max = 80;
  it_fixed = 1;
  sample = 0;
  sample_num = 10000;
  sample_string = 'uniform';
  seed = 123456789;
  r = [];

  seed_init = seed;

  [ r, seed, it_num, it_diff, energy ] = cvt ( dim_num, n, batch, init, ...
    sample, sample_num, it_max, it_fixed, seed, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dimension DIM_NUM =        %12d\n', dim_num );
  fprintf ( 1, '  Number of points N =       %12d\n', n );
  fprintf ( 1, '  Initial SEED =             %12d\n', seed_init );
  fprintf ( 1, '  Current SEED =             %12d\n', seed );
  fprintf ( 1, '  INIT =                    "%s".\n', init_string );
  fprintf ( 1, '  Max iterations IT_MAX =    %12d\n', it_max );
  fprintf ( 1, '  IT_FIXED (fixed samples) = %12d\n', it_fixed );
  fprintf ( 1, '  Iterations IT_NUM =        %12d\n', it_num );
  fprintf ( 1, '  Difference IT_DIFF =       %14f\n', it_diff );
  fprintf ( 1, '  CVT ENERGY =               %14f\n', energy );
  fprintf ( 1, '  SAMPLE =                  "%s".\n', sample_string );
  fprintf ( 1, '  Samples SAMPLE_NUM    =    %12d\n', sample_num );
  fprintf ( 1, '  Sampling BATCH size =      %12d\n', batch );
  fprintf ( 1, '  EPSILON (unit roundoff) =  %12e\n', eps );

  r8mat_transpose_print ( dim_num, n, r, '  Generators (rows):' );

  return
end


text:
filename='BridgeCool_Quadrilateral_Bilinear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI'; incrementFactor = 1;
filterType = 'P1';
line_search_initiator = 'STANDARD';

HJiter0 = 1;
e2 = 100;
N_holes = [24 11 11];
R_holes = 0.4;
phase_holes = [0 0 0];
nsteps = 5;

Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 0.4;
optimality_initial = 5e-2;
constr_initial = 5e-2;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% DEMSPGP1DGP5 Do a simple 1-D regression after Snelson & Ghahramani's example.

% GP

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'spgp1d';
experimentNo = 5;

% load data
[X, y] = mapLoadData(dataSetName);

% Set up model
options = gpOptions('dtcvar');
options.numActive = 9;
options.optimiser = 'scg';

% use the deterministic training conditional.
q = size(X, 2);
d = size(y, 2);

model = gpCreate(q, d, X, y, options);
model.X_u = randn(9, 1)*0.25 - 0.75;
params = gpExtractParam(model);
model = gpExpandParam(model, params);

% Optimise the model.
iters = 1000;
display = 1;

model = gpOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

demSpgp1dPlot

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function s = ismatch(x, y)

% ISMATCH returns true if x is a member of array y, regardless of the class
% of x and y, if y is a string, or a cell-array of strings, it can contain
% the wildcard '*'

if isempty(x) || isempty(y)
  s = false;
elseif ischar(x) && ischar(y)
  y = sprintf('%s%s%s', '^', regexptranslate('wildcard',y), '$');
  s = ~isempty(regexp(x, y, 'once'));
elseif isnumeric(x) && isnumeric(y)
  s = ismember(x, y);
elseif ischar(x) && iscell(y)
  y = y(strcmp(class(x), cellfun(@class, y, 'UniformOutput', false)));
  s = ismember(x, y);

  % one or more of the elements in y can contain a wildcard, only proceed if
  % s=false
  if ~s && any(contains(y, '*'))
    y = y(contains(y, '*'));
    for i = 1:numel(y)
      tmpy = sprintf('%s%s%s', '^', regexptranslate('wildcard',y{i}), '$');
      s = ~isempty(regexp(x, tmpy, 'once'));
      if s, return; end
    end
  end
elseif isnumeric(x) && iscell(y) && all(cellfun(@isnumeric, y))
  s = false;
  for i=1:numel(y)
    s = s || ismember(x, y{i});
  end
else
  s = false;
end

author: spm


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% We will show how to acquire images from movie files.

% First we create an object to connect us to the image source which in
% this example is a movie file, a standard MPG format movie file
camera = Movie('traffic_sequence.mpg');
about camera
% and we see that the result is a Movie object

% The displayed value of this object
camera
% shows the frame dimension, frame rate, number of frames etc.

% To grab a frame from the movie
im = camera.grab();
% and then display it
idisp(im)

% We can display subsequent frames from the movie (note the traffic moving)
camera.grab();
camera.grab();
camera.grab();
camera.grab();
% Note that if we call the grab() method without a return argument the image
% is automatically displayed.

% Finally we destroy the camera object which deallocates resources
clear camera

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ o, c, e ] = lpp_to_polynomial ( m, l, o_max )

%*****************************************************************************80
%
%% LPP_TO_POLYNOMIAL writes a Legendre Product Polynomial as a polynomial.
%
%  Discussion:
%
%    For example, if
%      M = 3,
%      L = ( 1, 0, 2 ),
%    then
%      L(1,0,2)(X,Y,Z)
%      = L(1)(X) * L(0)(Y) * L(2)(Z)
%      = X * 1 * ( 3Z^2-1)/2
%      = - 1/2 X + (3/2) X Z^2
%    so
%      O = 2 (2 nonzero terms)
%      C = -0.5
%           1.5
%      E = 4    <-- index in 3-space of exponent (1,0,0)
%          15   <-- index in 3-space of exponent (1,0,2)
%
%    The output value of O is no greater than
%      O_MAX = product ( 1 <= I <= M ) (L(I)+2)/2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer L(M), the index of each Legendre product polynomial factor.
%    0 <= L(*).
%
%    Input, integer O_MAX, an upper limit on the size of the output arrays.
%      O_MAX = product ( 1 <= I <= M ) (L(I)+2)/2.
%
%    Output, integer O, the "order" of the polynomial product.
%
%    Output, real C(O), the coefficients of the polynomial product.
%
%    Output, integer E(O), the indices of the exponents of the
%    polynomial product.
%
  o1 = 1;
  c1 = 1.0;
  e1 = 1;
%
%  Implicate one factor at a time.
%
  for i = 1 : m

    [ o2, c2, f2 ] = lp_coefficients ( l(i) );

    o = 0;

    for j2 = 1 : o2
      for j1 = 1 : o1
        o = o + 1;
        c(o) = c1(j1) * c2(j2);

        if ( 1 < i )
          p = mono_unrank_grlex ( i - 1, e1(j1) );
        end
        p(i) = f2(j2);
        e(o) = mono_rank_grlex ( i, p );

      end
    end

    [ c, e ] = polynomial_sort ( o, c, e );
    [ o, c, e ] = polynomial_compress ( o, c, e );

    o1 = o;
    c1 = c;
    e1 = e;

  end

  return
end


text:
function [ cnt, S, bad_idx ] = paramStats( annotation_struct, options, verbose )
%PARAMSTATS Collects statistics about possible positions of each landmark
%point (component)
%   Detailed explanation goes here
%
% INPUT:
%   image      ... array of cells, each cell containes bbox found by detector and filename
%   bw         ... vector [w; h] contains size of base window.
%   components ... matrix of size 2xM, M is count of components [w0 ... wM-1; h0 ... hM-1]. Order of components is important!
%   bw_margin  ... vector [w; h] base window margin in percents of original bbox.
%   image_path ... path to directory containing image database. Must contain subdir /mgt with xml annotations.
%   verbose    ... true for image output with landmarks and components depicted
%
% OUTPUT:
%   cnt     ... count of images that passed test
%   S       ... matrix 4xM [x0_min; y0_min; x0_max; y0_max; ... xM-1_min; yM-1_min; xM-1_max; yM-1_max] bbox for each component
%   bad_idx ... indices to image array of images that does not passed this
%               test
%
% 16-07-10 Michal Uricar
% 10-04-11 Michal Uricar, estimation of S0 fixed (face center intead of nose)
% 12-07-11 Michal Uricar, corners dataset
% 21-03-12 Michal Uricar, LFW annotation in one file

    if (nargin < 3)
        verbose = false;
    end;

    bad_idx = [];

    % count of components
    M = size(options.components, 2);

    cnt = 0;
    S = [inf(2, M); -inf(2, M)];

    % for each image in db
    N = annotation_struct.N;
    for i = 1 : N
        [Iframe, Annotation, I, itmp.bbox, bbox, OrigPoints] = getImageFrame(options, i, annotation_struct);

        % !!! Annotation is empty at i == 578 !!!
        if (isempty(Annotation))
            bad_idx = [bad_idx i];
            continue;
        end;

        Points = prepareS0gt(Annotation.P, options);    % transform nose to the center of face
        Points = Points(:, options.comselect);          % extract relevant points only (name list in options.compnames)
        Points(:, M) = Annotation.P(:, 10);             % copy back original nose position

        %% Check for each component if it fits in normalized frame
        flag = true;
        for j = 1 : M
            if ( ((Points(1, j) - options.components(1, j)/2) < 0) || ((Points(1, j) + options.components(1, j)/2) > options.bw(1)) || ...
                 ((Points(2, j) - options.components(2, j)/2) < 0) || ((Points(2, j) + options.components(2, j)/2) > options.bw(2)) )
                flag = false;
            end;
        end;

        % increase number of valid images if flag is true
        if (flag)
            cnt = cnt + 1;
            fprintf('%.2f%% - Passed: tested image no.%d file %s...\n', i*100/N, i, Annotation.image.filename);

            %% recompute AABB for S
            for j = 1 : M
                bb = [Points(1, j) - options.components(1, j)/2 Points(2, j) - options.components(2, j)/2 ...
                      Points(1, j) + options.components(1, j)/2 Points(2, j) + options.components(2, j)/2 ];
                S(1, j) = min(S(1, j), bb(1));
                S(2, j) = min(S(2, j), bb(2));
                S(3, j) = max(S(3, j), bb(3));
                S(4, j) = max(S(4, j), bb(4));
            end;

        else
            fprintf('%.2f%% - Not passed: tested image no.%d file %s...\n', i*100/N, i, Annotation.image.filename);
            bad_idx = [bad_idx i];
        end;

        %% Visualization
        if (verbose)
            aabb = makeAABB(itmp.bbox);
            aabb2 = makeAABB(bbox);

            % show original image with annotation
            figure(1);
            subplot(2, 2, 1);
            imshow(I, []); hold on;
            % show ground truth points
%             for k = 1 : M
            plot(OrigPoints(1, :), OrigPoints(2, :), 'r.', 'LineWidth', 2, 'MarkerSize', 10);
%             end;
            % show aabb found by detector
            line([aabb(1,:) aabb(1,1)], [aabb(2,:) aabb(2,1)], 'color', 'b');
            line([aabb2(1,:) aabb2(1,1)], [aabb2(2,:) aabb2(2,1)], 'color', 'y');
            hold off;

            % show normalized frame
            subplot(2, 2, 2)
            imshow(Iframe, []); hold on;
            for k = 1 : M
                plot(Points(1, k), Points(2, k), 'r.', 'LineWidth', 2, 'MarkerSize', 10);
                bb = makeAABB([Points(1, k) - options.components(1, k)/2 Points(2, k) - options.components(2, k)/2 ...
                               Points(1, k) + options.components(1, k)/2 Points(2, k) + options.components(2, k)/2 ]);
                line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', 'y');
            end;
            hold off;

            % show recomputed S also
            if (max(S(1, :)) < inf && min(S(3, :)) > -inf)
                subplot(2, 2, 3)
                imshow(Iframe, []); hold on;
                for k = 1 : M
                    plot(Points(1, k), Points(2, k), 'r.', 'LineWidth', 2, 'MarkerSize', 10);
                    bb = makeAABB(S(:, k));
                    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', 'b');
                    bb = makeAABB([Points(1, k) - options.components(1, k)/2 Points(2, k) - options.components(2, k)/2 ...
                                   Points(1, k) + options.components(1, k)/2 Points(2, k) + options.components(2, k)/2 ]);
                    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', 'y');
                end;
                hold off;
            end;

%             saveas(gcf, ['./img/' 'bw_' Annotation.image.filename]);
%             close gcf;
        end;
    end;
end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
classdef PowerNetChargeRateTermCondition < AbstractEventTerminationCondition
    %PowerNetChargeRateTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        netChargeRate(1,1) double = 0; %km
        initialStateLogEntry LaunchVehicleStateLogEntry
    end

    methods
        function obj = PowerNetChargeRateTermCondition(netChargeRate)
            obj.netChargeRate = netChargeRate;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.initialStateLogEntry = initialStateLogEntry;
        end

        function name = getName(obj)
            name = sprintf('Power Net Charge Rate (%.3f EC)', obj.netChargeRate);
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = true;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Net Charge Rate';
            params.paramUnit = 'EC/s';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = obj.netChargeRate;
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = NetChargeRateOptimizationVariable(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = PowerNetChargeRateTermCondition(paramValue);
        end
    end

    methods(Access=private)
        function [value,isterminal,direction] = eventTermCond(obj, t,y)
            initStateLogEntry = obj.initialStateLogEntry;

            powerStorageStates = initStateLogEntry.getAllActivePwrStorageStates();

            numTankStates = initStateLogEntry.getNumActiveTankStates();
            numPwrStorageStates = initStateLogEntry.getNumActivePwrStorageStates();
            [~, ~, ~, ~, storageSoCs] = AbstractPropagator.decomposeIntegratorTandY(t,y, numTankStates, numPwrStorageStates);

            stgStates = initStateLogEntry.stageStates;

            ut = initStateLogEntry.time;
            rVect = initStateLogEntry.position(:);
            vVect = initStateLogEntry.velocity(:);
            bodyInfo = initStateLogEntry.centralBody;

            steeringModel = initStateLogEntry.steeringModel;

            storageRates = LaunchVehicleStateLogEntry.getStorageChargeRatesDueToSourcesSinks(storageSoCs, powerStorageStates, stgStates, ut, rVect, vVect, bodyInfo, steeringModel);
            actualNetStorageRate = sum(storageRates);

            value = actualNetStorageRate - obj.netChargeRate;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [max,maxrest] = getMaxInARow(stimList,varargin)
	if nargin>1,dojitter = varargin{1},else dojitter = 0;,end

	max = 1; counter = 1;maxrest = 1;counterrest = 1;
        if dojitter, restnum = max(stimList);,else restnum = Inf;,end
	for i = 2:size(stimList,1)
		if stimList(i,1) == stimList(i-1,1)& ~stimList(i,1) == 0,
			counter = counter + 1;
			if counter > max, max = counter;,end
		else counter = 1;
                end

                if dojitter
                     if stimList(i,1) == restnum & stimList(i-1,1) == restnum
                        counterrest = counterrest+1;
                        if counterrest > maxrest,maxrest = counterrest;,end
                     else counterrest = 1;
		     end
                end

	end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function W = unfold( U, dir )
    %UNFOLD Left/right-unfold a 3D array.
    %   W = UNFOLD(U, DIR) unfolds the 3-dim. tensor U in direction DIR, where
    %   DIR is either 'LEFT' or 'RIGHT' (case insensitive).
    %
    %   See also MATRICIZE, TENSORIZE, TENSORPROD_TTEMPS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    d = size(U);
    % pad with 1 for the last dim (annoying)
    if length(size(U)) == 2
        d = [d, 1];
    end

    if strcmpi(dir, 'left')
        W = reshape( U, [d(1)*d(2), d(3)] );
    elseif strcmpi(dir, 'right')
        W = reshape( U, [d(1), d(2)*d(3)] );
    else
        error('Unknown direction specified. Choose either LEFT or RIGHT')
    end
end

author: NicolasBoumal

repo: manopt


text:
function F = transpose(F)
% .'   Transpose of a SPHEREFUNV.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

F.isTransposed = ~F.isTransposed;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [fx,J,dfdp] = f_embed0(Xt,Theta,ut,in)
% evolution function for dynamical system's delay embedding
% function [fx,J,dfdp] = f_embed(Xt,Theta,ut,in)

% First form delayed state vector and input
iX = (1:in.dim.n) + (in.options.delays(:)'.*in.dim.n);
dX = Xt(iX,:);

% Evaluate evolution function at the delayed state vector:
[opt,dim] = getOptions4EvalFun(in);
[fx0,J0,dfdp0] = VBA_evalFun('f',dX,Theta,ut,opt,dim);

% Construct full embedding flow and gradients:
Xe = Xt(1:in.dim.n_embed);
fx = [fx0;Xe];
J = zeros(in.dim.n_embed+in.dim.n);
J(iX,1:in.dim.n) = J0;
J(1:in.dim.n_embed,in.dim.n+1:in.dim.n+in.dim.n_embed) = ...
    eye(in.dim.n_embed);
dfdp = [dfdp0,zeros(in.dim.n_theta,in.dim.n_embed)];

function [opt,dim] = getOptions4EvalFun(in)
opt = in.options;
opt.f_fname = in.f_fname;
opt.f_nout = in.f_nout;
opt.checkGrads = 0;
dim.n = in.dim.n;
dim.n_theta = in.dim.n_theta;

author: MBB-team

repo: VBA-toolbox


text:
function [sFile, ChannelMat] = in_fopen_itab(DataFile)
% IN_FOPEN_ITAB: Open a ITAB raw MEG file.
%
% USAGE:  [sFile, ChannelMat] = in_fopen_itab(DataFile)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2018

%% ===== READ HEADER =====
% Locate header file
HeaderFile = [DataFile '.mhd'];
if ~file_exist(HeaderFile)
    error(['Header file was not found: ' HeaderFile]);
end
% Read header
hdr = read_itab_mhd(HeaderFile);
% Get endianness
switch (hdr.data_type)
    case {0,1,2},  byteorder = 'b';
    case {3,4,5},  byteorder = 'l';
    otherwise,     error('Data type not supported.');
end

%% ===== CREATE CHANNEL FILE =====
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'ITAB channels';
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.nchan]);
iChannels = 1:hdr.nchan;
% For each channel
for i = 1:hdr.nchan
    ChannelMat.Channel(i).Name    = hdr.ch(iChannels(i)).label;
    ChannelMat.Channel(i).Comment = '';
    % Position
    for iCoil = 1:hdr.ch(iChannels(i)).ncoils
        ChannelMat.Channel(i).Loc(:,iCoil)    = hdr.ch(iChannels(i)).position(iCoil).r_s' ./ 1000;
        ChannelMat.Channel(i).Orient(:,iCoil) = hdr.ch(iChannels(i)).position(iCoil).u_s' ./ 1000;
        ChannelMat.Channel(i).Weight(1,iCoil) = hdr.ch(iChannels(i)).wgt(iCoil);
    end
    % Type
    switch (hdr.ch(iChannels(i)).type)
        case 1,    ChannelMat.Channel(i).Type = 'EEG';
        case 2,    ChannelMat.Channel(i).Type = 'MEG';
        case 4,    ChannelMat.Channel(i).Type = 'EEG';
        case 8,    ChannelMat.Channel(i).Type = 'MEG REF';
        case 16,   ChannelMat.Channel(i).Type = 'Misc';     %AUX
        case 32,   ChannelMat.Channel(i).Type = 'Misc';     %PARAM
        case 64,   ChannelMat.Channel(i).Type = 'Misc';     %DIGIT
        case 128,  ChannelMat.Channel(i).Type = 'Misc';     %FLAG
        otherwise, ChannelMat.Channel(i).Type = 'Misc';     %Misc
    end

    % Type: everything below the "_"
%     iUnder = find(ChannelMat.Channel(i).Name == '_');
%     if (length(iUnder) == 1) && (iUnder > 1)
%         Type = ChannelMat.Channel(i).Name(1:iUnder-1);
%         % Rename some known types
%         switch (Type)
%             case 'MAG',   ChannelMat.Channel(i).Type = 'MEG';
%             case 'REF',   ChannelMat.Channel(i).Type = 'MEG REF';
%             case 'ELE',   ChannelMat.Channel(i).Type = 'EEG';
%             otherwise,    ChannelMat.Channel(i).Type = Type;
%         end
%     else
%         ChannelMat.Channel(i).Type = 'Misc';
%     end
end

% Channel flag
ChannelFlag = double([hdr.ch(iChannels).flag] == 0);
ChannelFlag(ChannelFlag == 0) = -1;

% %% ===== ADD MEG POSITIONS =====
% iChan = channel_find(ChannelMat.Channel, 'MEG');
% % Position / orientation
% if ~isempty(SensorDat)
%     % Split the matrix in sensor and
%     LocAll    = SensorDat(1:size(SensorDat,1)/2, :);
%     OrientAll = SensorDat((size(SensorDat,1)/2+1):end, :);
%     % Normalize orientation vector
%     OrientAll = bst_bsxfun(@rdivide, OrientAll, sqrt(sum(OrientAll.^2,2)));
%     % Apply to every sensor
%     for i = 1:length(iChan)
%         ind = str2double(ChannelMat.Channel(i).Name);
%         if ~isempty(ind) && ~isnan(ind) && (ind <= size(LocAll,1))
%             ChannelMat.Channel(i).Loc    = LocAll(ind,:)' ./ 1000;
%             ChannelMat.Channel(i).Orient = OrientAll(ind,:)';
%         end
%     end
% end
% % Add definition of sensors
% ChannelMat.Channel = ctf_add_coil_defs(ChannelMat.Channel, 'KRISS');

%% ===== EXTRA HEAD POINTS =====
% Get extra head points
HpLoc = hdr.marker;
iGood = find(~all(HpLoc == 0, 1));
HpLoc = HpLoc(:,iGood) ./ 1000;
nPoints = size(HpLoc,2);
% Save points in channel structure
if ~isempty(HpLoc) && (nPoints > 3)
    % All positions
    ChannelMat.HeadPoints.Loc = HpLoc;
    % All types and labels
    ChannelMat.HeadPoints.Type  = repmat({'EXTRA'}, 1, nPoints);
    ChannelMat.HeadPoints.Label = repmat({'EXTRA'}, 1, nPoints);
    % First three points: NAS/LPA/RPA
    ChannelMat.HeadPoints.Label{1} = 'NAS';
    ChannelMat.HeadPoints.Label{2} = 'RPA';
    ChannelMat.HeadPoints.Label{3} = 'LPA';
    ChannelMat.HeadPoints.Type{1} = 'CARDINAL';
    ChannelMat.HeadPoints.Type{2} = 'CARDINAL';
    ChannelMat.HeadPoints.Type{3} = 'CARDINAL';
    % Force re-alignment on the new set of NAS/LPA/RPA (switch from coil-based to SCS anatomical-based coordinate system)
    ChannelMat = channel_detect_type(ChannelMat, 1, 0);
end

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder  = byteorder;
sFile.filename   = DataFile;
sFile.format     = 'ITAB';
sFile.device     = 'ITAB MEG';
sFile.header     = hdr;
% Comment: short filename
[tmp__, sFile.comment, tmp__] = bst_fileparts(DataFile);
% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq   = hdr.smpfq;
sFile.prop.times   = [0, hdr.ntpdata - 1] ./ sFile.prop.sfreq;
sFile.prop.nAvg    = 1;
% No info on bad channels
sFile.channelflag = ChannelFlag;
% Acquisition date
sFile.acq_date = str_date(hdr.date);

%% ===== EVENTS =====
if (hdr.nsmpl >= 1)
    % Get events
    allStart = [hdr.smpl(1:hdr.nsmpl).start];
    allType  = [hdr.smpl(1:hdr.nsmpl).type];
    uniqueType = unique(allType);
    % Initialize list of events
    events = repmat(db_template('event'), 1, length(uniqueType));
    % Format list
    for iEvt = 1:length(uniqueType)
        % Find list of occurences of this event
        iOcc = find(strcmpi(allType, uniqueType(iEvt)));
        % Fill events structure
        events(iEvt).label      = num2str(uniqueType(iEvt));
        events(iEvt).color      = [];
        events(iEvt).reactTimes = [];
        events(iEvt).select     = 1;
        events(iEvt).times      = allStart(iOcc) ./ sFile.prop.sfreq;  % Get time and samples  (considering that samples are zero-based)
        events(iEvt).epochs     = ones(1, length(events(iEvt).times));  % Epoch: set as 1 for all the occurrences
        events(iEvt).channels   = [];
        events(iEvt).notes      = [];
    end
    % Import this list
    sFile = import_events(sFile, [], events);
end

author: brainstorm-tools


text:
function [planC, leksellInfo] = importLeksellPlan(planDir)
%"importLeksellPlan"
%   Imports a Leksell gamma knife plan into CERR format.  The directory
%   structure of the Leksell plan must be intact.
%
%By JRA 06/24/05
%
%LM: KRK, 05/31/07, all readLeksell<...>File functions are now complete and
%                   included in the importLeksellPlan function
%    KRK, 06/01/07, corrected the transform matrices of the scans
%    KRK, 06/04/07, added missing code in the structure import section to
%                   account for non associated scans
%    KRK, 06/05/07, added a second output variable (leksellInfo) to hold all
%                   information from the plan directory for easy access into
%                   the variables not currently used in CERR (metadata)
%    KRK, 06/07/07, the code for converting the normalized percentages in
%                   the dose matrices to units of Gy is now complete
%    KRK, 06/08/07, rescaled the imported scan images' intensities (MR only)
%                   so that they can now be viewed with the CERR CT Window
%                   preset "Head"
%    DK , 07/01/08, Rewrote the way directories are read. Change in file
%                   names and missing directories causing errors. Typo
%                   errors in converting mm to cm. Also typo error in
%                   reading CT Loop.
%
%Usage:
%   planC = importLeksellPlan(planDir)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

LekC = [];

files = dir(planDir);
lekSellDirStr = {files.name};

%%%%%%%%%%%%%%%%% Importing directories %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1: length(lekSellDirStr)
    if files(i).isdir
        if strfind(upper(lekSellDirStr{i}), 'MR')
            LekC = readLeksellDir(planDir,'MR', LekC, lekSellDirStr{i});
        elseif strfind(lekSellDirStr{i},'CT')
            LekC = readLeksellDir(planDir,'CT', LekC, lekSellDirStr{i});
        else
            %             try
            LekC = readLeksellDir(planDir, lekSellDirStr{i}, LekC, lekSellDirStr{i});
            %             catch
            %                 disp(['No Import Options for ' lekSellDirStr{i} ' ...'])
            %                 continue
            %             end
        end
    end
end

%--------Initialize planC-------------------------------------------------%
planC   = initializeCERR;
indexS  = planC{end};
planC{indexS.CERROptions} = CERROptions;

%--------Initialize scan information--------------------------------------%
if isfield(LekC, 'mrScanS')
    for i=1:length(LekC.mrScanS)
        planC{indexS.scan}(i).scanType = LekC.mrScanS(i).scanInfo(1).scanType;
        bool = ismember({LekC.studyS.modality}, LekC.mrScanS(i).scanInfo(1).scanType);
        matches = find(bool);
        if isempty(matches)
            planC{indexS.scan}(i).transM = [];
        else
            planC{indexS.scan}(i).transM = LekC.studyS(matches(1)).rcsToxyzTransM;
        end

        planC{indexS.scan}(i).scanArray = LekC.mrScanS(i).scanArray;
        planC{indexS.scan}(i).scanInfo = LekC.mrScanS(i).scanInfo;

        planC{indexS.scan}(i).scanUID = createUID('scan');

        scanTypeC = {planC{indexS.scan}(i).scanInfo(:).scanType};
        matchC = strfind(upper(scanTypeC), 'COR');
        toDelete = [];
        for indMatch = 1:length(matchC)
           if ~isempty(matchC{indMatch})  && matchC{indMatch} == 1
              toDelete = [toDelete indMatch];
           end
        end

        planC{indexS.scan}(i).scanInfo(toDelete) = [];

        planC{indexS.scan}(i).scanArray(:,:,toDelete) = [];

        for indSI = 1:length(planC{indexS.scan}(i).scanInfo)
           planC{indexS.scan}(i).scanInfo(indSI).imageType = 'MR';
        end

    end
end

if isfield(LekC, 'ctScanS')

    for i = 1:length(LekC.ctScanS)
        indX = length(planC{indexS.scan}) + 1;
        planC{indexS.scan}(indX).scanType = LekC.ctScanS(i).scanInfo(1).scanType;

        bool = ismember({LekC.studyS.modality}, LekC.ctScanS(i).scanInfo(1).scanType);
        matches = find(bool);
        planC{indexS.scan}(indX).transM = LekC.studyS(matches(1)).rcsToxyzTransM;

        planC{indexS.scan}(indX).scanArray = LekC.ctScanS(i).scanArray;
        planC{indexS.scan}(indX).scanInfo = LekC.ctScanS(i).scanInfo;

        planC{indexS.scan}(indX).scanUID = createUID('scan');

        for indSI = 1:length(planC{indexS.scan}(indX).scanInfo)
           planC{indexS.scan}(indX).scanInfo(indSI).imageType = 'CT';
        end

    end
end

%Remove transM if importing only one scan
if length(planC{indexS.scan}) == 1
    planC{indexS.scan}(1).savedtransM.transM = planC{indexS.scan}(1).transM;
    planC{indexS.scan}(1).savedtransM.name = 'GammaKnifeImport';
    planC{indexS.scan}(1).transM = [];
end

%Set the uniformized scan data in the plan (planC{indexS.scan})
planC = setUniformizedData(planC);
%The original transformation matrices read in by the readLeksellStudy
%function contained scaling as well as transformation/rotation.  Since CERR
%is compatible with only transformation/rotation matrices and scales
%according to its 3 voxel sizes, the scaling needs to be removed from the
%current transformation matrices.  To do this, multiply the original transM
%by the inverse of the pure scaling matrix (multiplied by 10 since the
%voxel sizes are in cm and must be mm).
for i=1:length(planC{indexS.scan})
    if ~isempty(planC{indexS.scan}(i).transM)
        planC{indexS.scan}(i).transM = planC{indexS.scan}(i).transM / [planC{indexS.scan}(i).uniformScanInfo.grid1Units*10 0 0 0; 0 planC{indexS.scan}(i).uniformScanInfo.grid2Units*10 0 0; 0 0 planC{indexS.scan}(i).uniformScanInfo.sliceThickness*10 0; 0 0 0 1];
        planC{indexS.scan}(i).transM(3,4) = planC{indexS.scan}(i).transM(3,4) - planC{indexS.scan}(i).uniformScanInfo.firstZValue;
    end
end

%-----------------Load structure(s)---------------------------------------%
structNum = 1;
nonAssocCounter = 0;
for i=1:length(LekC.volS)
    %Check if the scan is either MR or CT, if not, disregard
    regTo = LekC.volS(i).registeredTo;
    bool = ismember({planC{indexS.scan}.scanType}, regTo);
    associatedScan = find(bool);
    if isempty(associatedScan)
        %Continue if not associated with any scan, eg Leksell Skull.
        nonAssocCounter = nonAssocCounter + 1; %keep a count of the non associated scans to index associated scans correctly
        continue;
    end

    %Check for contours done on fused images (slicenums are therefore out of
    %the range CT and MR images).  There may be a way to map these fusion
    %contours to MR or CT, possibly using the information from
    %LekC.fusionS, but it hasn't been completed yet.
    fusionContour = 0;
    for j=1:length(LekC.volS(i).contour)
        for k=1:length(planC{indexS.scan})
            if(LekC.volS(i).contour(j).sliceNum > length(planC{indexS.scan}(k).scanInfo))
                fusionContour = 1;
                break;
            end
        end
    end
    if(fusionContour ~= 0)
        nonAssocCounter = nonAssocCounter + 1;
        continue;
    end

    %After checking if its associated, put in the structure information
    newStructS = newCERRStructure(associatedScan, planC);
    planC{indexS.structures} = dissimilarInsert(planC{indexS.structures},newStructS,i - nonAssocCounter);
    planC{indexS.structures}(i - nonAssocCounter).structureName = LekC.volS(i).structName;
    contour = LekC.volS(i).contour;
    for j=1:length(contour)
        sliceNum = contour(j).sliceNum + 1; %+ 1 to match Matlab's indexing (1-n+1 rather than 0-n)
        %scale the xy points based on the voxel sizes
        xyPoints = contour(j).contour;
        xyPoints(:,1) = planC{indexS.scan}(associatedScan).uniformScanInfo.grid1Units * contour(j).contour(:,1);
        xyPoints(:,2) = planC{indexS.scan}(associatedScan).uniformScanInfo.grid2Units * contour(j).contour(:,2);
        %set z values for the current x,y vertices
        xyPoints = [xyPoints ones(size(xyPoints(:,1)))*planC{indexS.scan}(associatedScan).scanInfo(sliceNum).zValue];
        %store the xyz points of each vertex
        planC{indexS.structures}(i - nonAssocCounter).contour(sliceNum).segments.points = xyPoints;
    end
    structNum = structNum + 1; %keep a count of all of the associated structures
    planC{indexS.structures}(i - nonAssocCounter).strUID = createUID('structure');
    planC{indexS.structures}(i - nonAssocCounter).assocScanUID = planC{indexS.scan}(associatedScan).scanUID;

    %Remove the contours from the second return variable to save space.
    %Only do this if the contour is read into CERR. In other words, don't
    %remove the contours from this return variable if they are from a non
    %associated scan (any contour not done on MR or CT, ex: fused).
    LekC.volS(i).contour = [];
end
%Calculate how to display the contoured segments inside of CERR
planC = getRasterSegs(planC);

%------------Load dose information----------------------------------------%
normalizingConst = -1; %used to find the global maximum normalizing constant in the dose matrix (see readLeksellDoseGroupFile for more information)
%Put dose info into the CERR data structure.

if isfield(LekC,'shotS')
    for i=1:length(LekC.shotS.doseMatrixInfo)
        dMI = LekC.shotS.doseMatrixInfo(i);

        planC{indexS.dose}(i).fractionGroupID = dMI.name;
        planC{indexS.dose}(i).horizontalGridInterval = dMI.gridSize;
        planC{indexS.dose}(i).verticalGridInterval = dMI.gridSize;
        planC{indexS.dose}(i).depthGridInterval = dMI.gridSize;

        try
            planC{indexS.dose}(i).doseArray = LekC.doses{i};
        catch
            planC{indexS.dose}(i) = [];
            continue
        end

        % Temporarily hold the max of each dose matrix to see if it is larger
        % than the other dose matrices' max value. If it is, store the new
        % larger value in normalizingConst (this only affects plans with multiple DMs).
        tmpNorm = max(planC{indexS.dose}(i).doseArray(:));
        if(tmpNorm > normalizingConst)
            normalizingConst = tmpNorm;
        end

        planC{indexS.dose}(i).coord1OFFirstPoint = LekC.shotS.doseMatrixInfo(i).minExtent(1);
        planC{indexS.dose}(i).coord2OFFirstPoint = LekC.shotS.doseMatrixInfo(i).minExtent(2);
        planC{indexS.dose}(i).coord3OfFirstPoint = LekC.shotS.doseMatrixInfo(i).minExtent(3);

        siz = size(planC{indexS.dose}(i).doseArray);
        planC{indexS.dose}(i).sizeOfDimension1 = siz(2);
        planC{indexS.dose}(i).sizeOfDimension2 = siz(1);
        planC{indexS.dose}(i).sizeOfDimension3 = siz(3);

        planC{indexS.dose}(i).zValues = linspace(LekC.shotS.doseMatrixInfo(i).minExtent(3), LekC.shotS.doseMatrixInfo(i).maxExtent(3), siz(3));

        planC{indexS.dose}(i).doseUID = createUID('dose');
        planC{indexS.dose}(i).doseUnits = 'Gy';
        planC{indexS.dose}(i).assocScanUID = [];
    end
end
%Normalize the dose arrays and convert them to Gy using previously found max.
%This converts all dose matrix values into CERR's native radiation units, Gy.
for i=1:length(planC{indexS.dose})
    planC{indexS.dose}(i).doseArray = (LekC.superS.globalMaxDose)*(planC{indexS.dose}(i).doseArray)/normalizingConst;
end

%Log end of import

%-------------------Store metadata----------------------------------------%
%Put all information obtained from the readLeksell<...>File functions into
%the other return variable to view and/or possibly incorporate into CERR.
%Also, get rid of some of the larger pieces of data (doses, pixel arrays)
%that have already been read into the CERR variable to save space.  See the
%"Load Structures" section of this function to see how the contours are
%removed from this variable.

LekC.doses = [];
try
    for i=1:length(LekC.ctScanS)
        LekC.ctScanS(i).scanArray = [];
    end
end
try
    for i=1:length(LekC.mrScanS)
        LekC.mrScanS(i).scanArray = [];
    end
end
leksellInfo = LekC;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% READ THE DIRECTORY STRUCTURE OF LEKSELL %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LekC = readLeksellDir(planDir,dirName, LekC, mrctdir)

switch dirName

    case 'Fusions'
        %Get Fusion info. (not used in CERR currently)
        disp('Reading Fusion info...')
        % fusionsFilename = fullfile(planDir, 'Fusions', 'Fusion.1');
        fusionsFilename = getLeksellFilesNames(fullfile(planDir, 'Fusions'), 'Fusion');
        for i = 1:length(fusionsFilename)
            LekC(i).fusionS = readLeksellFusionFile(fullfile(planDir, 'Fusions', fusionsFilename{i}));
        end

    case 'Orientations'
        %Get Orientations info. (not used in CERR currently)
        disp('Reading Orientations info...')
        % oriFilename = fullfile(planDir, 'Orientations', 'Orientation.1');
        oriFilename = getLeksellFilesNames(fullfile(planDir, 'Orientations'), 'Orientation');
        for i = 1:length(oriFilename)
            LekC(i).oriS = readLeksellOrientationFile(fullfile(planDir, 'Orientations', oriFilename{i}));
        end

    case 'SRdata'
        %Get Dosegroup data.
        disp('Reading Dosegroup info...')
        % doseFilename = fullfile(planDir, 'SRdata', ['DoseGroup' num2str(doseNum) '.1']);

        doseFilename = getLeksellFilesNames(fullfile(planDir, 'SRdata'), 'DoseGroup');

        for i = 1:length(doseFilename)
            LekC.doses{i} = readLeksellDoseGroupFile(fullfile(planDir, 'SRdata', doseFilename{i}));
        end

        %Get skull info. (not used in CERR currently)
        disp('Reading Skull info...')
        % skullFilename = fullfile(planDir, 'SRdata', 'Skull.1');
        skullFilename = getLeksellFilesNames(fullfile(planDir, 'SRdata'), 'Skull');
        for i = 1:length(skullFilename)
            LekC(i).skullS = readLeksellSkullFile(fullfile(planDir, 'SRdata', skullFilename{i}));
        end

        %Get super info.
        disp('Reading Super info...')

        % superFilename = fullfile(planDir, 'SRdata', 'Super.1');
        superFilename = getLeksellFilesNames(fullfile(planDir, 'SRdata'), 'Super');
        for i = 1:length(superFilename )
            LekC(i).superS = readLeksellSuperFile(fullfile(planDir, 'SRdata', superFilename{i}));
        end

        %Get target frame info. (not used in CERR currently)
        disp('Reading TargetFrame info...')
        % targetframeFilename = fullfile(planDir, 'SRdata', 'TargetFrame.1');
        targetframeFilename = getLeksellFilesNames(fullfile(planDir, 'SRdata'), 'TargetFrame');
        for i = 1:length(targetframeFilename)
            LekC(i).targetframeS = readLeksellTargetframeFile(fullfile(planDir, 'SRdata', targetframeFilename{i}));
        end

    case 'Shots'
        %Get the shot info.
        disp('Reading Shot info...')

        % shotFilename = fullfile(planDir, 'Shots', 'Shot.1');
        shotFilename = getLeksellFilesNames(fullfile(planDir, 'Shots'), 'Shot');

        if isempty(shotFilename)
            shotFilename = getLeksellFilesNames(fullfile(planDir, 'Shots'), 'Group');
        end

        for i = 1:length(shotFilename)
            LekC(i).shotS = readLeksellShotsFile(fullfile(planDir, 'Shots', shotFilename{i}));
        end

    case 'Studyregs'
        %Get the studyReg info. (not used in CERR currently)
        disp('Reading StudyReg info...')
        % srFilename = fullfile(planDir, 'Studyregs', 'Studyreg.1');
        srFilename = getLeksellFilesNames(fullfile(planDir, 'Studyregs'), 'studyreg');
        for i = 1:length(srFilename)
            LekC(i).studyRegS = readLeksellStudyregFile(fullfile(planDir, 'Studyregs', srFilename{i}));
        end

    case 'Studys'
        %Get the study info.
        disp('Reading Study info...')
        % studyFilename = fullfile(planDir, 'Studys', 'Study.1');
        studyFilename = getLeksellFilesNames(fullfile(planDir, 'Studys'), 'Study');
        for i = 1:length(studyFilename)
            LekC(i).studyS = readLeksellStudyFile(fullfile(planDir, 'Studys', studyFilename{i}));
        end

    case 'Volumes'
        %Get volume info.
        disp('Reading Volume info...')
        % volFilename = fullfile(planDir, 'Volumes', 'Volume.1');
        volFilename  = getLeksellFilesNames(fullfile(planDir, 'Volumes'), 'Volume');
        for i = 1:length(volFilename)
            LekC(i).volS = readLeksellVolumeFile(fullfile(planDir, 'Volumes', volFilename{i}));
        end

    case 'Plugs'
        %Get plug info. (not used in CERR currently)
        disp('Reading Plug info...')

        % plugFilename = fullfile(planDir, 'Plugs', 'Plug.1');
        plugFilename = getLeksellFilesNames(fullfile(planDir, 'Plugs'), 'Plug');
        for i = 1:length(plugFilename)
            LekC(i).plugS = readLeksellPlugFile(fullfile(planDir, 'Plugs', plugFilename{i}));
        end

    case 'Sortings'
        %Get sorting info. (not used in CERR currently)
        disp('Reading Sorting info...')
        % sortingFilename = fullfile(planDir, 'Sortings', 'Sorting.1');
        sortingFilename = getLeksellFilesNames(fullfile(planDir, 'Sortings'), 'Sorting');
        for i = 1:length(sortingFilename)
            LekC(i).sortingS = readLeksellSortingFile(fullfile(planDir, 'Sortings', sortingFilename{i}));
        end

    case 'Thresholds'
        %Get threshold info. (not used in CERR currently)
        disp('Reading Threshold info...')
        % thresholdFilename = fullfile(planDir, 'Thresholds', 'Threshold.1');
        thresholdFilename = getLeksellFilesNames(fullfile(planDir, 'Thresholds'), 'Threshold');
        for i = 1:length(thresholdFilename)
            LekC(i).thresholdS = readLeksellThresholdFile(fullfile(planDir, 'Thresholds', thresholdFilename{i}));
        end

    case 'MR'
        disp('Reading MR scan info...')

        if isfield(LekC,'mrScanS')
            indX = length(LekC.mrScanS) + 1;
        else
            indX = 1;
        end

        imageDir = fullfile(planDir, mrctdir);
        [mat, scanInfo] = readLeksellImageset(imageDir,'MR');
        if isempty(mat)
            return
        end
        LekC.mrScanS(indX).scanArray = mat;
        LekC.mrScanS(indX).scanInfo = scanInfo;
        %Put the pixel intensities in the range of 1000-1700. This code works, but
        %is not generally correct.  Modify the commented code under this to correct
        %the MR pixel intensity problem.
        %LekC.mrScanS(indX).scanArray = LekC.mrScanS(indX).scanArray - median(LekC.mrScanS(indX).scanArray(:)) + 1000;
        LekC.mrScanS(indX).scanArray = LekC.mrScanS(indX).scanArray; % - median(LekC.mrScanS(indX).scanArray(:)) + 1000;
        %Match the scan to the study, then subtract the lower of the ranges to
        %standardize the pixel intensity ranges (make it go from 0-722 rather than
        %its original values in the range of 32000-33000).
        %    for j=1:length(LekC.studyS)
        %        if(strmatch(LekC.studyS(j).modality, LekC.mrScanS(indX).scanInfo(1).scanType))
        %            LekC.mrScanS(indX).scanArray = LekC.mrScanS(indX).scanArray -
        %            LekC.studyS(j).pixelIntensityRange(1);
        %            break;
        %        end
        %    end

    case 'CT'
        if isfield(LekC,'ctScanS')
            indX = length(LekC.ctScanS) + 1;
        else
            indX = 1;
        end
        disp('Reading CT scan info...')

        imageDir = fullfile(planDir, mrctdir);

        [mat, scanInfo] = readLeksellImageset(imageDir,'CT');
        if isempty(mat)
            return
        end

        LekC.ctScanS(indX).scanInfo = scanInfo;

        LekC.ctScanS(indX).scanArray = mat;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% (this is just an interface)
function varargout = mydatehelp (varargin)
    [varargout{1:nargout}] = mydatestd_aux (varargin{:});
end

%!test
%! % just check whether it runs:
%! mydatehelp(1);

author: Sable

repo: mcbench-benchmarks


text:
% Simple Example #3 -- products of gaussian mixtures
%
%
fprintf('KDE Example #3 : Product sampling methods (single, anecdotal run)\n');
rand('state',0);
randn('state',0);

p = kde([.1,.45,.55,.8],.05);  % create a mixture of 4 gaussians for testing
q = kde([.25,.35,.7,.75],.05);  % and another
p = ksize(resample(p,100),'rot');
q = ksize(resample(q,150),'rot'); fprintf('.');
pd = hist(p,100,[1],[0,1]);
qd = hist(q,100,[1],[0,1]);

%pEx = productExact(p,{p,q},{},{}); fprintf('.');
pEx = pd.*qd; pEx = 100*pEx./sum(pEx);
dummy = kde(rand(1,200),1);  % placeholder for productApprox call
pSm = productApprox(dummy,{p,q},{},{},'exact'); fprintf('.');
pEp = productApprox(dummy,{p,q},{},{},'epsilon',1e-3); fprintf('.');
pGi = productApprox(dummy,{p,q},{},{},'gibbs1',100); fprintf('.');
pGM = productApprox(dummy,{p,q},{},{},'gibbsMS1',5); fprintf('.\n');

figure(1); hold off;
plot(linspace(0,1,100),pEx,'r-');
hold on;
plot(pSm,'k-');
plot(pEp,'g-');
plot(pGi,'c-');
plot(pGM,'b-');

legend('Discretized','Exact Sampling','Epsilon-Exact','Gibbs1','GibbsMS1');

author: ShapeNet

repo: RenderForCNN