text:
classdef AMICO_ACTIVEAX

properties
    id, name                % id and name of the model
    dPar                    % parallel diffusivity of the tensors [units of mm^2/s]
    dIso                    % isotropic diffusivity [units of mm^2/s]
    IC_Rs                   % radii of the axons [units of 1E-6 (micrometers)]
    IC_VFs                  % volume fractions of the axons
    OUTPUT_names            % suffix of the output maps
    OUTPUT_descriptions     % description of the output maps
end

methods

    % =================================
    % Setup the parameters of the model
    % =================================
    function obj = AMICO_ACTIVEAX()
        global CONFIG

        % set the parameters of the model
        obj.id        = 'ACTIVEAX';
        obj.name      = 'ActiveAx';
        obj.dPar      = 0.6 * 1E-3;
        obj.dIso      = 2.0 * 1E-3;
        obj.IC_Rs     = [0.01 linspace(0.5,10,20)];
        obj.IC_VFs    = [0.3:0.1:0.9];

        obj.OUTPUT_names        = { 'v', 'a', 'd' };
        obj.OUTPUT_descriptions = {'Intra-cellular volume fraction', 'Mean axonal diameter', 'Axonal density'};

        % set the parameters to fit it
        CONFIG.OPTIMIZATION.SPAMS_param.mode    = 2;
        CONFIG.OPTIMIZATION.SPAMS_param.pos     = true;
        CONFIG.OPTIMIZATION.SPAMS_param.lambda  = 0.25; % l1 regularization
        CONFIG.OPTIMIZATION.SPAMS_param.lambda2 = 4;    % l2 regularization
    end

    % ==================================================================
    % Generate high-resolution kernels and rotate them in harmonic space
    % ==================================================================
    function GenerateKernels( obj, ATOMS_path, schemeHR, AUX, idx_IN, idx_OUT )
        global CONFIG AMICO_data_path CAMINO_path

        % check if high-resolution scheme has been created
        schemeHrFilename = fullfile(ATOMS_path,'protocol_HR.scheme');
        if ~exist( schemeHrFilename, 'file' )
            error( '[AMICO_GenerateKernels_ACTIVEAX] File "protocol_HR.scheme" not found in folder "%s"', ATOMS_path )
        end

        filenameHr = [tempname '.Bfloat'];
        progress = ProgressBar( numel(obj.IC_Rs) + numel(obj.IC_VFs) + 1 );

        % Restricted
        % ==========
        for R = obj.IC_Rs
            % generate
            if exist( filenameHr, 'file' ), delete( filenameHr ); end
            CMD = sprintf( '%s/datasynth -synthmodel compartment 1 CYLINDERGPD %E 0 0 %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.dPar*1e-6, R*1e-6, schemeHrFilename, filenameHr );
            [status result] = system( CMD );
            if status>0
                disp(result)
                error( '[AMICO_ACTIVEAX.GenerateKernels] Problems generating the signal with datasynth' );
            end

            % rotate and save
            fid = fopen( filenameHr, 'r', 'b' );
            signal = fread(fid,'float');
            fclose(fid);
            delete( filenameHr );
            lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, false );
            save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )
            progress.update();
        end

        % Hindered
        % ========
        for ICVF = obj.IC_VFs
            % generate
            d_perp = obj.dPar * ( 1.0 - ICVF );
            if exist( filenameHr, 'file' ), delete( filenameHr ); end
            CMD = sprintf( '%s/datasynth -synthmodel compartment 1 ZEPPELIN %E 0 0 %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.dPar*1e-6, d_perp*1e-6, schemeHrFilename, filenameHr );
            [status result] = system( CMD );
            if status>0
                disp(result)
                error( '[AMICO_ACTIVEAX.GenerateKernels] problems generating the signal' );
            end

            % rotate and save
            fid = fopen( filenameHr, 'r', 'b' );
            signal = fread(fid,'float');
            fclose(fid);
            delete( filenameHr );
            lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, false );
            save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )

            progress.update();
        end

        % Isotropic
        % =========
        % generate
        if exist( filenameHr, 'file' ), delete( filenameHr ); end
        CMD = sprintf( '%s/datasynth -synthmodel compartment 1 BALL %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.dIso*1e-6, schemeHrFilename, filenameHr );
        [status result] = system( CMD );
        if status>0
            disp(result)
            error( '[AMICO_ACTIVEAX.GenerateKernels] problems generating the signal' );
        end

        % resample and save
        fid = fopen( filenameHr, 'r', 'b' );
        signal = fread(fid,'float');
        fclose(fid);
        delete( filenameHr );
        lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, true );
        save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )
        progress.update();

        progress.close();
    end

    % ==============================================
    % Project kernels from harmonic to subject space
    % ==============================================
    function ResampleKernels( obj, ATOMS_path, idx_OUT, Ylm_OUT )
        global CONFIG AMICO_data_path KERNELS

        % Setup the KERNELS structure
        % ===========================
        nIC = numel(obj.IC_Rs);
        nEC = numel(obj.IC_VFs);

        KERNELS = {};
        KERNELS.model    = 'ACTIVEAX';
        KERNELS.nS       = CONFIG.scheme.nS;
        KERNELS.nA       = nIC + nEC + 1; % number of atoms

        KERNELS.dPar     = obj.dPar;

        KERNELS.Aic      = zeros( [KERNELS.nS nIC 181 181], 'single' );
        KERNELS.Aic_R    = zeros( 1, nIC, 'single' );

        KERNELS.Aec      = zeros( [KERNELS.nS nEC 181 181], 'single' );
        KERNELS.Aec_icvf = zeros( 1, nEC, 'single' );

        KERNELS.Aiso     = zeros( [KERNELS.nS 1], 'single' );
        KERNELS.Aiso_d   = NaN;

        progress = ProgressBar( KERNELS.nA );

        % Restricted
        % ==========
        for i = 1:nIC
            load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
            KERNELS.Aic(:,i,:,:) = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, false );
            KERNELS.Aic_R(i)     = obj.IC_Rs(i);
            progress.update();
        end

        % Hindered
        % ========
        for i = 1:nEC
            load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
            KERNELS.Aec(:,i,:,:) = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, false );
            KERNELS.Aec_icvf(i)  = obj.IC_VFs(i);
            progress.update();
        end

        % Isotropic
        % =========
        load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
        KERNELS.Aiso   = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, true );
        KERNELS.Aiso_d = obj.dIso;
        progress.update();

        progress.close();
    end

    % ===========================
    % Fit the model to each voxel
    % ===========================
    function [ MAPs ] = Fit( obj, y, i1, i2 )
        global CONFIG KERNELS

        % prepare SIGNAL and DICTIONARY
        if numel(KERNELS.Aiso_d) > 0
            A = double( [ KERNELS.Aic(CONFIG.scheme.dwi_idx,:,i1,i2) KERNELS.Aec(CONFIG.scheme.dwi_idx,:,i1,i2) KERNELS.Aiso(CONFIG.scheme.dwi_idx) ] );
        else
            A = double( [ KERNELS.Aic(CONFIG.scheme.dwi_idx,:,i1,i2) KERNELS.Aec(CONFIG.scheme.dwi_idx,:,i1,i2) ] );
        end
        AA = [ ones(1,KERNELS.nA) ; A ];
        yy = [ 1 ; y(CONFIG.scheme.dwi_idx) ];

        % estimate coefficients
        x = full( mexLasso( yy, AA, CONFIG.OPTIMIZATION.SPAMS_param ) );

        % compute MAPS
        nIC = numel(obj.IC_Rs);
        nEC = numel(obj.IC_VFs);
        f1 = sum( x( 1:nIC ) );
        f2 = sum( x( (nIC+1):(nIC+nEC) ) );
        MAPs(1) = f1 / ( f1 + f2 + eps );                           % intra-cellular volume fraction (v)
        MAPs(2) = 2 * KERNELS.Aic_R * x(1:nIC) / ( f1 + eps );      % mean axonal diameter
        MAPs(3) = (4*MAPs(1)) / ( pi*MAPs(2)^2 + eps );             % axonal density
    end

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function script_rfcn_VOC0712_ResNet50_OHEM_ss()
% script_rfcn_VOC0712_ResNet50_OHEM_ss()
% RFCN training and testing with OHEM using ResNet50 model and selective
% search proposals
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

clc;
clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));
%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe_rfcn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);

% model
model                       = Model.ResNet50_for_RFCN_VOC0712_OHEM();
% cache name
opts.cache_name             = 'rfcn_VOC0712_ResNet50_OHEM_ss';
% config
conf                        = rfcn_config_ohem('image_means', model.mean_image);
% train/test data
fprintf('Loading dataset...')
dataset                     = [];
dataset                     = Dataset.voc0712_trainval_ss(dataset, 'train', conf.use_flipped);
dataset                     = Dataset.voc2007_test_ss(dataset, 'test', false);
fprintf('Done.\n');

% do validation, or not
opts.do_val                 = true;

%% -------------------- TRAINING --------------------

opts.rfcn_model        = rfcn_train(conf, dataset.imdb_train, dataset.roidb_train, ...
                                'do_val',           opts.do_val, ...
                                'imdb_val',         dataset.imdb_test, ...
                                'roidb_val',        dataset.roidb_test, ...
                                'solver_def_file',  model.solver_def_file, ...
                                'net_file',         model.net_file, ...
                                'cache_name',       opts.cache_name, ...
                                'caffe_version',    opts.caffe_version);
assert(exist(opts.rfcn_model, 'file') ~= 0, 'not found trained model');

%% -------------------- TESTING --------------------
                          rfcn_test(conf, dataset.imdb_test, dataset.roidb_test, ...
                                'net_def_file',     model.test_net_def_file, ...
                                'net_file',         opts.rfcn_model, ...
                                'cache_name',       opts.cache_name,...
                                'ignore_cache',     true);

end

repo: R-FCN


text:
%% COLLIDINGFLOW colliding flow in a square domain
%
% Stokes equations on the square [-1,1]^2.
%
% Reference: page 237 in Finite Elements and Fast Iterative Solvers with
% Applications in Incompressible Fluid Dynamics. by Howard C. Elman, David
% J. Silvester, and Andrew J. Wathen.
%
% See also squareStokes, Poiseuilieflow, StokesP2P1
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear all;
%% Parameters
maxIt = 3; N = zeros(maxIt,1);
erru = zeros(maxIt,1); errp = zeros(maxIt,1);

%% Generate an initial mesh
[node,elem] = squaremesh([-1 1 -1 1], 0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

%% PDE and options
pde = Stokesdata1;

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [u,p,edge,A] = StokesP2P0(node,elem,pde,bdFlag);
    N(k) = length(u)+length(p);
    if N(k) < 2e3 % show solution for small size
        figure(1);  showresult(node,elem,p);
    end
    % compute error
    uI = pde.exactu([node; (node(edge(:,1),:)+node(edge(:,2),:))/2]);
    erru(k) = sqrt((u-uI(:))'*A*(u-uI(:)));
    errp(k) = getL2error(node,elem,pde.exactp,p);
    % refine mesh
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
end

%% Plot convergence rates
figure(2); clf;
showrate2(N,erru,1,'-*','||Du_I-Du_h||',N,errp,1,'k-+','|| p - p_h||');

repo: Aero-matlab


text:
function [scimat, sigma] = scimat_estimate_bias_field(scimat, x, a, sigma, FASTINTERP)
% SCIMAT_ESTIMATE_BIAS_FIELD  Estimate MRI bias field.
%
%   This function provides an estimate of the bias field from a magnetic
%   resonance image (MRI).
%
%   The function samples a SCIMAT image in locations provided by the
%   user (typically, corresponding to background voxels), and then creates
%   another image that interpolates the sampled intensity values.
%
%   The locations can be selected e.g. using our Spline Tool extension to
%   the Seg3D platform (click the points, and then export the control
%   points).
%
%   Note that the interpolant is an interpolating thin-plate spline (TPS),
%   but instead of just sampling the selected voxels, a neighbourhood is
%   sampled, and used to locally low-pass filter the image with a Gaussian
%   filter. Thus, this function is somehow robust in the presence of noise.
%   Alternatively, this function could use an approximating TPS instead of
%   an interpolating one.
%
% SCIMAT2 = scimat_estimate_bias_field(SCIMAT, X)
%
%   SCIMAT is an image (see "help scimat" for details).
%
%   X is a 3-colum matrix with the real world coordinates of the sampling
%   points. Note that each point is rounded to the closest voxel centre.
%
%   SCIMAT2 is the estimated bias field in SCIMAT format.
%
% SCIMAT2 = scimat_estimate_bias_field(SCIMAT, X, A, SIGMA, FASTINTERP)
%
%   A is a scaling factor. Because using the TPS to interpolate all voxels
%   can be rather slow, and the bias field is anyway a slow varying field,
%   it's convenient to first quickly reduce the image size by A (using
%   bilinear interpolation), interpolate the bias field in the smaller
%   image with the TPS, and then expand to the original size. By default, A
%   = 1.0 and no rescaling is used. If A has 1 value, then the same scaling
%   is applied to all dimensions. If A has 3 values, each value is used to
%   scale one dimension.
%
%   Note: Small values of A imply a large SIGMA to avoid aliasing.
%   Internally, this attempts to create a large Gaussian filter, and the
%   computer may run out of memory. Larger values of A require a lot of
%   memory from the thin-plate spline, and this may give an out of memory
%   error too.
%
%   SIGMA is the standard deviation of the Gaussian filter. By default,
%   sigma = 2*sqrt(-2 ln(alpha)) / a, alpha=1/sqrt(2)*ones(1,3), so that
%   the Gaussian filter has a 3dB drop at the cut-off frequency in each
%   dimension (see Note 1 for details).
%
%   FASTINTERP is a boolean to decide whether fast (more memory consuming)
%   thin-plate spline interpolation should be used. By default,
%   FASTINTERP=true. If you run out of memory, try setting it to false.
%
% [SCIMAT2, SIGMA] = ...
%
%   SIGMA as an output is the sigma value (in voxel units) used for the
%   Gaussian filter.
%
% ======
% Note 1
% ======
%
% The Gaussian filter in the frequency domain is
%
%   F(u) = 1/sqrt(2*pi) * exp(-u^2 *sigma^2 / 2)
%
% So if we want our filter to have a drop alpha at cut-off frequency B'
%
%   exp(-B'^2 * sigma^2 / 2) = alpha
%
% then
%
%   sigma = sqrt(-2 ln(alpha)) / B'
%
% Because of the Nyquist theorem, in order to prevent aliasing, we require
% the cut-off frequency to be
%
%   B' < fs' / 2 = a / (2*dx)
%
% then the standard deviation for the Gaussian must be
%
%   sigma > 2*sqrt(-2 ln(alpha)) / fs
%
%   sigma > 2*sqrt(-2 ln(alpha)) * dx / a
%
% There are two free parameters:
%
%   * alpha: filter drop at cut-off frequency (alpha=1/sqrt(2) makes a 3 dB
%            cut-off point)
%
%   * a: scaling factor for the image (a=.5, image becomes half size)
%
%  and a parameter that depends on the image
%
%   * dx: voxel size, image resolution; because sigma is given in pixel
%     units, dx = 1 for any image.
%
% ======
% Note 2
% ======
%
% The Gaussian filter in the space domain is
%
%   f(x) = 1/(sqrt(2*pi)*sigma) * exp(-x^2 / (2*sigma^2))
%
% The Gaussian tails tend asymptotically to 0, so we need to truncate the
% filter at some point. For sigma = 1.0, the filter tail at
% f(4)=1.3383e-04, which is a good compromise between having a small
% numerical error and a small filter.

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

% check arguments
narginchk(2, 5);
nargoutchk(0, 2);

% defaults
alpha = 1/sqrt(2)*ones(1,3);
if (nargin < 3 || isempty(a))
    a = 1.0;
end
if (length(a) == 1)
    a = a([1 1 1]);
elseif (length(a) ~= 3)
    error('A must have 1 or 3 values')
end
if (nargin < 4 || isempty(sigma))
    sigma = 2*sqrt(-2*log(alpha)) ./ a;
end
if (nargin < 5 || isempty(FASTINTERP))
    FASTINTERP = true;
end

% squeeze volume
scimat = scimat_squeeze(scimat);

% size of input volume
nin = [scimat.axis.size];

% recompute the scaling factor so that we obtain an integer number of
% voxels in the output volume
nout = round(nin .* a);
a = nout ./ nin;

% Seg3D provides continuous point coordinates over the whole image volume.
% But in fact, intensity values correspond to the voxel centre, so we have
% to convert the coordinates to indices and round them
x = scimat_world2index(x, scimat);
x = round(x);

% compute low-pass anti-aliasing filter
% sigma = 2 * sqrt(-2 * log(alpha)) ./ a .* dx; % dx not needed, sigma
% given in pixel units
halfsz = round(4 * sigma);
h = fspecial3('gaussian', halfsz*2+1, sigma);

% loop each sampling point. For each sampling point we want to sample a
% whole cube around it, so that we can apply a low-pass filter without
% having to filter the whole image volume
v = zeros(size(x, 1), 1);
for I = 1:size(x, 1)

    % get indices of the area around the sampling point to sample. If the
    % sampling point is too close to the edge, we have to be careful to not
    % overflow
    idxr = max(1, x(I, 1)-halfsz(1)):min(nin(1), x(I, 1)+halfsz(1));
    idxc = max(1, x(I, 2)-halfsz(2)):min(nin(2), x(I, 2)+halfsz(2));
    idxs = max(1, x(I, 3)-halfsz(3)):min(nin(3), x(I, 3)+halfsz(3));

    % if necessary, crop low-pass filter so that it has the same size as
    % the sampling area
    hbox = h(...
        max(1, idxr - x(I, 1) + halfsz(1) + 1), ...
        max(1, idxc - x(I, 2) + halfsz(2) + 1), ...
        max(1, idxs - x(I, 3) + halfsz(3) + 1));

    % extract image area around the sampling point
    im = scimat.data(idxr, idxc, idxs);

    % compute low-pass filtered intensity of current voxel
    v(I) = sum(hbox(:) .* im(:)) / sum(hbox(:));

%     v(I) = sum(hbox(:) .* im(:)); % this is how Matlab's imfilter() does
%     it, which is incorrect near the edges because it reduces the value of
%     the pixels

end

% clear memory
clear h hbox im

% convert back to coordinates
x = scimat_index2world(x, scimat);

% compute coordinates of the two extreme voxels that define the image
% volume
cmin = scimat_index2world([1 1 1], scimat);
cmax = scimat_index2world(nin, scimat);

% if we have point coordinates with values like 100, 400, matrix L for the
% thin-plate spline weight computation is badly scaled. Thus, we make use
% of the invariability of thin-plate splines to scaling and make all the
% coordinate values <= 5.0 (if we make them <= 1.0, interpn gives an error
% saying that they grid is not monotonic)
K = max(cmax)/5;

% the warp will be defined as from the xyz-points to the corresponding
% intensity values

% compute what would be the size of the image if we downsample it
nmid = round(nin .* a);

% compute weights for thin-plate spline interpolation
w = pts_tps_weights(x/K, v);

% generate grid
[gx, gy, gz] = meshgrid(...
    linspace(cmin(1)/K, cmax(1)/K, nmid(2)), ...
    linspace(cmin(2)/K, cmax(2)/K, nmid(1)), ...
    linspace(cmin(3)/K, cmax(3)/K, nmid(3)));

% clear memory
scimat.data = [];

% interpolate intensity values for each point in the grid
scimat.data = pts_tps_map(x/K, v, [gx(:) gy(:) gz(:)], w, FASTINTERP, false);

% clear memory
clear gx gy gz

% reshape the interpolated values to go from a vector to an image volume
scimat.data = reshape(single(scimat.data), nmid);

% recover original size
scimat.data = tformarray(scimat.data, ...
    maketform('affine', [1/a(1) 0 0 0; 0 1/a(2) 0 0; 0 0 1/a(3) 0; 0 0 0 1]), ...
    makeresampler('linear', 'replicate'), ...
    1:length(scimat.axis), 1:length(scimat.axis), ...
    nin, [], []);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout = progress(varargin)

% This function is a backward compatibility wrapper. It allows existing
% MATLAB scripts that do not use the new FieldTrip ft_xxx function naming
% scheme to work with recent versions of the FieldTrip toolbox.
%
% Please look in ft_xxx for the help of the function that you are looking
% for, where xxx is the name of the function that you were looking for.

% Copyright (C) 2009, Robert Oostenveld

prefix    = 'ft_';
funname   = mfilename;
warning([upper(mfilename) ' is only a compatibility wrapper, which will soon be removed. Please instead call ' upper(prefix) upper(funname) '.']);
funhandle = str2func([prefix funname]);
[varargout{1:nargout}] = funhandle(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ stack, W_t ] = initialize_weights( eI )
%INITIALIZE_WEIGHTS Random weight structures for a network architecture
%   eI describes an RNN via the fields layerSizes, inputDim and
%   temporalLayer
%
%   This uses Xavier's weight initialization tricks for better backprop
%   See: X. Glorot, Y. Bengio. Understanding the difficulty of training
%        deep feedforward neural networks. AISTATS 2010.

%% initialize hidden layers
stack = cell(1, numel(eI.layerSizes));
for l = 1 : numel(eI.layerSizes)
    if l > 1
        prevSize = eI.layerSizes(l-1);
    else
        prevSize = eI.inputDim;
    end;
    curSize = eI.layerSizes(l);
    % Xaxier's scaling factor
    s = sqrt(6) / sqrt(prevSize + curSize);
    % Ilya suggests smaller scaling for recurrent layer
    if l == eI.temporalLayer
        s = sqrt(6) / sqrt(prevSize + 2*curSize);
    end;
    stack{l}.W = rand(curSize, prevSize)*2*s - s;
    stack{l}.b = zeros(curSize, 1);
end
%% weight tying
% default weight tying to false
if ~isfield(eI, 'tieWeights')
    eI.tieWeights = 0;
end;
% overwrite decoder layers for tied weights
if eI.tieWeights
    decList = [(numel(eI.layerSizes)/2)+1 : numel(eI.layerSizes)-1];
    for l = 1:numel(decList)
        lDec = decList(l);
        lEnc = decList(1) - l;
        assert( norm(size(stack{lEnc}.W') - size(stack{lDec}.W)) == 0, ...
            'Layersizes dont match for tied weights');
        stack{lDec}.W = stack{lEnc}.W';
    end;
end;
%% initialize temporal weights if they should exist
W_t = [];
if eI.temporalLayer
    % assuems temporal init type set
    if strcmpi(eI.temporalInit, 'zero')
        W_t = zeros(eI.layerSizes(eI.temporalLayer));
    elseif strcmpi(eI.temporalInit, 'rand')
        % Ilya's modification to Xavier's update rule
        s = sqrt(6) / sqrt(3*eI.layerSizes(eI.temporalLayer));
        W_t = rand(eI.layerSizes(eI.temporalLayer))*2*s - s;
    elseif strcmpi(eI.temporalInit, 'eye')
        W_t = eye(eI.layerSizes(eI.temporalLayer));
    else
        error('unrecognized temporal initialization: %s', eI.temporalInit);
    end;
end;
%% init short circuit connections
% default short circuits to false
if ~isfield(eI, 'shortCircuit')
    eI.shortCircuit = 0;
end;
if eI.shortCircuit
    %padSize = (eI.winSize-1) / 2;
    %stack{end}.W_ss = [zeros(eI.featDim, padSize*eI.featDim), eye(eI.featDim),...
    %    zeros(eI.featDim, padSize*eI.featDim)];
    % use random init since input might contain noise estimate
    s = sqrt(6) / sqrt(eI.inputDim + eI.layerSizes(end));
    stack{end}.W_ss = rand(eI.inputDim, eI.layerSizes(end))*2*s - s;
end;

author: amaas

repo: rnn-speech-denoising


text:
function test_triangulation_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests P00_TEST_NUM, P00_TITLE and P00_HEADER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  P00_TEST_NUM reports the number of problems.\n' );
  fprintf ( 1, '  P00_TITLE returns a title for each problem.\n' );
  fprintf ( 1, '  P00_HEADER prints some information about each problem.\n' );

  test_num = p00_test_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of tests available = %d\n', test_num );

  for test = 1 : test_num

    title = p00_title ( test );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Test number: %d\n', test );
    fprintf ( 1, '  Title:      "%s"\n', title );

    p00_header ( test );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = gpr(hyper, cov, x, y, xs)

% gpr - Gaussian process regression, with a named covariance function. Two
% modes are possible: training and prediction: if no test data are given, the
% function returns minus the log likelihood and its partial derivatives with
% respect to the hyperparameters; this mode is used to fit the hyperparameters.
% If test data are given, then (marginal) Gaussian predictions are computed,
% whose mean and variance are returned. Note that in cases where the covariance
% function has noise contributions, the variance returned in s2 is for noisy
% test targets; if you want the variance of the noise-free latent function, you
% must substract the noise variance.
%
% usage: [nlZ dnlZ] = gpr(hyp, cov, x, y)
%    or: [mu s2]    = gpr(hyp, cov, x, y, xs)
%
% where:
%
%   hyp      is a (column) vector of log hyperparameters
%   cov      is the covariance function
%   x        is a n by D matrix of training inputs
%   y        is a (column) vector (of size n) of targets
%   xs       is a ns by D matrix of test inputs
%   nlZ      is the returned value of the negative log marginal likelihood
%   dnlZ     is a (column) vector of partial derivatives of the negative
%                 log marginal likelihood wrt each log hyperparameter
%   mu       is a (column) vector (of size nn) of prediced means
%   s2       is a (column) vector (of size nn) of predicted variances
%
% For more help on covariance functions, see "help covFunctions".
%
% Copyright (c) 2010 Carl Edward Rasmussen and Hannes Nickisch 2010-06-18.

err = 'we need cov = {''covSum'', {cov1, ..,''covNoise'', .., covP} }; to map to gp.m';
if ~strcmp(cov{1},'covSum'), error(err), end
id = 0; nhyp = zeros(length(cov{2}),1);
for i=1:length(cov{2})
  if strcmp(cov{2}{i},'covNoise'), id = i; break, end
  nhyp(i) = eval(feval(cov{2}{i}));
end
if id==0, error(err), else nhyp = sum(nhyp(1:id-1)); end
cov{2} = cov{2}([1:id-1,id+1:end]);
hyp.lik = hyper(nhyp+1); hyp.cov = hyper([1:nhyp,nhyp+2:end]);

% Note, this function is just a wrapper provided for backward compatibility,
% the functionality is now provided by the more general gp function.
lik = @likGauss; inf = @infExact; mean = @meanZero;

varargout = cell(nargout, 1);    % allocate the right number of output arguments
if nargin==4
  [varargout{:}] = gp(hyp,inf,mean,cov,lik,x,y);
  if nargout>1, varargout{2} = [varargout{2}.lik; varargout{2}.cov]; end
else
  [varargout{:}] = gp(hyp,inf,mean,cov,lik,x,y,xs);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% get_zc.m
% Purpose: returns Zc for given alpha
% Change History:
%	021100 tdr created

function zc = get_zc(alpha)

zc=norminv(1-(alpha/2),0,1);

author: Sable

repo: mcbench-benchmarks


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function showSyms(obj)
for i = 1 : length(obj)
    assignin('caller', char(obj(i)), obj(i));
end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef TestCvtColor
    %TestCvtColor

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_1
            img = imread(TestCvtColor.im);
            [h,w,~] = size(img);
            result = cv.cvtColor(img, 'RGB2GRAY');
            validateattributes(result, {class(img)}, {'2d', 'size',[h,w]});
        end

        function test_2
            img = imread(TestCvtColor.im);
            result = cv.cvtColor(img, 'RGB2GRAY', 'DstCn',0);
        end

        function test_3
            img = cv.imread(TestCvtColor.im, 'Grayscale',true);
            result = cv.cvtColor(img, 'GRAY2RGB');
            validateattributes(result, {class(img)}, ...
                {'3d', 'size',[NaN NaN 3]});
        end

        function test_4
            codes = {'RGB2BGR', 'RGB2XYZ', 'RGB2YCrCb', 'RGB2YUV', ...
                'RGB2HSV', 'RGB2HLS', 'RGB2Lab', 'RGB2Luv'};
            img = imread(TestCvtColor.im);
            [h,w,~] = size(img);
            for i=1:numel(codes)
                result = cv.cvtColor(img, codes{i});
                validateattributes(result, {class(img)}, ...
                    {'3d', 'size',[h,w,3]});
            end
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [Population,W,B] = WeightUpdate(Population,W,Archive,Z,T,Global)
% Weight Update

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Routine to find undeveloped individuals (correspondingly their weights) in the archive set
    % Normalisation
    N_arc         = length(Archive);
    fmin_arc      = min(Archive.objs);
    fmax_arc      = max(Archive.objs);
    Archiveobjs   = (Archive.objs - repmat(fmin_arc,N_arc,1) )./repmat(fmax_arc - fmin_arc,N_arc,1);
    Populaionobjs = (Population.objs - repmat(fmin_arc,Global.N,1) )./repmat(fmax_arc - fmin_arc,Global.N,1);
    % Euclidean distance between individuals in the archive set and individuals in the Population
    dis1 = pdist2(Archiveobjs,Populaionobjs);
    dis1 = sort(dis1,2);
    % Euclidean distance between any two individuals in the archive set
    dis2 = pdist2(Archiveobjs,Archiveobjs);
    dis2 = sort(dis2,2);
    % Calculate the niche size(median of the distances from their closest solution in the archive )
    niche_size = median(dis2(:,2));
    % Find undeveloped
    Archive_und = Archive(dis1(:,1) >= niche_size);
    N_und = length(Archive_und);

    %% If the undeveloped individuals are promising then add them into the evolutionary Population
    % Obtain their corresponding weights.
	if ~isempty(Archive_und)
        W1 = (Archive_und.objs - repmat(Z,N_und,1))./repmat( sum(Archive_und.objs,2)-repmat(sum(Z),N_und,1), 1, Global.M );
        for i = 1 : size(W1,1)
            W_all = [W;W1(i,:)];
            B1 = pdist2(W_all,W_all);
            B1(logical(eye(length(B1)))) = inf;
            [~,B1] = sort(B1,2);
            B1 = B1(:,1:T);

            Population1 = [Population,Archive_und(i)];
            Population2 = Population1(B1(end,:));

            Value_Tche_all = max(abs(Population2.objs-repmat(Z,T,1))./repmat(W1(i,:),T,1),[],2);
            Value_Tche     = max(abs(Archive_und(i).obj -    Z     )./W1(i,:),[],2);
            index = find(Value_Tche_all<Value_Tche, 1);

            if isempty(index)
                % Put the wight into the W, as well as the corresponding solution
                W = [W;W1(i,:)];
                Population = [Population Archive_und(i)];

                % Update neighbour solutions after adding a weight
                P = B1(end,:);
                g_old = max( abs( Population(P).objs - repmat(Z,T,1) )./W(P,:),[],2 );
                g_new = max( abs( repmat(Archive_und(i).obj,T,1) - repmat(Z,T,1) )./W(P,:),[],2 );
                Population(P(g_old > g_new)) = Archive_und(i);
            end
        end
    end

    %% Delet the poorly performed weights until the size of W is reduced to N
    % find out the solution that is shared by the most weights in the population
    while length(Population) > Global.N
        [~,ai,bi] = unique(Population.objs,'rows');
        if length(ai) == length(bi)   % If every solution in the population corresponds to only one weight
            % Normalisation
            fmax  = max(Population.objs,[],1);
            fmin  = min(Population.objs,[],1);
            PCObj = (Population.objs-repmat(fmin,length(Population),1))./repmat(fmax-fmin,length(Population),1);
            % Determine the radius of the niche
            d  = pdist2(PCObj,PCObj);
            d(logical(eye(length(d)))) = inf;
            sd = sort(d,2);
            num_obj = size(Population.objs,2);
            r  = median(sd(:,min(num_obj,size(sd,2))));
            R  = min(d./r,1);
            % Delete solution one by one
            while length(Population) > Global.N
                [~,worst]  = max(1-prod(R,2));
                Population(worst)  = [];
                R(worst,:) = [];
                R(:,worst) = [];
                W(worst,:) = [];
            end
        else
            Index = find(bi==mode(bi));
            Value_Tche2 = max(abs(Population(Index).objs-repmat(Z,size(Index,1),1))./W(Index,:),[],2);
            Index_max= find(Value_Tche2 == max(Value_Tche2));
            Population(Index(Index_max(1)))=[];
            W(Index(Index_max(1)),:)=[];
        end
    end
    % Update the neighbours of each weight
    B = pdist2(W,W);
    [~,B] = sort(B,2);
    B = B(:,1:T);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function truncated_normal_b_cdf_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_B_CDF_TEST tests TRUNCATED_NORMAL_AB_CDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_B_CDF_TEST\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_B_CDF evaluates the upper Truncated Normal CDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean = mu\n' );
  fprintf ( 1, '    standard deviation = sigma\n' );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval (-oo,b]\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                                                  Stored         Computed\n' );
  fprintf ( 1, '       X        Mu         S         B             CDF             CDF\n' );
  fprintf ( 1, '\n');

  n_data = 0;

  while ( 1 )

    [ n_data, mu, sigma, b, x, cdf1 ] ...
      = truncated_normal_b_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    cdf2 = truncated_normal_b_cdf ( x, mu, sigma, b );

    fprintf( 1, '  %8.1f  %8.1f  %8.1f  %8.1f  %14g  %14g\n', x, mu, sigma, b, cdf1, cdf2 );

  end

  return
end


text:
 function out = subsref(ob, arg)
%function out = subsref(ob, arg)
%	evaluate call of the form ob.ref or ob(ref,:)

%
%	ob.?
%
if arg.type == '.'
	out = getfield(struct(ob), arg.subs);

%
%	ob(?)
%
elseif arg.type == '()'

	if ~isempty(ob.index1) | ~isempty(ob.index2)
		error 'subscripts after subscripts perhaps not done'
	end

	sub = arg.subs;

	%
	%	G(:,:)
	%
	if length(sub) == 2 & sub{1} == ':' & sub{2} == ':'
		out = ob.G;

	%
	%	G(:,j)
	%
	elseif length(sub) == 2 & sub{1} == ':'
		jj = col(sub{2});
		if islogical(jj)
			if length(jj) ~= ob.dims(2)
				error 'bad column logical length'
			end
		else
			bad = jj < 1 | jj > ob.dims(2);
			if any(bad)
				printf('bad column indeces:')
				jj(bad)'
				error 'subsref problem'
			end
		end

		%
		%	handle G(:,mask(:)) as a special case
		%
		if length(jj) == ob.nx*ob.ny & all(jj == ob.mask(:)) ...
				& ~ob.is_masked & ~ob.is_transpose
			ob.is_masked = true;
			ob.dims(2) = sum(ob.mask(:));
			ob.G = ob.G(:,ob.mask(:));
			out = ob;

		else
%			G = ob.G;
			out = subsref(ob.G, arg);
		end

	%
	%	G(i,:)
	%
	elseif length(sub) == 2 & sub{2} == ':'
		out = subsref(ob.G, arg);

	else
		arg.subs
		error todo
	end

else
	error(sprintf('type %s notdone', arg.type))
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function cv = combine_visitors(varargin)
% COMBINE_VISITORS Generate a new visitor by combining existing visitors
%
% cv = combine_visitors(v1, v2, ...) generates a new algorithm visitor that
% works by calling the functions in v1 followed by the functions in v2, and
% so on.
%
% The value returned by the combined visitor function is the bitwise & of
% all the individual return values.  So, if any visitor requests the
% algorithm to halt, then the algorithm will halt.
%
% Note: using a combined visitor is somewhat slower than writing a custom
% combined visitor yourself.
%
% Example:
%    vis1 = struct();
%    vis1.examine_vertex = @(u) fprintf('vis1: examine_vertex(%i)\n', u);
%    vis2 = struct();
%    vis2.examine_vertex = @(u) fprintf('vis2: examine_vertex(%i)\n', u);
%    combined_vis = combine_visitors(vis1, vis2);
%    load graphs/bfs_example.mat
%    breadth_first_search(A,1,combined_vis);

% David Gleich
% Copyright, Stanford University, 2006-2008

% trivial output
if isempty(varargin)
    error('matlab_bgl:invalidParameter', 'combine_visitors requires at least one argument');
end

if length(varargin) == 1, cv = varargin{1}; return; end

cv_fn = struct();

for ii = 1:length(varargin)
    fn = fieldnames(varargin{ii});

    for jj = 1:length(fn)
        if ~isfield(cv_fn,fn{jj})
            cv_fn.(fn{jj}) = 1;
        else
            cv_fn.(fn{jj}) = cv_fn.(fn{jj}) + 1;
        end
    end
end

    function cfunc=combine_function(name,count)

        used_visitors = cell(count,1);
        cur_used = 1;

        for kk = 1:length(varargin)
            if isfield(varargin{kk}, name)
                used_visitors{cur_used} = varargin{kk}.(name);
                cur_used = cur_used + 1;
            end
        end

        function rval = combine_function_impl(varargin)
            rval = 1;
            for ll=1:length(used_visitors)
                try
                    stop = feval(used_visitors{ll}, varargin{:});
                    rval = rval & double(stop);
                catch
                    [lastmsg,lastid] = lasterr;
                    if ~strcmp(lastid, 'MATLAB:TooManyOutputs')
                        rethrow(lasterr);
                    else
                        feval(used_visitors{ll}, varargin{:});
                    end
                end
            end
        end

        cfunc = @combine_function_impl;

    end

cv = struct();

% specify only the visitors functions that only occured once
for ii = 1:length(varargin)
    fn = fieldnames(varargin{ii});

    for jj = 1:length(fn)
        if (cv_fn.(fn{jj}) == 1)
            cv.(fn{jj}) = varargin{ii}.(fn{jj});
        end
    end
end

% specify all visitors that occured more than once
fn = fieldnames(cv_fn);
for jj=1:length(fn)
    if (cv_fn.(fn{jj}) > 1)
        cv.(fn{jj}) = combine_function(fn{jj}, cv_fn.(fn{jj}));
    end
end

% overall function end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% print_n_updates - example hook function for each update
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%

function err = print_n_updates (R, args)

if R.verbose == 1
    fprintf(2, '# updates: %d\n', R.iteration.n_updates);
end

err = 0;

author: kyunghyuncho

repo: deepmat


text:
function [ ret ] = setPin2On( t_Kuka )
%% About:
% This function is used to set on the Pin 2 of the KUKA iiwa 7 R 800.

%% Syntax:
% [ ret ] = setPin2On( t_Kuka )

%% Arreguments:
% t_Kuka: is the TCP/IP connection

%% Return value:
% ret: a boolean variable
% true: if the message is treated successfully
% false: if an error has happened

% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='pin2on';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
ret=checkAcknowledgment(message);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function npts = plotObj(prob,xb,data)
%PLOTOBJ Plot the objective function contour
%   plotObj(prob,xb,data)

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Contour Colour
dkg = [0.4 0.4 0.4];
%Gather Commonly Used Inputs
npts = data.npts;
scale = data.scale;
idx = data.idx;
fval = [];

%Detail based on problem type, only if not specified by the user
if(isempty(npts))
    switch(lower(prob.type))
        case {'lp','milp','bilp'}
            npts = 5;
        case {'qp','qcqp','miqp','miqcqp','sdp','misdp'}
            npts = 30;
        case {'nls','uno','nlp','minlp'}
            npts = 50;
        otherwise
            npts = 50;
    end
end

%Determine if we have a 1D or 2D problem
if((~isempty(prob.lb) && length(prob.lb) == 1) || (~isempty(xb) && length(xb) == 1))
    %Generate vector Based On Mode
    switch(data.mode)
        case {'normal','usex0','bounded_scale'}
            if(length(scale) == 2)
                x = linspace(scale(1),scale(2),npts);
            elseif(length(scale) == 1)
                x = linspace(xb(1)-scale,xb(1)+scale,npts);
            else
                error('Input ''scale'' should be a vector with two inputs [xmin,xmax] for 1D problems');
            end
        case {'multi','bounded'}
            x = linspace(prob.lb(1),prob.ub(1),npts);
    end
    obj = zeros(npts,1);
    %Get Objective (created in buildConfig)
    fun = prob.objective;
    %Create Objective Line
    for i = 1:npts
    	obj(i) = fun(x(i));
    end
    %Do Log Plot if asked
    if(data.dolog)
        obj = log(obj);
    end
    %Draw Plot
    plot(x,obj.*prob.sense,'-','color',dkg);
    yl = ylim; axis([x(1) x(end) yl]);
    %Get Min
    if(~isempty(xb))
        fval = fun(xb).*prob.sense;
    end

else %2D or higher Problem
    %Generate Grid Based On Mode
    switch(data.mode)
        case {'normal','usex0','bounded_scale'}
            if(length(scale) == 4)
                [x1,x2] = meshgrid(linspace(scale(1),scale(2),npts),linspace(scale(3),scale(4),npts));
            elseif(length(scale) == 1)
                [x1,x2] = meshgrid(linspace(xb(idx(1))-scale,xb(idx(1))+scale,npts),linspace(xb(idx(2))-scale,xb(idx(2))+scale,npts));
            else
                error('Input ''scale'' should be a vector with four inputs [x1min, x1max, x2min, x2max] for 2D problems');
            end
        case {'multi','bounded'}
            [x1,x2] = meshgrid(linspace(prob.lb(idx(1)),prob.ub(idx(1)),npts),linspace(prob.lb(idx(2)),prob.ub(idx(2)),npts));
    end
    nox = size(x1);
    noy = size(x2);
    obj = zeros(nox(1),noy(2));
    %Get Objective (created in buildConfig)
    fun = prob.objective;
    x = data.fixval;
    %Create Objective Surface
    for npts = 1:nox(1)
        for m = 1:noy(2)
            x(idx) = [x1(npts,m) x2(npts,m)]';
            obj(npts,m) = fun(x);
        end
    end
    %Do Log Plot if asked
    if(data.dolog)
        obj = log(obj);
    end
    %Draw Contour Plot
    [hc,hl] = contour(x1,x2,obj.*prob.sense,':','color',dkg);
    if(~isempty(hc))
        clabel(hc,hl);
        hold on;
        %Plot Minimum Contour
        if(~isempty(xb))
            fval = fun(xb).*prob.sense;
            contour(x1,x2,obj.*prob.sense,':','color',dkg,'levellist',fval);
        end
    else
        optiwarn('OPTI:EmptyContour','Objective Contour Data is Empty - Cannot Plot Objective!');
        %Manually set axis of interest
        try
            axis([min(min(x1)) max(max(x1)) min(min(x2)) max(max(x2))]);
        catch
        end
    end

    xlabel(sprintf('x_%d',idx(1))); ylabel(sprintf('x_%d',idx(2)));
end

%Plot Optimum + Title
data.fval = fval;
title(plotTitle(prob,xb,data));
hold off;
end


text:
function z=vecr(fit)

z = fit.fit_points.kappa;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [xnew, Vnew, loglik, VVnew] = kalman_update(A, C, Q, R, y, x, V, varargin)
% KALMAN_UPDATE Do a one step update of the Kalman filter
% [xnew, Vnew, loglik] = kalman_update(A, C, Q, R, y, x, V, ...)
%
% INPUTS:
% A - the system matrix
% C - the observation matrix
% Q - the system covariance
% R - the observation covariance
% y(:)   - the observation at time t
% x(:) - E[X | y(:, 1:t-1)] prior mean
% V(:,:) - Cov[X | y(:, 1:t-1)] prior covariance
%
% OPTIONAL INPUTS (string/value pairs [default in brackets])
% 'initial' - 1 means x and V are taken as initial conditions (so A and Q are ignored) [0]
% 'u'     - u(:) the control signal at time t [ [] ]
% 'B'     - the input regression matrix
%
% OUTPUTS (where X is the hidden state being estimated)
%  xnew(:) =   E[ X | y(:, 1:t) ]
%  Vnew(:,:) = Var[ X(t) | y(:, 1:t) ]
%  VVnew(:,:) = Cov[ X(t), X(t-1) | y(:, 1:t) ]
%  loglik = log P(y(:,t) | y(:,1:t-1)) log-likelihood of innovatio

% set default params
u = [];
B = [];
initial = 0;

args = varargin;
for i=1:2:length(args)
  switch args{i}
   case 'u', u = args{i+1};
   case 'B', B = args{i+1};
   case 'initial', initial = args{i+1};
   otherwise, error(['unrecognized argument ' args{i}])
  end
end

%  xpred(:) = E[X_t+1 | y(:, 1:t)]
%  Vpred(:,:) = Cov[X_t+1 | y(:, 1:t)]

if initial
  if isempty(u)
    xpred = x;
  else
    xpred = x + B*u;
  end
  Vpred = V;
else
  if isempty(u)
    xpred = A*x;
  else
    xpred = A*x + B*u;
  end
  Vpred = A*V*A' + Q;
end

e = y - C*xpred; % error (innovation)
n = length(e);
ss = length(A);
S = C*Vpred*C' + R;
Sinv = inv(S);
ss = length(V);
loglik = gaussian_prob(e, zeros(1,length(e)), S, 1);
K = Vpred*C'*Sinv; % Kalman gain matrix
% If there is no observation vector, set K = zeros(ss).
xnew = xpred + K*e;
Vnew = (eye(ss) - K*C)*Vpred;
VVnew = (eye(ss) - K*C)*A*V;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function boundary = p06_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P06_BOUNDARY_NEAREST returns a nearest boundary point in problem 06.
%
%  Discussion:
%
%    The given input point need not be inside the region.
%
%    In some cases, more than one boundary point may be "nearest",
%    but only one will be returned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  r1 = 1.0;
  r2 = 0.5;

  for j = 1 : n

    x = point(1,j);
    y = point(2,j);

    if ( x == 0.0 & y == 0.0 )

      boundary(1:2,j) = [ r2, 0.0 ]';
%
%  Determine the angle formed by (0,0) and the point.
%
    else

      t = atan4 ( y, x );
%
%  Find the nearest point on the superellipse x^4 + y^4 = 1^4.
%
      t1 = t - pi / 4.0;
      t2 = t + pi / 4.0;

      status = 0;
      dstar1 = 0.0;

      while ( 1 )

        [ t1, t2, tstar1, status ] = fmin_rc ( t1, t2, status, dstar1 );

        if ( status == 0 )
          break
        end

        cm = abs ( cos ( tstar1 ) );
        cs = r8_sign ( cos ( tstar1 ) );
        sm = abs ( sin ( tstar1 ) );
        ss = r8_sign ( sin ( tstar1 ) );

        dstar1 = ( x - r1 * cs * sqrt ( cm ) ).^2 ...
               + ( y - r1 * ss * sqrt ( sm ) ).^2;

      end

      boundary(1,j) = r1 * cs * sqrt ( cm );
      boundary(2,j) = r1 * ss * sqrt ( sm );
%
%  Find the nearest point on the superellipse x^4 + y^4 = 1/2^4.
%
      t1 = t - pi / 4.0;
      t2 = t + pi / 4.0;

      status = 0;
      dstar2 = 0.0;

      while ( 1 )

        [ t1, t2, tstar2, status ] = fmin_rc ( t1, t2, status, dstar2 );

        if ( status == 0 )
          break
        end

        cm = abs ( cos ( tstar2 ) );
        cs = r8_sign ( cos ( tstar2 ) );
        sm = abs ( sin ( tstar2 ) );
        ss = r8_sign ( sin ( tstar2 ) );

        dstar2 = ( x - r2 * cs * sqrt ( cm ) ).^2 ...
               + ( y - r2 * ss * sqrt ( sm ) ).^2;

      end
%
%  Because of some MATLAB idiocy I don't have the patience to figure out,
%  involving SQRT having the wrong number of arguments(!),
%  I must replace the line
%    boundary(1:2,j) = [ r2 * cs * sqrt ( cm ), r2 * ss * sqrt ( sm ) ]';
%
      if ( dstar2 < dstar1 )
        boundary(1,j) = r2 * cs * sqrt ( cm );
        boundary(2,j) = r2 * ss * sqrt ( sm );
      end

    end

  end

  return
end


text:
function onoff = emgonoff(rawemg, fs, ws, sd)
% EMGONOFF - Find on/off times and indicies of raw EMG data.
% Calculate average(mean) value of resting EMG
% Define "on" EMG as the sample where average value of EMG in a given window
% range around the sample is a given # of std. dev. above avg. resting EMG.
%
% onoff = emgonoff(rawemg, fs, ws, sd)
%
% Use the mouse to select two ranges of "resting" EMG from a graph of the
% full-wave rectified EMG data.  Click four times: start and end of 1st
% resting range, and start and end of 2nd resting range.  Mouse clicks need
% to be consecutive and in order of increasing time (i.e. left-to-right on
% the graph).  The first range should precede the EMG burst associated with
% the muscle contraction under consideration.  The sedond range should be
% the resting EMG data immediately following the EMG burst.
%
% rawemg = input file raw emg data (1-column vector)
% fs = sampling rate of raw EMG data in Hz
% ws = window size in milliseconds (50ms @ 2400Hz = 120 samples)
% sd = number of std. deviations above resting rms emg to trigger an "ON"
% Default values:
%   ws = 50ms
%   sd = 1

% Algorythm for EMG onset & offset taken from
% Hodges, P.W. and B.H. Bui, _A comparison of computer-based methods for
% the determination of onset of muscle contraction using electromyography._
% Electroencephalography & Clinical Neurophysiology, 1996. 101(6): p. 511-9
%
% Created by: Kieran A. Coghlan, BSME, MSES
% SUNY at Buffalo, New York
% <kc_news@sonic.net>
% Last modified: 1 May, 2006

%% Check inputs for defaults
if nargin < 2, error('Not enough inputs. Type "help emgonoff" for help.'); end
if nargin < 3, ws = 50; end
if nargin < 4, sd = 1; end;
%% Full-Wave-Rectify the raw data
fwlo = abs(rawemg(:,1));
%% prepare for loop
% Get two ranges for resting emg (before & after burst) using ginput
R = input('\nUse the mouse to select FOUR points to define the begining and \nend of two data ranges that will be used to calculate average\nresting EMG values before and after the EMG burst (muscle contraction)\nPress [RETURN] to begin: ');
clear R;
f1 = figure;
plot(fwlo);
[x y] = ginput(4); %click four times: two for start/end of resting emg before burst two for resting emg after burst
x = round(x);
clear y;
%close(f1);clear f1; % Leave commented if you want to keep EMG graph up to
%                    % do a visual QA of on/off results.
%% preallocate arrays
mvgav = zeros(x(4)-x(1),1);
onoff(1,1) = 0;
i=0;
restav = mean(fwlo(x(1):x(2))); %average value of rest EMG before ON
reststd = std(fwlo(x(1):x(2))); %std. dev. of rest EMG before ON
restav2 = mean(fwlo(x(3):x(4))); %average value of rest EMG after OFF
reststd2 = std(fwlo(x(3):x(4))); %std. dev. of rest EMG after OFF
%% window size (in samples) = ws*fs e.g. 50ms*2400Hz = 120 samples
sws2 = fs*(0.001*ws);
sws = 0.5*(sws2);
sws = round(sws);
%% find "ON" index:
% for xi, change from x(1) to x(2) if you want to ignore any "blips"
% within the resting range.
%xi = x(1);
xi = x(2);
xi = round(xi);
for n = 2:length(mvgav);
    mvgav(n,1) = mean(fwlo((xi-sws):(xi+sws)));
    if mvgav(n) > restav+sd*reststd;
        i = i+1;
        onoff(i,1) = xi;
        break
    end
    xi = xi+1;
end
%% find "OFF" index:
clear n xi i
mvgav2=zeros(x(4)-x(1),1);
i=0;
xi=onoff(1,1)+(1/2)*(x(3)-onoff(1,1)); %start OFF search approx. 1/2 way through ON burst.
%% OFF loop:
xi=round(xi);
for n=2:length(mvgav2);
    mvgav2(n,1)=mean(fwlo((xi-sws):(xi+sws)));
    if mvgav2(n)<restav2+sd*reststd2;
        i=i+1;
        onoff(i,2)=xi;
        break
    end
    xi=xi+1;
end


text:
% CVX: Internal functions and scripts.
%    This directory contains code that is meant for internal use by
%    the CVX system itself. Documentation of the functions in this
%    directory is more sparse and not intended for end users.

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function out = poly(f)
%POLY   Polynomial coefficients of a CHEBTECH.
%   C = POLY(F) returns the polynomial coefficients of F so that
%           F(x) = C(1)*x^N + C(2)*x^(N-1) + ... + C(N)*x + C(N+1)
%
%   Note that unlike the MATLAB POLY command, CHEBTECH/POLY can operate on
%   array-valued CHEBTECH objects, and hence produce a matrix output. In such
%   instances, the rows of C correspond to the columns of F = [F1, F2, ...].
%   That is,
%           F1(x) = C(1,1)*x^N + C(1,2)*x^(N-1) + ... + C(1,N)*x + C(1,N+1)
%           F2(x) = C(2,1)*x^N + C(2,2)*x^(N-1) + ... + C(2,N)*x + C(2,N+1).
%   This strange behaviour is a result of MATLAB's decision to return a row
%   vector from the POLY command, even for column vector input.
%
% See also CHEBCOEFFS, TRIGCOEFFS, LEGCOEFFS.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Mathematical reference]: Section 3.3 Mason & Handscomb, "Chebyshev
% Polynomials". Chapman & Hall/CRC (2003).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Deal with empty case:
if ( isempty(f) )
    out = [];
    return
end

% Flip the Chebyshev coefficients to match Matlab ordering:
coeffs = flipud(f.coeffs);
[n, m] = size(coeffs);

% Coefficients on the unit interval:
if ( n == 1 )
    % Constant case:
    out = coeffs;

elseif ( n == 2 )
    % Linear case:
    out = coeffs.';

else
    % General case:

    % Initialise working vectors:
    tn = zeros(m, n);
    tnold1 = tn;
    tnold1(:,2) = 1;
    tnold2 = tn;
    tnold2(:,1) = 1;
    zer = zeros(m, 1);

    % Transpose coeffs (as output is in this form):
    coeffs = coeffs.';

    % Initialise output:
    out = zeros(m, n);

    % Initial step:
    out(:,[1,2]) = [ zer coeffs(:,end).*tnold2(:,1) ] + ...
        bsxfun(@times, coeffs(:,end-1), tnold1(:,[2,1]));

    % Recurrence:
    for k = 3:n
        tn(:,1:k)  = [ zer 2*tnold1(:,1:k-1) ] - [ tnold2(:,1:k-2), zer, zer ];
        out(:,1:k) = bsxfun(@times, coeffs(:,end-k+1), tn(:,k:-1:1)) + ...
                        [ zer, out(:,1:k-1) ];
        tnold2 = tnold1;
        tnold1 = tn;
    end

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%
% Test for c-transform and Laguerre cells in 1D.

addpath('../toolbox/');

rep = 'results/c-transform-1d/';
[~,~] = mkdir(rep);

SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1]);

% # input points
N = 7;
% Grid size
q = 512;
rng('default'); rng(646);

M = rand(1,N);

% grid
t = linspace(0,1,q)';

% potentials
u = 0.1*randn(N,1);

% exponent/scale for the transform
if not(exist('p'))
    p = 2;
end
rho = 2;
if not(exist('epsilon'))
    epsilon = .01;
end

%%
% study influence of p

epsilon = 0;
p_list = linspace(.5,2,13);
lw = 2;
clf; hold on;
Col = [];
for ip=1:length(p_list)
    p = p_list(ip);
    s = (ip-1)/(length(p_list)-1);
    col = [0,s,1-s];  Col(end+1,:) = col;
    % v(x,y) = min_i c((x,y),m_i) - u_i
    D = abs( repmat(M(:), [1 q]) - repmat(t(:)', [N 1])  );
    S = rho * D.^p - repmat(u,[1 q]);
    %
    [v,I] = min(S);
    if epsilon>0
        S = S - repmat(v, [size(S,1) 1]);
        v = - epsilon * log( sum(  exp( -S/epsilon ) ) ) + v;
    end
    %
    plot(t, v, 'LineWidth', lw, 'Color', col);
end
%
plot(M, -u, 'r.', 'MarkerSize', 25);
axis tight; box on;
colormap(Col); caxis([p_list(1),p_list(end)]);
set(gca, 'FontSize', 20);
colorbar;
SetAR(1/2);
saveas(gcf, [rep 'c-transform-p.png'], 'png');
saveas(gcf, [rep 'c-transform-p.eps'], 'epsc');

%%
% study influence of epsilon

p = 1;
eps_list = linspace(0,.3,13);
lw = 2;
clf; hold on;
Col = [];
for ieps=1:length(eps_list)
    epsilon = eps_list(ieps);
    s = (ieps-1)/(length(eps_list)-1);
    col = [s,s,1-s];  Col(end+1,:) = col;
    % v(x,y) = min_i c((x,y),m_i) - u_i
    D = abs( repmat(M(:), [1 q]) - repmat(t(:)', [N 1])  );
    S = rho * D.^p - repmat(u,[1 q]);
    %
    [v,I] = min(S);
    if epsilon>0
        S = S - repmat(v, [size(S,1) 1]);
        v = - epsilon * log( sum(  exp( -S/epsilon ) ) ) + v;
    end
    %
    plot(t, v, 'LineWidth', lw, 'Color', col);
end
%
plot(M, -u, 'r.', 'MarkerSize', 25);
axis tight; box on;
colormap(Col); caxis([eps_list(1),eps_list(end)]);
set(gca, 'FontSize', 20);
colorbar;
SetAR(1/2);
saveas(gcf, [rep 'c-transform-eps.png'], 'png');
saveas(gcf, [rep 'c-transform-eps.eps'], 'epsc');

author: optimaltransport


text:
function divdif_test21 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST21tests DIF_BASIS_DERIV;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 June 2013
%
%  Author:
%
%    John Burkardt
%
  nd = 3;
  xd = [ -2.0, 1.0, 5.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST21\n' );
  fprintf ( 1, '  DIF_BASIS_DERIV computes difference tables for\n' );
  fprintf ( 1, '  the first derivative of each Lagrange basis.\n' );

  [ xdp, ddp ] = dif_basis_deriv ( nd, xd );
%
%  Because the difference tables were shifted to all 0 abscissas,
%  they contain the polynomial coefficients.
%
  r8mat_transpose_print ( nd - 1, nd, ddp, ...
    '  Lagrange basis derivative polynomial coefficients:' );

  c = dif_to_r8poly ( nd - 1, xdp, ddp(1:nd-1,1) );
  r8poly_print ( nd - 1, c, '  P1''=-(2x-6)/21' );

  c = dif_to_r8poly ( nd - 1, xdp, ddp(1:nd-1,2) );
  r8poly_print ( nd - 1, c, '  P2''=-(2x-3)/12' );

  c = dif_to_r8poly ( nd - 1, xdp, ddp(1:nd-1,3) );
  r8poly_print ( nd - 1, c, '  P3''=(2x+1)/28' );

  return
end


text:
function [volCc,dims] = mrAnatGetBrainVolume(brainMask)
%
% [volCc,dims] = mrAnatGetbrainVolume([brainMask=uigetfile])
%
% Returns the brain volume (in cubic centimeters) for the specified brain
% mask. Optionally also returns the left-right, anterior-posterior and
% superior-inferior dimensions (in centimeters) of the brain.
%
% If you don't have a brain mask, use mrAnatExtractBrain.
%
% 2008.08.14 RFD wrote it.
% 2008.09.09 RFD fixed x,y bug in dims- the first two dims were flipped.

if(~exist('brainMask','var')||isempty(brainMask))
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a brain mask NIFTI file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    brainMask = fullfile(p,f);
end

if(ischar(brainMask))
    % It's a nifti file
    brainMask = niftiRead(brainMask);
end

mm = brainMask.pixdim(1:3);
ccPerPixel = prod(mm)/1000;
volCc = numel(find(brainMask.data)) * ccPerPixel;

if(nargout>1)
    tmp = sum(brainMask.data,3);
    y = find(sum(tmp,1));
    y = [y(1) y(end)];
    x = find(sum(tmp,2));
    x = [x(1) x(end)];
    tmp = squeeze(sum(brainMask.data,1));
    z = find(sum(tmp,1));
    z = [z(1) z(end)];
    dims = [diff(x)*mm(1) diff(y)*mm(2) diff(z)*mm(3)];
    dims = dims./10;
end

return;

bd = '/biac3/wandell4/data/reading_longitude/dti_y1';
d = dir(fullfile(bd,'*04*'));
n = 0;
for(ii=1:numel(d))
    f = fullfile(bd,d(ii).name,'t1','t1_mask.nii.gz');
    if(exist(f,'file'))
        n = n+1;
        [volCc(n),dims(n,:)] = mrAnatGetBrainVolume(f);
        sc{n} = d(ii).name;
    end
end

[behData,colNames] = dtiGetBehavioralData(sc);

sexInd = strmatch('Sex (1=male)',colNames);
boys = behData(:,sexInd)==1;
fprintf('boys = %0.1fcc (%0.2f), girls = %0.1fcc (%0.2f)\n', mean(volCc(boys)), std(volCc(boys)), mean(volCc(~boys)), std(volCc(~boys)));

brInd = strmatch('Basic reading (W-J)',colNames);
necrosisControls = ~boys & behData(:,brInd)>=90;
mn = mean(volCc(necrosisControls));
sd = std(volCc(necrosisControls));

[sVol, sDims] = mrAnatGetBrainVolume('/biac1/wandell/data/radiationNecrosis/dti/al060406/t1/t1_mask.nii.gz');
% Compute the percentile:
[hcnt,hcent] = hist(volCc(necrosisControls),50);
pctile = sum(hcnt(find(sVol>hcent)))/sum(hcnt)*100;
fprintf('Controls: mean = %0.0fcc (%0.1f); S = %0.0fcc (%0.0f %%tile)\n', mn, sd, sVol, pctile);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef PTKLungSurface < PTKPlugin
    % PTKLungSurface. Plugin for finding points around the surface of the lungs
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKLungSurface runs the library function PTKGetSurfaceFromSegmentation
    %     to find the lung surface and returns this as an image.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lung Surface'
        ToolTip = 'Segment the lung surface'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, ~)
            lung_mask = dataset.GetResult('PTKLeftAndRightLungs');
            results = lung_mask.BlankCopy;
            lungs = PTKGetSurfaceFromSegmentation(uint8(lung_mask.RawImage > 0));
            results.ChangeRawImage(lungs);
            results.ImageType = PTKImageType.Colormap;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit