text:
function val=logOnePlusXMinusX(x)
%%LOGONEPLUSXMINUSX Evaluate log(1+x)-x reducing the loss of precision that
%            occurs for values of x that are very close to 0 for real
%            values of x.
%
%INPUTS: x A matrix of real values>=-1.
%
%OUTPUTS: val The values log(1+x)-x.
%
%The function first computes z=log(1+x) using log1p(x) to reduce finite
%precision errors. This value is put into the expSlope2 function, which
%computes 2*(exp(x)-1-x)/x^2 without a low loss of finite precision.
%Multiplying the result by (-1/2)*z^2, one gets the expression log(1+x)-x.
%This approach ovoids a direct subtraction of x, which is where the finite
%precision problems occur.
%
%EXAMPLE:
%Here, we show the improvement that this algorithm offers to a direct
%evaluation.
% x=1e-150;
% valAccurate=logOnePlusXMinusX(x)
% valInaccurate=log1p(x)-x
%One will see that valAccurate is about -5e-301, whereas valInaccurate is
%just zero. The value -5e-301 can be verified using greatly extended
%precision arithmetic.
%
%October 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(any(x(:)<-1)||~isreal(x))
       error('This function only supports real values of x >=-1.');
    end

    val=zeros(size(x));

    selNeg=x<-0.99;
    xNeg=x(selNeg);
    xPos=x(~selNeg);
    if(~isempty(xNeg))
        valNeg=log1p(x)-x;
    else
        valNeg=[];
    end

    if(~isempty(xPos))
        z=log1p(x);
        e2Val=expSlope2(z);
        valPos=-(1/2)*e2Val.*z.^2;
    else
        valPos=[];
    end

    val(selNeg)=valNeg;
    val(~selNeg)=valPos;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [f,g,h] = fminunc_wrapper(x,F,G,H)
% [f,g,h] = fminunc_wrapper( x, F, G, H )
% for use with Matlab's "fminunc"
f = F(x);
if nargin > 2 && nargout > 1
    g = G(x);
end
if nargin > 3 && nargout > 2
    h = H(x);
end

repo: Deep-ADMM-Net


text:
function openInplaneWindowforAlign3
%
% Opens a new inplane window and initializes the corresponding data structure.
%
% djh 8/29/01, modified from openInplaneWindow

% Make sure the global variables exist
global mrSESSION INPLANE

disp('Initializing Inplane view')

% Set name, viewType, & subdir
INPLANE.name='INPLANE';
INPLANE.viewType='Inplane';
INPLANE.subdir='Inplane';

% Refresh function, gets called by refreshScreen
INPLANE.refreshFn = 'refreshView';

% Initialize slot for anat
INPLANE.anat = [];

% Initialize slots for co, amp, and ph
INPLANE.co = [];
INPLANE.amp = [];
INPLANE.ph = [];
INPLANE.map = [];
INPLANE.mapName = '';
INPLANE.spatialGrad = [];

% Initialize slots for tSeries
INPLANE.tSeries = [];
INPLANE.tSeriesScan = NaN;
INPLANE.tSeriesSlice = NaN;

% Initialize ROIs
INPLANE.ROIs = [];
INPLANE.selectedROI = 0;

% Initialize curDataType
INPLANE.curDataType = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make displayModes and color maps %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize displayModes
INPLANE=resetDisplayModes(INPLANE);
INPLANE.ui.displayMode='anat';

%%%%%%%%%%%%%%%%%%%
% Open the window %
%%%%%%%%%%%%%%%%%%%

% Figure number for inplane view window
INPLANE.ui.figNum=figure('MenuBar','none');

% Handle for inplane view window
INPLANE.ui.windowHandle = gcf;

% Handle for main axis of inplane view
INPLANE.ui.mainAxisHandle = gca;

% Adjust position and size of main axes so there's room for
% buttons, colorbar, and sliders
set(INPLANE.ui.mainAxisHandle,'position',[0.1 0.1 0.725 0.7]);

% Set minColormap property so there's potentially room for 128
% colors
set(INPLANE.ui.windowHandle,'minColormap',128)

% Sharing of colors seems like it might be OK, but I'm turning it
% off just to be sure (djh, 1/26/98).
set(gcf,'sharecolors','off');

%%%%%%%%%%%%%%%%%%%%%%%%%
% Add Annotation String %
%%%%%%%%%%%%%%%%%%%%%%%%%

INPLANE = makeAnnotationString(INPLANE);

%%%%%%%%%%%%%%%%%
% Add Color Bar %
%%%%%%%%%%%%%%%%%

% Make color bar and initialize it to 'off'
INPLANE.ui.colorbarHandle=makeColorBar(INPLANE);
setColorBar(INPLANE,'off');
INPLANE.ui.cbarRange = [];

%%%%%%%%%%%%%%%%%%%
% Add popup menus %
%%%%%%%%%%%%%%%%%%%

disp('Attaching popup menus')

INPLANE = makeROIPopup(INPLANE);
INPLANE = makeDataTypePopup(INPLANE);

%%%%%%%%%%%%%%%
% Add sliders %
%%%%%%%%%%%%%%%

disp('Attaching sliders')

% scan slider
w = 0.12; h = 0.04; l = 0; b = 0.95;
INPLANE = makeSlider(INPLANE,'scan',[],[l b w h]);
INPLANE = initScanSlider(INPLANE,1);
INPLANE = selectDataType(INPLANE,INPLANE.curDataType);

% slice slider
w = 0.12; h = 0.04; l = 0; b = 0.85;
INPLANE = makeSlider(INPLANE,'slice',[],[l b w h]);
INPLANE = initSliceSlider(INPLANE);

% correlation threshold:
INPLANE = makeSlider(INPLANE,'cothresh',[0,1],[.85,.85,.15,.05]);
setCothresh(INPLANE,0);

% phase window:
INPLANE = makeSlider(INPLANE,'phWinMin',[0,2*pi],[.85,.75,.15,.05]);
INPLANE = makeSlider(INPLANE,'phWinMax',[0,2*pi],[.85,.65,.15,.05]);
setPhWindow(INPLANE,[0 2*pi]);

% parameter map window:
INPLANE = makeSlider(INPLANE,'mapWinMin',[0,1],[.85,.55,.15,.05]);
INPLANE = makeSlider(INPLANE,'mapWinMax',[0,1],[.85,.45,.15,.05]);
setMapWindow(INPLANE,[0 1]);

% anatClip: determines clipping of the anatomy base-image
%           values to fill the range of available grayscales.
INPLANE = makeSlider(INPLANE,'anatMin',[0,1],[.85,.2,.15,.05]);
INPLANE = makeSlider(INPLANE,'anatMax',[0,1],[.85,.1,.15,.05]);
setAnatClip(INPLANE,[0 .5]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize display parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize image field
INPLANE.ui.image = [];

% Show all ROIs
INPLANE.ui.showROIs = 2;

%%%%%%%%%%%%%%%%%%%%%%%%%
% Load user preferences %
%%%%%%%%%%%%%%%%%%%%%%%%%

INPLANE = loadPrefs(INPLANE);

%%%%%%%%%%%%%%%%%%
% Refresh screen %
%%%%%%%%%%%%%%%%%%

INPLANE=refreshScreen(INPLANE);

disp('Done initializing Inplane view')

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function bnet = bnet_from_engine(engine)
% BNET_FROM_ENGINE Return the bnet structure stored inside the engine (inf_engine)
% bnet = bnet_from_engine(engine)

bnet = engine.bnet;

% We cannot write 'engine.bnet' without writing a 'subsref' function,
% since engine is an object with private parts.
% The bnet field should be the only thing external users of the engine should need access to.
% We do not pass bnet as a separate argument, since it could get out of synch with the one
% encoded inside the engine.

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
fprintf('\nHere we train a two-layer Deep Belief Network\non the MNIST dataset.');

% LOAD THE MINIST DATASET
load mnistSmall.mat

trainLabels = single(trainLabels);
testLabels = single(testLabels);
[nObs,nIn]= size(trainData);
[~,nOut]=size(trainLabels);

nHid = 256;

% DEFINE MULTI-LAYERED ARCHITECTURE
arch.size = [nIn nHid nHid];
arch.lRate = [.1 .1];
arch.nEpoch = [200 200];

% GLOBAL OPTIONS
arch.opts = {'inputType', 'binary', ...
			 'classifier',true};

% INITIALIZE NETWORK
d = dbn(arch);

% TRIAN & TEST THE DBN
fprintf('\n\nTraining a two-layer Deep Belief Network...\n');
d = d.train(trainData,trainLabels);

[~,classErr1,misClass] = d.classify(testData,testLabels);
fprintf('\nClassification error: %1.2f%%',classErr1*100);

misClass = testData(misClass,:);
close all
figure;
visWeights(misClass',0,[0 1]); title(sprintf('Missclassified -- Error=%1.2f %%',classErr1*100));
drawnow

fprintf('\n\nNow we fine-tune the DBN using backprop\n')

% USE WEIGHTS TO INITIALIZE NEURAL NETWORK
ftArch.lRate = 0.01;
ftArch.opts = {'costFun','xent', 'nEpoch',30,'sparsity',0.01};

n = d.fineTune(trainData,trainLabels,ftArch);
[~,classErr2] = n.test(testData,testLabels,'classerr');
fprintf('\nFine-tuned classification error: %1.2f%%',classErr2*100);
fprintf('\n...an improvement of %1.2f%%',(classErr1-classErr2)*100);

figure;
nVis = 100;
cLims = [-1.5 1.5];
subplot(131);
visWeights(d.rbmLayers{1}.W(:,1:nVis),cLims);
title('Features Selected from 1st Layer RBM')
drawnow

subplot(132);
visWeights(n.layers{1}.W(1:nVis,:)',cLims);
title('Fine-Tuned First-layer Weights')

subplot(133);
visWeights(d.rbmLayers{1}.W(:,1:nVis)-n.layers{1}.W(1:nVis,:)',cLims);
title('Fine-Tuned First-layer Features')
title('Difference')

author: dustinstansbury

repo: medal


text:
classdef Hyperplane < handle
    properties
        p_;    % point on the plane
        n_;    % Normal of the plane, directional,
        local_bbox_;
    end
    methods
        function obj = Hyperplane(p, n)
            obj.p_ = p;
            obj.n_ = n;
        end
        % Calculate the signed distance from point
        function dis = signed_dist(obj, pt)
            dis = dot(obj.n_, pt - obj.p_);
        end
        % Calculate the distance from point
        function dis = dist(obj, pt)
            dis = abs(obj.signed_dist(pt));
        end
    end
end

author: LenaShengzhen

repo: AerialRobotics


text:
function [Cxy, f] = mycohere(varargin)
%COHERE Coherence function estimate.
%   Cxy = COHERE(X,Y,NFFT,Fs,WINDOW) estimates the coherence of X and Y
%   using Welch's averaged periodogram method.  Coherence is a function
%   of frequency with values between 0 and 1 that indicate how well the
%   input X corresponds to the output Y at each frequency.  X and Y are
%   divided into overlapping sections, each of which is detrended, then
%   windowed by the WINDOW parameter, then zero-padded to length NFFT.
%   The magnitude squared of the length NFFT DFTs of the sections of X and
%   the sections of Y are averaged to form Pxx and Pyy, the Power Spectral
%   Densities of X and Y respectively. The products of the length NFFT DFTs
%   of the sections of X and Y are averaged to form Pxy, the Cross Spectral
%   Density of X and Y. The coherence Cxy is given by
%       Cxy = (abs(Pxy).^2)./(Pxx.*Pyy)
%   Cxy has length NFFT/2+1 for NFFT even, (NFFT+1)/2 for NFFT odd, or NFFT
%   if X or Y is complex. If you specify a scalar for WINDOW, a Hanning
%   window of that length is used.  Fs is the sampling frequency which does
%   not effect the cross spectrum estimate but is used for scaling of plots.
%
%   [Cxy,F] = COHERE(X,Y,NFFT,Fs,WINDOW,NOVERLAP) returns a vector of freq-
%   uencies the same size as Cxy at which the coherence is computed, and
%   overlaps the sections of X and Y by NOVERLAP samples.
%
%   COHERE(X,Y,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none',
%   specifies a detrending mode for the prewindowed sections of X and Y.
%   DFLAG can take the place of any parameter in the parameter list
%   (besides X and Y) as long as it is last, e.g. COHERE(X,Y,'mean');
%
%   COHERE with no output arguments plots the coherence in the current
%   figure window.
%
%   The default values for the parameters are NFFT = 256 (or LENGTH(X),
%   whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2,
%   P = .95, and DFLAG = 'none'.  You can obtain a default parameter by
%   leaving it off or inserting an empty matrix [], e.g.
%   COHERE(X,Y,[],10000).
%
%   See also PSD, CSD, TFE.
%   ETFE, SPA, and ARX in the Identification Toolbox.

%   Author(s): T. Krauss, 3-31-93
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 1998/06/03 14:42:19 $

narginchk(2,7);
x = varargin{1};
y = varargin{2};
%[msg,nfft,Fs,window,noverlap,p,dflag]=psdchk(varargin(3:end),x,y);
%error(msg)
nfft = 256;
Fs =8000;
window=hanning(nfft);
noverlap = 0.75*nfft;
p = .95;
dflag = 'none';

% compute PSD and CSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x , y if length is less than the window length
    x(nwind)=0;
    y(nwind)=0;
    n=nwind;
end
x = x(:);		% Make sure x is a column vector
y = y(:);		% Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;

Pxx = zeros(nfft,1); Pxx2 = zeros(nfft,1);
Pyy = zeros(nfft,1); Pyy2 = zeros(nfft,1);
Pxy = zeros(nfft,1); Pxy2 = zeros(nfft,1);
for i=1:k
    if strcmp(dflag,'none')
        xw = window.*x(index);
        yw = window.*y(index);
    elseif strcmp(dflag,'linear')
        xw = window.*detrend(x(index));
        yw = window.*detrend(y(index));
    else
        xw = window.*detrend(x(index),0);
        yw = window.*detrend(y(index),0);
    end
    index = index + (nwind - noverlap);
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft);
    Xx2 = abs(Xx).^2;
    Yy2 = abs(Yy).^2;
    Xy2 = Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pxx2 = Pxx2 + abs(Xx2).^2;
    Pyy = Pyy + Yy2;
    Pyy2 = Pyy2 + abs(Yy2).^2;
    Pxy = Pxy + Xy2;
    Pxy2 = Pxy2 + Xy2.*conj(Xy2);
end

% Select first half
if ~any(any(imag([x y])~=0)),   % if x and y are not complex
    if rem(nfft,2),    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end
    Pxx = Pxx(select);
    Pxx2 = Pxx2(select);
    Pyy = Pyy(select);
    Pyy2 = Pyy2(select);
    Pxy = Pxy(select);
    Pxy2 = Pxy2(select);
else
    select = 1:nfft;
end
%Coh = (abs(Pxy).^2)./(Pxx.*Pyy);             % coherence function estimate
Coh = Pxy./sqrt(Pxx.*Pyy);
freq_vector = (select - 1)'*Fs/nfft;

% set up output parameters
if (nargout == 2),
   Cxy = Coh;
   f = freq_vector;
elseif (nargout == 1),
   Cxy = Coh;
elseif (nargout == 0),   % do a plot
   newplot;
   plot(freq_vector,Coh), grid on
   xlabel('Frequency'), ylabel('Coherence Function Estimate');
end

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
function [results] = hyperPlmf(M, t, windowSize)
% HYPERPLMF Performs the PCA local matched filter (PLMF) target detection algorithm
%   Performs the PCA local matched filter (PLMF) target detection algorithm.
%
% Usage
%   [results] = hyperPlmf(M, target, windodwSize)
% Inputs
%   M - dd matrix of HSI data (m x n x p)
%   t - target of interest (p x 1)
% Outputs
%   results - vector of detector output (m x n)
%
% References
%   Sofa, Geva, Rotman. "Improved covariance matrices for point target detection in hyperspectral
% data."  IEEE International Conference on Microwaves, Communications, Antennas and Electronics
% Systems, 2009. COMCAS 2009.

% windowSize must be odd number
if ~mod(windowSize,2)
    error('windowSize must be an odd number.')
end

if (length(size(M)) ~= 3)
    error('M must be 3-dimensional matrix.')
end

[h,w,p] = size(M);
N = h*w;

% Remove mean from the target
M = hyperConvert2d(M);
u = mean(M.').';

[Mpca,V,lambda] = hyperPct(M,p);
t_pct = V.'*(t-u);

% Create map to get neighbors
map = 1:N;
map = reshape(map,h,w);

R_hat = hyperCov(Mpca);
G = inv(R_hat);

results = zeros(h,w);
s = floor(windowSize/2)+1;
for k=s:(h-s)
    for kk=s:(w-s)
        midIdx = map(k,kk);
        neighborhoodIdx = map((k-s+1):(k+s-1),(kk-s+1):(kk+s-1));

        Mlocal = M(:,neighborhoodIdx(:));
        [~,~,lambdaLocal] = hyperPct(Mlocal,p);

        y = Mpca(:,midIdx);
        results(k,kk) = sum((t_pct.*y)./max(lambdaLocal,lambda));
    end
end

author: davidkun

repo: HyperSpectralToolbox


text:
load -ascii heart

[N, m] = size(heart)

class = N

%rand('state',0); randn('state',0);
%heart = heart(:,randperm(m));

Napp = 150
Ntest = m-Napp

app  = heart(:,1:Napp);size(app)
test = heart(:,Napp+1:end);size(test)

unique(app(class,:))
unique(test(class,:))

ns = max(heart')
clear heart

% N, ns(class), Napp, Ntest, mean(ns),


text:
classdef PrincipalStressDirections < handle

    properties (GetAccess = public, SetAccess = private)
        principalStressDir
        principalStress
    end

    properties (Access = private)
        stress
    end

    methods (Access = public)

        function obj = PrincipalStressDirections(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.computeStressBase();
            obj.normalizePrincipalDirection();
            obj.transformBaseToElementalBase();
        end
    end

    methods (Access = private)

        function init(obj,cParams)
            obj.stress = cParams.stress;
        end

        function computeStressBase(obj)
            s = obj.stress;
            S = [s(1) s(3);s(3) s(2)];
            [V,D] = eig(S);
            obj.principalStressDir = V;
            obj.principalStress = D;
        end

    end

    methods (Access = private)

        function normalizePrincipalDirection(obj)
            V = obj.principalStressDir;
            V(:,1) = V(:,1)/norm(V(:,1));
            V(:,2) = V(:,2)/norm(V(:,2));
            obj.principalStressDir = V;
        end

        function transformBaseToElementalBase(obj)
            V = obj.principalStressDir;
            V(:,2) = V(:,2)*det(V);
            obj.principalStressDir = V;
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function results = CCR3_enterFunc(gamma, inputPath)
     %% ======================================================================
    %%STEP 1: load the data
    fprintf('start load the data...\n');
    load(inputPath);
%     TrainLabel(TrainLabel==0) = -1;

    if strcmp(class(TrainData),'cell')
        for x = 1:size(TrainData,2)
            TrainData{1,x} = TrainData{1,x}';
            tem = TrainLabel{1,x};
            TrainLabel{1,x}(tem==0) = -1;
        end
        traindata = TrainData;
        trainlabel = TrainLabel;
    else
        TrainData = TrainData';
        TrainLabel(TrainLabel==0) = -1;
        traindata{1,1} = TrainData;
        trainlabel{1,1} = TrainLabel;
    end
    if strcmp(class(TestData),'cell')
        tem = TestLabel{1,1};
        TestLabel{1,1}(tem==0) = -1;
        testdata = TestData;
        testlabel = TestLabel;
    else
        TestLabel(TestLabel==0) = -1;
        testdata{1,1} = TestData;
        testlabel{1,1} = TestLabel;
    end

     %% ======================================================================
    %%STEP 2: train the model
    results=CCR(traindata, testdata, trainlabel, testlabel, gamma);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_exp_to_ss(magnSigma2, lengthScale)
% CF_EXP_TO_SS - Exponential covariance functions to state space
%
% Syntax:
%   [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_exp_to_ss(magnSigma2, lengthScale)
%
% In:
%   magnSigma2  - Magnitude scale parameter (default: 1)
%   lengthScale - Length scale parameter (default: 1)
%
% Out:
%   F           - Feedback matrix
%   L           - Noise effect matrix
%   Qc          - Spectral density of white noise process w(t)
%   H           - Observation model matrix
%   Pinf        - Covariance of the stationary process
%   dF          - Derivatives of F w.r.t. parameters
%   dQc         - Derivatives of Qc w.r.t. parameters
%   dPinf       - Derivatives of Pinf w.r.t. parameters
%   params      - Input and output parameter information
%
% Description:
%   This function converts the exponential (Ornstein-Uhlenbeck) covariance
%   function to state space models. The covariance function
%   parametrization is as follows
%
%     k(tau) = magnSigma2 exp(-|tau|/lengthScale),
%
%   where magnSigma2 is the magnitude scale parameter, lengthScale the
%   distance scale parameter, and tau the time difference between states,
%   tau = t-t'.
%     This function takes the covariance function parameters as inputs and
%   outputs the corresponding state space model matrices. The state space
%   model is given as follows in terms of a stochastic differential
%   equation
%
%      df(t)/dt = F f(t) + L w(t),
%
%   where w(t) is a white noise process with spectral denisty Qc. The
%   observation model for discrete observation y_k of f(t_k) at step k,
%   is as follows
%
%      y_k = H f(t_k) + r_k, r_k ~ N(0, R),
%
%   where r_k is the Gaussian measurement noise with covariance R.
%     Pinf is the stationary covariance, where the value of Pinf(i,j),
%   is defined as follows
%
%      Pinf(i,j) = E[(f_i(t)-E[f_i(t)])(f_j(t)-E[f_j(t)])],
%
%   where f_i(t) is component i of state vector f(t).
%     Derivatives: All have same form. For example, dF has the following
%   form:
%
%       dF(:,:,1) = dF/d(magnSigma2 = input parameter_1),
%       dF(:,:,i) = dF/d(input parameter_i).
%
% References:
%
%   [1] Simo Sarkka, Arno Solin, Jouni Hartikainen (2013).
%       Spatiotemporal learning via infinite-dimensional Bayesian
%       filtering and smoothing. IEEE Signal Processing Magazine,
%       30(4):51-61.
%
%   [2] Jouni Hartikainen and Simo Sarkka (2010). Kalman filtering and
%       smoothing solutions to temporal Gaussian process regression
%       models. Proceedings of IEEE International Workshop on Machine
%       Learning for Signal Processing (MLSP).
%
% See also:
%   COV_EXP, SPEC_EXP
%
% Copyright:
%   2012-2014   Arno Solin
%   2013-2014   Jukka Koskenranta
%
%  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.

%% Apply defaults

  % Check if magnSigm2 is given
  if nargin < 1 || isempty(magnSigma2), magnSigma2 = 1; end

  % Check if lengthScale is given
  if nargin < 2 || isempty(lengthScale), lengthScale = 1; end

%% Form state space model

  % Feedback matrix
  F  = -1/lengthScale;

  % Noise effect matrix
  L  = 1;

  % Spectral density
  Qc = 2*magnSigma2/lengthScale;

  % Observation model
  H  = 1;

%% Stationary covariance

  % Calculate Pinf only if requsted
  if nargout > 4,
    Pinf = magnSigma2;
  end

%% Calculate derivatives

  % Calculate derivatives only if requested
  if nargout > 5

    % Derivative of F w.r.t. parameter magnSigma2
    dFmagnSigma2 = 0;

    % Derivative of F w.r.t parameter lengthScale
    dFlengthScale = 1/lengthScale^2;

    % Derivative of Qc w.r.t. parameter magnSigma2
    dQcmagnSigma2 = 2/lengthScale;

    % Derivative of Qc w.r.t. parameter lengthScale
    dQclengthScale = -2*magnSigma2/lengthScale^2;

    % Derivative of Pinf w.r.t. parameter magnSigma2
    dPinfmagnSigma2 = 1;

    % Derivative of Pinf w.r.t. parameter lengthScale
    dPinflengthScale = 0;

    % Stack all derivatives
    dF = zeros(1,1,2);
    dQc = zeros(1,1,2);
    dPinf = zeros(1,1,2);

    dF(:,:,1) = dFmagnSigma2;
    dF(:,:,2) = dFlengthScale;
    dQc(:,:,1) = dQcmagnSigma2;
    dQc(:,:,2) = dQclengthScale;
    dPinf(:,:,1) = dPinfmagnSigma2;
    dPinf(:,:,2) = dPinflengthScale;

  end

%% Return parameter names

  % Only return if requested
  if nargout > 8

    % Stationarity
    p.stationary = true;

    % Input parameter information
    p.in{1}.name = 'magnSigma2';  p.in{1}.default = 1; p.in{1}.opt = true;
    p.in{2}.name = 'lengthScale'; p.in{2}.default = 1; p.in{2}.opt = true;

    % Return parameter setup
    params = p;

  end

author: gpstuff-dev

repo: gpstuff

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


text:
function [rbm,dae,allZero,allOne] = ModelTraining(Mask,Dec,REAL)
% Training RBM and DAE

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Determine the size of hidden layers
    allZero = all(~Mask,1);
    allOne  = all(Mask,1);
    other   = ~allZero & ~allOne;
    K       = sum(mean(abs(Mask(:,other).*Dec(:,other))>1e-6,1)>rand(1,sum(other)));
    K       = min(max(K,1),size(Mask,1));

    %% Train RBM and DAE
    rbm = RBM(sum(other),K,10,1,0,0.5,0.1);
    rbm.train(Mask(:,other));
    if REAL
        dae = DAE(size(Dec,2),K,10,size(Dec,1),0.5,0.5,0.1);
        dae.train(Dec);
    else
        dae = [];
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function data = dcm2ml_Element(el, tag, varargin)
%"dcm2ml_element"
%   Convert a Java SimpleDicomElement object into a Matlab datatype.
%
%   TO DO: Define VRs labeled "Needs implementation."  All elements require
%   testing on various plans.
%
%   Discription for all the VR types
%     FD => 'double',
%     FL => 'float',
%     OB => 'int8u',
%     OF => 'float',
%     OW => 'int16u',
%     SL => 'int32s',
%     SS => 'int16s',
%     UL => 'int32u',
%     US => 'int16u',
%
%JRA 6/1/06
%   DK
%       Added support for multiple VR types.
%
%Usage:
%   data = dcm2ml_element(Java SimpleDicomElement)
%
% 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/>.

%Set cache buffer to true.
%{
buf = 1;

%Set character set to [], default.
cs = [];

%Get the tag value as a char array.
try
    tag = char(org.dcm4che2.util.TagUtils.toString(el.tag)); %removed el.tag
catch
    data = '';
    return;
end

%Get the VR, cast to ML char array.
vr = char(el.vr.toString);
%}
%modalityTag = '00080060';

vr = char (el.getVR(hex2dec(tag)));
% vr = cell(el.getVR(hex2dec(tag)));
% vr = vr{1}; % speedup

buf = 1;
cs = [];

switch upper(vr)
    case 'AE'
        %Needs implementation
        data = '';
    case 'AS'
        %Needs implementation
        data = '';
    case 'AT'
        data = dec2hex(el.getInts(buf));
    case {'CS', 'LO', 'SH', 'ST'}
        data = el.getStrings(cs, buf);
        %If more than one string, put in cell array.
        if numel(data) > 1
            data = cell(data);
        else
            data = char(data);
        end
    case 'DA'
        %Date string format: YYYYMMDD
        data = char(el.getString(cs, buf));
    case 'DS'
        data = el.getDoubles(buf);
    case 'DT'
        data = el.getDate(buf);

    case 'FL'
        %Needs implementation
        %wy
        %data =  float(el.getFloat(buf));
        data =  el.getFloats(buf);

    case 'FD'
        data = el.getDoubles(buf);
    case 'IS'
        data = el.getInts(buf);
    case 'LT'
        data = char(el.getString(cs, buf));
    case 'OB'
        %%%%%% Modified to import compressed data AI 02/06/17 %%%%%%%
        txSyntax = varargin{1};
        switch txSyntax
            case {'1.2.840.10008.1.2.4.50'
                    '1.2.840.10008.1.2.4.57'
                    '1.2.840.10008.1.2.4.70'
                    '1.2.840.10008.1.2.4.90'
                    '1.2.840.10008.1.2.4.91'
                    }
                %Decompress JPEG frame
                nElements = el.countItems;
                %if nElements>2
                %    warning(' dcm2ml_Element does not support multiple fragments');
                %    return
                %else
                    data = [];
                    for iFrag = 1:nElements
                    fragment = typecast(el.getFragment(iFrag-1),'uint8');
                    if ~isempty(fragment)
                    fileName = getTempName;
                    fid = fopen(fileName,'w');
                    if (fid < 0)
                        error('dcm2ml_element:Could not create temp file');
                    end
                    fwrite(fid,fragment,'uint8');
                    fclose(fid);
                    %tmp = onCleanup(@() delete(fileName));
                    %dataTmpV = imread(fileName).';
                    dataTmpV = imread(fileName);
                    if ndims(dataTmpV) == 3
                        dataTmpV = rgb2gray(dataTmpV); % temp for SM modality
                    end
                    dataTmpV = permute(dataTmpV,[2,1,3]);
                    data = [data;dataTmpV(:)];
                    end
                    end
                    delete(fileName)
                %end

            case '1.2.840.10008.1.2.5'

                % To do: Decompress RLE frame

            case {'1.2.840.10008.1.2'  %Implicit VR Little Endian (default)
                    '1.2.840.10008.1.2.1' %Explicit VR Little Endian
                    '1.2.840.10008.1.2.2' %Explicit VR Big Endian
                    '1.2.840.113619.5.2' %Implicit VR Big Endian (GE pvt)
                    '1.3.46.670589.33.1.4.1'}%Explicit VR Little Endian (Philips pvt)

                data = el.getBytes;

            otherwise
                error('dc2ml_Element : Encoding not supported');
        end
        %%%%%%%%%%%%%%%%%%% End Modified %%%%%%%%%%%%%%%%%%%%%

    case 'OF'
        data = el.getFloats(buf);
    case 'OW'
        %OW contains 16 bit words.  Conversion of this data into meaningful
        %values is the responsibility of the calling function.

        %wy it should be uint16 or int16, which depends on the value in data
        %representation fields, but the data conversion b/w matlab and java is int32.

        %data = uint16(el.getInts(buf));
        data = el.getInts(buf);
    case 'PN'
        nameObj = org.dcm4che3.data.PersonName(el.getString(cs, buf));

        %The # in get(#) as defined by dcm4che2, PersonName class.
        data.FamilyName = char(nameObj.get(0));
        data.GivenName  = char(nameObj.get(1));
        data.MiddleName = char(nameObj.get(2));
        data.NamePrefix = char(nameObj.get(3));
        data.NameSuffix = char(nameObj.get(4));
    case 'SL'
        %
         data = el.getInt(buf);
    case 'SQ'
        nElements = el.countItems;
        data = [];
        for i=0:nElements-1
            data.(['Item_' num2str(i+1)]) = dcm2ml_Object(el.getDicomObject(i));
        end
    case 'SS'
        %Needs implementation
        data = '';
    case 'TM'
        %Time string format: HHMMSS.ss where "ss" is fraction of a second.
        data = char(el.getString(cs, buf));
    case 'UI'
        data = char(el.getString(cs, buf));
    case 'UL'
        data = el.getInts(buf);
    case 'UN'
        data = el.getBytes;
    case 'US'
        data = el.getInt(buf);
    case 'UT'
        %Needs implementation
        data = '';
    otherwise
        error('Unrecognized VR type.'); %%Consider more gracious exit.
end

%DEBUGGING: remove this once all DICOM VRs are implemented and fully
%tested.  Until then, reaching this point in the code indicates that a VR
%MUST be defined for proper functioning of a called module.
if ~exist('data', 'var');
    disp(['DEBUGGING: ' vr ' is not defined.  Implement it in dcm2ml_Element.m']);
    data = '';
else
    %Handle empty data situations -- this needs to be tailored to individual
    %VR values if matching MATLAB's dicominfo function output is desired.
    if isempty(data);
        data = '';
    end
end

end % end of function

%% Function to generate temp filename
function tempfile = getTempName()
nFiles = 1;
persistent tempfiles counter;
if (isempty(tempfiles))
    for p = 1:nFiles
        tempfiles{p} = tempname;
    end

    counter = 0;
end
counter = rem(counter, nFiles) + 1;
tempfile = tempfiles{counter};

end % end of function

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function openCalcParamWindowGK(obj)
		%Here a window for the calculation parameters should be created
		%I suggest using the inputsdlg function

		Title = 'Gardner-Knopoff Decluster ';
		Prompt={'Window Method:', 'GK_Method';...
			'Force MainShock:', 'GK_SetMain'};

		LabelList2 =  {	'Gardener & Knopoff (table), 1974';...
				'Gardener & Knopoff (formula), 1974'; ...
				'Gruenthal pers. communication'; ...
				'Urhammer, 1986'};

		%Method
		Formats(1,1).type='list';
		Formats(1,1).style='popupmenu';
		Formats(1,1).items=LabelList2;

		%SetMain
		Formats(2,1).type='check';

		%%%% SETTING DIALOG OPTIONS
		Options.WindowStyle = 'modal';
		Options.Resize = 'on';
		Options.Interpreter = 'tex';
		Options.ApplyButton = 'off';

		%open the dialog window
		[NewParameter,Canceled] = inputsdlg(Prompt,Title,Formats,obj.CalcParameter,Options);

		obj.CalcParameter=NewParameter;

		if Canceled==1
			obj.StartCalc=false;
		else
			obj.StartCalc=true;
		end

end

author: CelsoReyes


text:
function chebcoeffs = legcoeffs2chebcoeffs(legcoeffs)
%CHEB2LEG  Convert Legendre coefficients to Chebyshev coefficients.
%   C_CHEB = LEG2CHEB(C_LEG) converts the vector C_LEG of Legendre coefficients
%   to a vector C_CHEB of Chebyshev coefficients such that
%       C_CHEB(1)*T0 + ... + C_CHEB(N)*T{N-1} = ...
%           C_LEG(N)*P0 + ... + C_LEG(1)*P{N-1},
%   where P{k} is the degree k Legendre polynomial normalized so that max(|P{k}|
%   = 1.

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

% This command is a wrapper for leg2cheb.
chebcoeffs = leg2cheb(legcoeffs);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Ycon,ln] = concatenate_files(files,filename,output_type)

% concatenate a list of files into a single file
% INPUTS
% files:        list of a files as a struct array. Name of each file is in files(i).name
% filename:     output filename
% type:         output type ('mat' array loaded in RAM (default), 'hdf5' and 'tif' are currently supported)

% OUTPUTS
% Ycon:         concantenated file as an array if output_type == 'mat'
% ln:           length of each file (in frames)

if ~exist('output_type','var')
    output_type = 'mat';
end

if ~exist('filename','var') || isempty(filename)
    filename = ['concatenated_file.',output_type];
end

numFiles = length(files);
ln = zeros(numFiles,1);

Y1 = read_file(fullfile(files(1).folder, files(1).name));
sizY = size(Y1);
T1 = sizY(end);
data_type = class(Y1);
nd = max(ndims(Y1)-1,2);
sizY = sizY(1:nd);

switch lower(output_type)
    case 'mat'
        Ycon = cast([],data_type); %zeros([sizY,T1],data_type);
    case {'hdf5','h5'}
        Ycon = filename;
        h5create(filename,'/mov',[sizY,Inf],'Chunksize',[sizY,1],'Datatype',data_type);
    case {'tif','tiff'}
        Ycon = filename;
        opts_tiff.append = true;
        opts_tiff.big = true;
    otherwise
        error('This filetype is currently not supported')
end

for i = 1:numFiles
    Y_temp = read_file(fullfile(files(i).folder, files(i).name)); %cat(nd+1,Ycon,Y_temp);
    ln(i) = size(Y_temp,ndims(Y_temp));

    switch lower(output_type)
        case 'mat'
            Ycon = cat(nd+1,Ycon,Y_temp);
        case {'hdf5','h5'}
            h5write(filename,'/mov',Y_temp,[ones(1,nd),sum(ln(1:i-1))+1],size(Y_temp));
        case {'tif','tiff'}
            saveastiff(cast(Y_temp,data_type),filename,opts_tiff);
    end
end

author: flatironinstitute

repo: NoRMCorre


text:
%  Figure 3.40      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.40
fh=@(ki,k) 6+3*k-ki;
ezplot(fh)
hold on;
f=@(ki,k) ki;
ezplot(f);

% add shading
fillvert = [0, -2
	        0,  7
			7, 7
			7, 0];
fcolor = [ .1 .1 .1 .1];
patch(fillvert(:,1)', fillvert(:,2)', fcolor)

xlabel('K_I');
ylabel('K');
title('Allowable region for stability');
nicegrid;
hold off

author: Sable

repo: mcbench-benchmarks


text:
function J = lotkaObjectiveFCN(u,x,Ts,N,xref,u0,p,Q,R,Ru)
%% Cost function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%   xref:   state references, constant from time k+1 to k+N
%   u0:     previous controller output at time k-1
%
% Output:
%   J:      objective function cost
%

%% Nonlinear MPC design parameters

%% Cost Calculation
% Set initial plant states, controller output and cost
xk = x;
uk = u(1);
J = 0;
% Loop through each prediction step
for ct=1:N
    % Obtain plant state at next prediction step
    xk1 = rk4u(@sparseGalerkinControl_Discrete,xk,uk,Ts,1,[],p);

    % Accumulate state tracking cost from x(k+1) to x(k+N)
    J = J + (xk1-xref)'*Q*(xk1-xref);
    % Accumulate MV rate of change cost from u(k) to u(k+N-1)
    if ct==1
        J = J + (uk-u0)'*R*(uk-u0) + uk'*Ru*uk;
    else
        J = J + (uk-u(ct-1))'*R*(uk-u(ct-1)) + uk'*Ru*uk;
    end
    % Update xk and uk for the next prediction step
    xk = xk1;
    if ct<N
        uk = u(ct+1);
    end
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function vUnique=uniqueVec(v)
%%UNIQUEVEC This function takes an mXn matrix of n m-dimensional vectors
%           and returns only the vectors that are unique. This is currently
%           suited for a small to moderate  numbers of vectors, because a
%           brute-force comparison algorithm is used (rather than a
%           sorting-based approach that might lessen the number of
%           comparisons).
%
%INPUTS: v An mXn matrix of n vectors.
%
%OUTPUTS: vUnique An mXnumUnique matrix of the unique vectors in v.
%
%EXAMPLE:
% v=[1, 1, 1, 1, 1, 12;
%    1, 2, 1, 2, 2, 24;
%    1, 3, 1, 4, 3, 36]
% uniqueVec(v)
%
%February 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(v))
    vUnique=v;
    return;
end

numRows=size(v,1);
numCols=size(v,2);

vUnique=zeros(numRows,numCols);

vUnique(:,1)=v(:,1);
numUnique=1;
for k=2:numCols
    isUnique=true;
    for kPrev=1:numUnique
        if(all(v(:,k)==vUnique(:,kPrev)))
            isUnique=false;
            break;
        end
    end

    if(isUnique)
        numUnique=numUnique+1;
        vUnique(:,numUnique)=v(:,k);
    end
end

%Size to fit
vUnique=vUnique(:,1:numUnique);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [month, day, year] = gdate (jdate)

% convert Julian date to Gregorian (calendar) date

% input

%  jdate = julian day

% output

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

%  note: day may include fractional part

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jd = jdate;

z = fix(jd + .5);
fday = jd + .5 - z;

if (fday < 0)
   fday = fday + 1;
   z = z - 1;
end

if (z < 2299161)
   a = z;
else
   alpha = floor((z - 1867216.25) / 36524.25);
   a = z + 1 + alpha - floor(alpha / 4);
end

b = a + 1524;
c = fix((b - 122.1) / 365.25);
d = fix(365.25 * c);
e = fix((b - d) / 30.6001);
day = b - d - fix(30.6001 * e) + fday;

if (e < 14)
   month = e - 1;
else
   month = e - 13;
end

if (month > 2)
   year = c - 4716;
else
   year = c - 4715;
end

author: goGPS-Project


text:
%% Machine Learning Online Class - Exercise 3 | Part 1: One-vs-all

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     lrCostFunction.m (logistic regression cost function)
%     oneVsAll.m
%     predictOneVsAll.m
%     predict.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Setup the parameters you will use for this part of the exercise
input_layer_size  = 400;  % 20x20 Input Images of Digits
num_labels = 10;          % 10 labels, from 1 to 10   ~k
                          % (note that we have mapped "0" to label 10)

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  You will be working with a dataset that contains handwritten digits.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('ex3data1.mat'); % training data stored in arrays X, y
m = size(X, 1);

% Randomly select 100 data points to display
%rand_indices = randperm(m);
%sel = X(rand_indices(1:100), :);

%displayData(sel);

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ============ Part 2: Vectorize Logistic Regression ============
%  In this part of the exercise, you will reuse your logistic regression
%  code from the last exercise. You task here is to make sure that your
%  regularized logistic regression implementation is vectorized. After
%  that, you will implement one-vs-all classification for the handwritten
%  digit dataset.
%

fprintf('\nTraining One-vs-All Logistic Regression...\n')

lambda = 0.1;
[all_theta] = oneVsAll(X, y, num_labels, lambda);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Predict for One-Vs-All ================
%  After ...
pred = predictOneVsAll(all_theta, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

author: scruel

repo: Notes-ML-AndrewNg


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [y,dy] = translation3D(w,x,varargin)
%
% computes y = Q*w = and the derivative wrt. w.
% x = reshape(x,[],2);
% Q = [1 0 0; 0 1 0;0 0 1];  dy = Q:
% if no arguments are given, the parameters for the identity map are returned.
%
% see also transformations/contents.m, trafo.m
%==============================================================================

function [y,dy] = translation3D(w,x,varargin)

% the persistent variable Q stores the matrix
% Q(x) = kron( I_2 , [x(:,1),x(:,2),1] );
persistent Q

if nargin == 0,
  runMinimalExample;
  return;
else
  y = mfilename('fullfile');
  dy = [0;0;0];         % parameterization of identity
  if ischar(w),  Q  = []; w = []; end; % reset Q
  if isempty(w), return;          end;
end;

if isempty(w) || (size(Q,1) ~= numel(x)),
  n = length(x)/3;
  Q = sparse(kron(speye(3),ones(n,1)));
  if nargout == 0, return; end;
end;
y  = x + Q*w;
dy = Q;

%------------------------------------------------------------------------------
function runMinimalExample
help(mfilename);
fprintf('%s: minimal example\n',mfilename)

omega = [0,10,0,8,0,9]; m = [8,7,6];
c = (omega(2:2:end)-omega(1:2:end))'/2;
w = c/10
x = getCellCenteredGrid(omega,m);
z = feval(mfilename,w(:),x);
FAIRfigure(1); clf;
plotGrid(x,omega,m,'color','r'); axis image; hold on;
plotGrid(z,omega,m,'color','b'); axis image; hold off;

fctn = @(w) feval(mfilename,w,x);
w = w + randn(size(w));
checkDerivative(fctn,w,'fig',2);
%==============================================================================


text:
%TEST 04: Algebraic Reconstruction Techniques
clear
clc
fprintf('TEST 04: Algebraic Reconstruction Techniques\rCompare ART and SART\r');
angles = 0:6:179;
%Phantom
im = imtest('s4',64);
siz = size(im);
figure('name','Original Image: Phantom')
imshow(im)
[W, p, ~, ~] = buildWeightMatrix(im,angles);
%Iteration 100
fprintf('ART method\niteration 100:\r')
im_rec1 = tomo_recon_myart(W,p,siz,100);
im_rec1 = uint8(imscale(im_rec1));
figure('name','ART')
imshow(im_rec1);
fprintf('SART method\niteration 100:\r')
im_rec2 = tomo_recon_sart(W,p,siz,100);
im_rec2 = uint8(imscale(im_rec2));
figure('name','SART')
imshow(im_rec2);
%Iteration 1000
fprintf('ART method\niteration 1000:\r')
im_rec1 = tomo_recon_myart(W,p,siz,1000);
im_rec1 = uint8(imscale(im_rec1));
figure('name','ART')
imshow(im_rec1);
fprintf('SART method\niteration 1000:\r')
%Code efficiency is too low!
%im_rec2 = tomo_recon_mysart(W,p,siz,1000);
%figure('name','ART')
%imshow(im_rec2);


text:
function SVF = normal
% the normal vector field
%

SVF = S2VectorFieldHandle(@(v) v);

end

author: mtex-toolbox

repo: mtex

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


text:
function tagS = content_identification_module_tags
%"content_identification_module_tags"
%   Return the tags used to represent content as specified by C.10.9 in
%   PS3.3 of 2006 DICOM specification. Tags are returned in a struct array
%   with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 06/06/06
%YWU Modified 03/01/08
%
%Usage:
%   tagS = patient_module_tags
%
% 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/>.

tagS = content_identification_macro_tags;

%
% %Initialize the tagS structure.
% tagS = struct('tag', {}, 'type', {}, 'children', {});
%
% %Create an empty tagS template for sequence creation.
% template = tagS;
%
% %Add tags based on PS3.3 attribute lists.
%
% %Instance Number
% tagS(end+1) = struct('tag', ['00200013'], 'type', ['1'], 'children', []);
%
% %Content Label
% tagS(end+1) = struct('tag', ['00700080'], 'type', ['1'], 'children', []);
%
% %Content Description
% tagS(end+1) = struct('tag', ['00700081'], 'type', ['2'], 'children', []);
%
% %Content Creator's Name
% tagS(end+1) = struct('tag', ['00700084'], 'type', ['2'], 'children', []);
%


text:
 function MyPerecptronExample
% Rosenblatt's Perecptron
% References :
%      Neural Networks for Pattern Recognition By C. Bishop
%      An Introduction to Support Vector Machines By Cristianini
%      Introduction to Machine Learning By Alpayd?n
%      Neural Networks. A Comprehensive Foundation By Haykin
% Code by By Ibraheem Al-Dhamari 2010

clc
%============================================
% Generate 2 dimensions linear separable data
%===========================================

% You may change the size of the data from here or input your own data
 %   note that the drawing is for two dimensions only, hence you need to
 %   modify the code for different data.
mydata = rand(500,2);
% Separate the data into two classes
acceptindex = abs(mydata(:,1)-mydata(:,2))>0.012;
acceptindex = abs(mydata(:,1)-mydata(:,2))>0.012;
mydata = mydata(acceptindex,:); % data
myclasses = mydata(:,1)>mydata(:,2); % labels
[m n]=size(mydata);
%training data
 x=mydata(1:400,:);   y=myclasses(1:400);
% test data
xt=mydata(401:m,:); yt=myclasses(401:m);
%=====================================
% Train the perceptron
%=====================================
[w,b,pass] = PerecptronTrn(x,y);
Iterations=pass

%=====================================
% Test
%=====================================
e=PerecptronTst(xt,yt,w,b);
disp(['Test_Errors=' num2str(e) '     Test Data Size= ' num2str(m-400)])

%=====================================
% Draw the result (sparating hyperplane)
%=====================================
 l=y==0;
 hold on
 plot(x(l,1),x(l,2),'k.' );
 plot(x(~l,1),x(~l,2),'g.');
 [l,p]=size(x);
 plot([0,1],[0,1],'r-')
 axis([0 1 0 1]), axis square, grid on
 drawnow


text:
function test_bug2994

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqstatistics ft_freqgrandaverage ft_selectdata

% Bug 2994 has been reported by Niels. Apparently once upon a time it was
% possible to call ft_freqstatistics with ft_freqgrandaveraged data in the
% input, where the dimord is 'subj_freq_time' (NOTE the lacking chan
% dimord) with the data containing apparently an average across channels.
% The input data contains a label field. ft_freqstatistics crashes when the
% low level function tries to cluster across space.

% Let's start this form the beginning, create a multisubject
% freq-structure, with data averaged across channels, and see whether it
% pulls through

tmp.freq   = [1 2 3];
tmp.time   = [1 2 3 4 5];
tmp.label  = {'a';'b'};
tmp.dimord = 'chan_freq_time';

freq = cell(1,10);
freq0 = cell(1,10);
for k = 1:10
  tmp.powspctrm = rand(2,3,5);
  freq{k} = tmp;
  tmp.powspctrm(:) = 0;
  freq0{k} = tmp;
end

%% ------------------------------------------------------------------------
% here, the following is done:
% -ft_freqgrandaverage (keepindividual = 'yes')
% -ft_selectdata (avgoverchan = 'yes')
% -ft_freqstatistics

cfg1 = [];
cfg1.keepindividual = 'yes';

cfg2 = [];
cfg2.avgoverchan = 'yes';

cfg3 = [];
cfg3.method    = 'montecarlo';
cfg3.statistic = 'ft_statfun_depsamplesT';
cfg3.design    = [ones(1,10) ones(1,10)*2;1:10 1:10];
cfg3.numrandomization = 100;
cfg3.correctm  = 'cluster';
cfg3.ivar      = 1;
cfg3.uvar      = 2;

freqall  = ft_selectdata(cfg2, ft_freqgrandaverage(cfg1, freq{:}));
freqall0 = ft_selectdata(cfg2, ft_freqgrandaverage(cfg1, freq0{:}));
assert(isequal(size(freqall.powspctrm),[10 1 3 5])); % there should be a singleton channel dimension
stat1    = ft_freqstatistics(cfg3, freqall, freqall0);

%% ------------------------------------------------------------------------
% here, the following is done:
% -ft_freqgrandaverage (keepindividual = 'yes')
% -ft_freqstatistics (avgoverchan = 'yes')

cfg1 = [];
cfg1.keepindividual = 'yes';

cfg3 = [];
cfg3.method    = 'montecarlo';
cfg3.statistic = 'ft_statfun_depsamplesT';
cfg3.avgoverchan = 'yes';
cfg3.design    = [ones(1,10) ones(1,10)*2;1:10 1:10];
cfg3.numrandomization = 100;
cfg3.correctm  = 'cluster';
cfg3.ivar      = 1;
cfg3.uvar      = 2;

freqall  = ft_freqgrandaverage(cfg1, freq{:});
freqall0 = ft_freqgrandaverage(cfg1, freq0{:});
assert(isequal(size(freqall.powspctrm),[10 2 3 5])); % there should be two channels
stat2    = ft_freqstatistics(cfg3, freqall, freqall0);

%% ------------------------------------------------------------------------
% here, the following is done:
% - NO ft_freqgrandaverage
% -ft_freqstatistics (avgoverchan = 'yes')

cfg3 = [];
cfg3.method    = 'montecarlo';
cfg3.statistic = 'ft_statfun_depsamplesT';
cfg3.avgoverchan = 'yes';
cfg3.design    = [ones(1,10) ones(1,10)*2;1:10 1:10];
cfg3.numrandomization = 100;
cfg3.correctm  = 'cluster';
cfg3.ivar      = 1;
cfg3.uvar      = 2;
cfg3.avgoverchan = 'yes';

stat3    = ft_freqstatistics(cfg3, freq{:}, freq0{:});

%% ------------------------------------------------------------------------
% here, the following is done:
% -ft_appendfreq
% -ft_freqstatistics (avgoverchan = 'yes')

cfg1 = [];
cfg1.parameter = 'powspctrm';
cfg1.appenddim = 'rpt';

cfg3 = [];
cfg3.method    = 'montecarlo';
cfg3.statistic = 'ft_statfun_depsamplesT';
cfg3.avgoverchan = 'yes';
cfg3.design    = [ones(1,10) ones(1,10)*2;1:10 1:10];
cfg3.numrandomization = 100;
cfg3.correctm  = 'cluster';
cfg3.ivar      = 1;
cfg3.uvar      = 2;
cfg3.avgoverchan = 'yes';

freqall  = ft_appendfreq(cfg1, freq{:});
freqall0 = ft_appendfreq(cfg1, freq0{:});

stat4    = ft_freqstatistics(cfg3, freqall, freqall0);

assert(isequal(stat1.stat,stat2.stat));
assert(isequal(stat1.stat,stat3.stat));
assert(isequal(stat1.stat,stat4.stat));
assert(isequal(stat2.stat,stat3.stat));
assert(isequal(stat2.stat,stat4.stat));
assert(isequal(stat3.stat,stat4.stat));

% the conclusion so far is that with the current state of the code, all
% seems to work fine.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Y = acsch(X)
    % Symbolic inverse hyperbolic cosecant.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['ArcCsch[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
classdef PTKWavefront < handle
    % PTKWavefront. A data structure representing a segmented airway tree
    %
    %     PTKWavefront is used as part of the airway region growing process in
    %     PTKAirwayRegionGrowingWithExplosionControl.
    %
    %     A root PTKTreeSegment is returned by
    %     PTKAirwayRegionGrowingWithExplosionControl. From this you
    %     can extract and analyse the resulting airway tree.
    %
    %     PTKTreeSegment is used in the construction and storage of
    %     the airway trees. A PTKTreeSegment stores an individual
    %     segment of the centreline tree, with references to the parent and child
    %     PTKTreeSegments, so that it is possible to reconstruct the entire
    %     tree from a single segment.
    %
    %     The way the airway voxels are stored in each segment is as follows:
    %         The wavefront only exists during the region growing. It is a thick
    %         layer of voxels which exists at the end of each segment
    %         which is currently growing. If the wavefront forms more than one
    %         connected component, the segment is ended and new child segments
    %         formed from the wavefront components. When new voxels are added to
    %         the front of the wavefront, old voxels from the back of the
    %         wavefront are pushed out into the Pending voxels of the segment.
    %         Pending voxels are 'accepted' or 'rejected' according to whether
    %         the heuristics have determined an explosion has occurred.
    %
    %     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 (SetAccess = private)
        CurrentBranch = []
    end

    properties (Access = private)

        % The wavefront is a thick layer of voxels which is used to detect
        % and process bifurcations in the airway tree before these voxels are
        % added to the segement's list of pending indices
        WavefrontVoxelIndices

        % Additional voxels that were not originally part of the segment but
        % were added later after a morphological closing operation
        ClosedPoints

        NumberOfVoxelsSkipped = 0

        WavefrontSize
        PermittedVoxelSkips = 0

        MinimumDistanceBeforeBifurcatingMm
        MinimumChildDistanceBeforeBifurcatingMm = 5

        FirstSegmentWavefrontSizeMm = 10
        ChildWavefrontSizeMm = 5
        VoxelSizeMm
        MinCoords
        MaxCoords

        % Generations greater than this are automatically terminated
        MaximumNumberOfGenerations

        MinimumNumberOfPointsThresholdMm3 = 6

        ExplosionMultiplier = 7

    end

    methods
        function obj = PTKWavefront(segment_parent, min_distance_before_bifurcating_mm, voxel_size_mm, maximum_generations, explosion_multiplier)
            if nargin > 0
                obj.WavefrontVoxelIndices = int32([]);

                obj.MinimumDistanceBeforeBifurcatingMm = min_distance_before_bifurcating_mm;
                max_voxel_size_mm = max(voxel_size_mm);
                obj.VoxelSizeMm = voxel_size_mm;
                obj.MaximumNumberOfGenerations = maximum_generations;
                obj.ExplosionMultiplier = explosion_multiplier;

                voxel_volume = voxel_size_mm(1)*voxel_size_mm(2)*voxel_size_mm(3);
                min_number_of_points_threshold = max(3, round(obj.MinimumNumberOfPointsThresholdMm3/voxel_volume));

                if ~isempty(segment_parent)
                    obj.WavefrontSize = ceil(obj.ChildWavefrontSizeMm/max_voxel_size_mm);
                    obj.CurrentBranch = PTKTreeSegment(segment_parent, min_number_of_points_threshold, explosion_multiplier);
                else
                    obj.WavefrontSize = ceil(obj.FirstSegmentWavefrontSizeMm/max_voxel_size_mm);
                    obj.CurrentBranch = PTKTreeSegment([], min_number_of_points_threshold, explosion_multiplier);
                end
            end
        end

        % Returns the very front layer of voxels at the wavefront
        function frontmost_points = GetFrontmostWavefrontVoxels(obj)
            frontmost_points = obj.WavefrontVoxelIndices{end};
        end

        % Returns the wavefront for this segment, which includes voxels that may
        % be separated into child segments
        function wavefront_voxels = GetWavefrontVoxels(obj)
            wavefront_voxels = obj.ConcatenateVoxels(obj.WavefrontVoxelIndices);
        end

        % Add new voxels to this segment, and returns a list of all segments
        % that require further processing (including this one, and any child
        % segments which have been created as a result of bifurcations)
        function segments_to_do = AddNewVoxelsAndGetNewSegments(obj, indices_of_new_points, image_size, reporting)

            % Check that indices are unique
            if (numel(indices_of_new_points) ~= numel(unique(indices_of_new_points)))
                reporting.Error('PTKWavefront:Duplicates', 'Algorithm error - some points have been duplicated');
            end

            if any(ismember(indices_of_new_points, obj.CurrentBranch.GetAcceptedVoxels))
                reporting.Error('PTKWavefront:Duplicates', 'Algorithm error - some points have been duplicated');
            end

            % First we move voxels at the rear of the wavefront into the
            % PendingVoxels
            if ~isempty(obj.WavefrontVoxelIndices)
                while length(obj.WavefrontVoxelIndices) > obj.WavefrontSize
                    obj.MoveVoxelsFromRearOfWavefrontToPendingVoxels;
                end
            end

            % Next add the new points to the front of the wavefront
            obj.WavefrontVoxelIndices{end + 1} = indices_of_new_points;

            % If an explosion has been detected then do not continue
            if obj.CurrentBranch.MarkedExplosion
                obj.MoveAllWavefrontVoxelsToPendingVoxels;
%                 obj.DeleteSegmentIfNoAcceptedVoxels;
                segments_to_do = PTKWavefront.empty; % This segment has been terminated
                return
            end

            obj.AdjustMaxAndMinForVoxels(indices_of_new_points, image_size);

            % Do not allow the segment to bifurcate until it is above a minimum
            % length
            if ~obj.MinimumLengthPassed
                segments_to_do = obj; % This segment is to continue growing
                return
            end

            % Do not allow the segment to bifurcate until it is above a minimum size
            if ~obj.WavefrontIsMinimumSize
                segments_to_do = obj; % This segment is to continue growing
                return
            end

            % Determine whether to continue growing the current segment, or to
            % split it into a new set of child segments

            % Find connected components from the wavefront (which is several voxels thick)
            [offset, reduced_image, reduced_image_size] = MimImageCoordinateUtilities.GetMinimalImageForIndices(int32(obj.GetWavefrontVoxels)', image_size);
            wavefront_connected_components = bwconncomp(reduced_image, 26);
            number_of_components = wavefront_connected_components.NumObjects;

            % If there is only one component, it will be growing, so there is no
            % need to do any further component analysis
            if number_of_components == 1
                segments_to_do = obj;
                return;
            end

            segments_to_do = PTKWavefront.empty;
            growing_branches = [];
            points_by_branches = [];

            % Iterate over the components and separate the wavefront voxels of
            % the current segment into a new branch for each growing component
            for component_number = 1 : wavefront_connected_components.NumObjects

                % Get voxel list, and adjust the indices to match those for the full image
                indices_of_component_points = wavefront_connected_components.PixelIdxList{component_number};
                indices_of_component_points = MimImageCoordinateUtilities.OffsetIndices(int32(indices_of_component_points), offset, reduced_image_size, image_size);
                points_by_branches{component_number} = indices_of_component_points;

                still_growing = obj.IsThisComponentStillGrowing(indices_of_component_points);
                if (still_growing)
                    growing_branches(end + 1) = component_number;
                end
            end

            if length(growing_branches) < 1
                reporting.Error('PTKWavefront:NoGrowingBranches', 'Algorithm error - no growing branches');
            end

            if length(growing_branches) == 1
                segments_to_do = obj;
                return
            end

            if length(growing_branches) > 1

                if ~isempty(obj.MaximumNumberOfGenerations)
                    % If the maximum permitted number of generations is exceeded
                    % then terminate this segment. This will discard any
                    % remaining wavefront voxels and mark the segment as
                    % incomplete (unless it is marked as exploded)
                    if obj.CurrentBranch.GenerationNumber >= obj.MaximumNumberOfGenerations
                        obj.CurrentBranch.EarlyTerminateBranch;
                        return;
                    end
                end

                for index = 1 : length(growing_branches)
                    segments_to_do(end + 1) = obj.SpawnChildFromWavefrontVoxels(points_by_branches{growing_branches(index)});
                end

                % If the branch has divided, there may be some unaccepted points
                % left over
                obj.CompleteThisSegment;

                if isempty(obj.CurrentBranch.GetAcceptedVoxels)
                    reporting.ShowWarning('PTKWavefront:EmptyBranch', 'Algorithm error - no points in final branch');
                end

            end

        end

        function CompleteThisSegment(obj)
            obj.MoveAllWavefrontVoxelsToPendingVoxels;
            obj.CurrentBranch.CompleteThisSegment
        end
    end

    methods (Access = private)

        function MoveAllWavefrontVoxelsToPendingVoxels(obj)
            while ~isempty(obj.WavefrontVoxelIndices)
                obj.MoveVoxelsFromRearOfWavefrontToPendingVoxels;
            end
        end

        function MoveVoxelsFromRearOfWavefrontToPendingVoxels(obj)
            % The wavefront may be empty after voxels have been divided
            % amongst child branches
            if ~isempty(obj.WavefrontVoxelIndices{1})
                obj.CurrentBranch.AddPendingVoxels(obj.WavefrontVoxelIndices{1});
            end
            obj.WavefrontVoxelIndices(1) = [];
        end

        function AdjustMaxAndMinForVoxels(obj, voxel_indices, image_size)
            [x, y, z] = ind2sub(image_size, voxel_indices);
            mins = [min(x), min(y), min(z)];
            maxs = [max(x), max(y), max(z)];
            if isempty(obj.MinCoords)
                obj.MinCoords = mins;
                obj.MaxCoords = maxs;
            else
                obj.MinCoords = min(mins, obj.MinCoords);
                obj.MaxCoords = max(maxs, obj.MaxCoords);
            end
        end

        function is_minimum_size = WavefrontIsMinimumSize(obj)
            is_minimum_size = (length(obj.WavefrontVoxelIndices) >= obj.WavefrontSize);
        end

        function passed_minimum_lengths = MinimumLengthPassed(obj)
            lengths = obj.MaxCoords - obj.MinCoords;
            lengths = single(lengths).*obj.VoxelSizeMm;
            max_length = max(lengths);

            passed_minimum_lengths = max_length >= obj.MinimumDistanceBeforeBifurcatingMm;
        end

        function new_segment = SpawnChildFromWavefrontVoxels(obj, voxel_indices)
            wavefront_voxels = [];
            for index = 1 : length(obj.WavefrontVoxelIndices)
                wavefront_voxels{index} = intersect(int32(voxel_indices), obj.WavefrontVoxelIndices{index});
                obj.WavefrontVoxelIndices{index} = setxor(wavefront_voxels{index}, obj.WavefrontVoxelIndices{index});
            end
            new_segment = PTKWavefront(obj.CurrentBranch, obj.MinimumChildDistanceBeforeBifurcatingMm, obj.VoxelSizeMm, obj.MaximumNumberOfGenerations, obj.ExplosionMultiplier);
            new_segment.WavefrontVoxelIndices = wavefront_voxels;
        end

        function still_growing = IsThisComponentStillGrowing(obj, voxel_indices)
            wavefront_voxels_end = intersect(int32(voxel_indices), obj.WavefrontVoxelIndices{end});
            still_growing = ~isempty(wavefront_voxels_end);
        end

    end

    methods (Static, Access = private)
        function concatenated_voxels = ConcatenateVoxels(voxels)
            concatenated_voxels = [];
            number_layers = length(voxels);
            for index = 1 : number_layers
                next_voxels = voxels{index};
                if ~isempty(next_voxels)
                    concatenated_voxels = cat(1, concatenated_voxels, next_voxels);
                end
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function PBI = CalPBI(PopObj,W,Region,Z,Sub)
% Calculate the PBI value between each solution and its associated weight
% vector

%------------------------------- 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".
%--------------------------------------------------------------------------

    Z   = repmat(Z,sum(Sub),1);
    NormW = sqrt(sum(W(Region(Sub),:).^2,2));
    d1  = abs(sum((PopObj(Sub,:)-Z).*W(Region(Sub),:),2))./NormW;
    d2  = sqrt(sum((PopObj(Sub,:)-(Z+W(Region(Sub),:).*repmat(d1./NormW,1,size(W,2)))).^2,2));
    PBI = zeros(1,size(PopObj,1));
    PBI(Sub) = d1 + 5*d2;
end


text:
% ReadSegy : Reads a SEG Y rev 1 formatted file, and forces Constant Trcae Length
%
% Call :
% [Data,SegyTraceHeaders,SegyHeader]=ReadSegyConstantTraceLength(filename);
%
% See ReadSegy for optional arguments
%
% See also : ReadSegy
%
function [Data,SegyTraceHeaders,SegyHeader]=ReadSegyConstantTraceLength(filename,varargin);

SegyHeader=ReadSegyHeader(filename);
SegyHeader.FixedLengthTraceFlag=1;

[Data,SegyTraceHeaders,SegyHeader]=ReadSegy(filename,varargin);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
%labels segments comparing their eventIdx property to event annotations.
%The eventIdx prperty is set if an eventSegmentation has been used
%to generate the segment
classdef EventSegmentsLabeler < Algorithm

    properties (Access = public)
        manualAnnotations;
    end

    methods (Access = public)

        function obj = EventSegmentsLabeler()
            obj.name = 'eventSegmentsLabeler';
            obj.inputPort = DataType.kSegment;
            obj.outputPort = DataType.kSegment;
        end

        function labeledSegments = compute(obj,segments)
            eventsLabeler = EventsLabeler();
            labels = eventsLabeler.labelEventIdxs([segments.eventIdx],obj.manualAnnotations.eventAnnotations);

            labeledSegments = Helper.LabelSegmentsWithValidLabels(segments,labels);
        end

        function str = toString(obj)
            str = sprintf('%s',obj.name);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function box = boundingBox(points)
%BOUNDINGBOX Bounding box of a set of points
%
%   BOX = boundingBox(POINTS)
%   Returns the bounding box of the set of points POINTS. POINTS can be
%   either a N-by-2 or N-by-3 array. The result BOX is a 1-by-4 or 1-by-6
%   array, containing:
%   [XMIN XMAX YMIN YMAX] (2D point sets)
%   [XMIN XMAX YMIN YMAX ZMIN ZMAX] (3D point sets)
%
%   Example
%   % Draw the bounding box of a set of random points
%     points = rand(30, 2);
%     figure; hold on;
%     drawPoint(points, '.');
%     box = boundingBox(points);
%     drawBox(box, 'r');
%
%   % Draw bounding box of a cubeoctehedron
%     [v e f] = createCubeOctahedron;
%     box3d = boundingBox(v);
%     figure; hold on;
%     drawMesh(v, f);
%     drawBox3d(box3d);
%     set(gcf, 'renderer', 'opengl')
%     axis([-2 2 -2 2 -2 2]);
%     view(3)
%
%   See also
%   polygonBounds, drawBox
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-04-01,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% HISTORY
% 2011-04-08 add example
% 2011-12-09 rename to boundingBox

% compute extreme x and y values
xmin = min(points(:,1));
xmax = max(points(:,1));
ymin = min(points(:,2));
ymax = max(points(:,2));

if size(points, 2) > 2
    % process case of 3D points
    zmin = min(points(:,3));
    zmax = max(points(:,3));

    % format as box 3D data structure
    box = [xmin xmax ymin ymax zmin zmax];
else
    % format as box data structure
    box = [xmin xmax ymin ymax];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [X,Y,indsort] = grid_communities(c)
% GRID_COMMUNITIES       Outline communities along diagonal
%
%   [X Y INDSORT] = GRID_COMMUNITIES(C) takes a vector of community
%   assignments C and returns three output arguments for visualizing the
%   communities. The third is INDSORT, which is an ordering of the vertices
%   so that nodes with the same community assignment are next to one
%   another. The first two arguments are vectors that, when overlaid on the
%   adjacency matrix using the PLOT function, highlight the communities.
%
%   Example:
%
%   >> load AIJ;                                % load adjacency matrix
%   >> [C,Q] = modularity_louvain_und(AIJ);     % get community assignments
%   >> [X,Y,INDSORT] = fcn_grid_communities(C); % call function
%   >> imagesc(AIJ(INDSORT,INDSORT));           % plot ordered adjacency matrix
%   >> hold on;                                 % hold on to overlay community visualization
%   >> plot(X,Y,'r','linewidth',2);             % plot community boundaries
%
%   Inputs:     C,       community assignments
%
%   Outputs:    X,       x coor
%               Y,       y coor
%               INDSORT, indices
%
%   Richard Betzel, Indiana University, 2012
%

%#ok<*AGROW>

nc = max(c);
[c,indsort] = sort(c);

X = [];
Y = [];
for i = 1:nc
    ind = find(c == i);
    if ~isempty(ind)
        mn = min(ind) - 0.5;
        mx = max(ind) + 0.5;
        x = [mn mn mx mx mn NaN];
        y = [mn mx mx mn mn NaN];
        X = [X, x];
        Y = [Y, y];
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ c, sigma ] = covariance_to_correlation ( n, k )

%*****************************************************************************80
%
%% COVARIANCE_TO_CORRELATION: correlation matrix from a covariance matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real K(N,N), the covariance matrix.
%
%    Output, real C(N,N), the correlation matrix.
%
%    Output, real SIGMA(N), the standard deviations.
%
  tol = sqrt ( eps );
%
%  K must be symmetric.
%
  error_frobenius = r8mat_is_symmetric ( n, n, k );

  if ( tol < error_frobenius )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COVARIANCE_TO_CORRELATION - Fatal error!\n' );
    fprintf ( 1, '  Input matrix K fails symmetry test with error %g\n', error_frobenius );
    error ( 'COVARIANCE_TO_CORRELATION - Fatal error!' );
  end
%
%  It must be the case that K(I,J)^2 <= K(I,I) * K(J,J).
%
  error_max = 0.0;
  for i = 1 : n
    for j = i + 1 : n
      error_max = max ( error_max, k(i,j)^2 - k(i,i) * k(j,j) );
    end
  end

  if ( tol < error_max )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COVARIANCE_TO_CORRELATION - Fatal error!\n' );
    fprintf ( 1, '  Input matrix K fails K(I,J)^2 <= K(I,I)*K(J,J)\n' );
    error ( 'COVARIANCE_TO_CORRELATION - Fatal error!' );
  end
%
%  Get the diagonal.
%
  sigma = zeros ( n, 1 );
  for i = 1 : n
    sigma(i) = k(i,i);
  end
%
%  Ensure the diagonal is positive.
%
  sigma_min = min ( sigma(1:n) );

  if ( sigma_min <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COVARIANCE_TO_CORRELATION - Fatal error!\n' );
    fprintf ( 1, '  Input matrix K has nonpositive diagonal entry = %g\n', sigma_min );
    error ( 'COVARIANCE_TO_CORRELATION - Fatal error!' );
  end
%
%  Convert from variance to standard deviation.
%
  sigma(1:n) = sqrt ( sigma(1:n) );
%
%  Form C.
%
  c(1:n,1:n) = diag ( 1.0 ./ sigma(1:n) ) * k(1:n,1:n) * diag ( 1.0 ./ sigma(1:n) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g_hofstadter_test ( )

%*****************************************************************************80
%
%% G_HOFSTADTER_TEST tests G_HOFSTADTER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'G_HOFSTADTER_TEST\n' );
  fprintf ( 1, '  G_HOFSTADTER evaluates Hofstadter''s recursive\n' );
  fprintf ( 1, '  G function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N   G(N)' );
  fprintf ( 1, '\n' );

  for i = 0 : 30
    g = g_hofstadter ( i );
    fprintf ( 1, '  %4d  %4d\n', i, g );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [vertices, sph_verts, faces, fvec, new_name]=align_FOE(filename, confs)

load(filename);
[path,name,ext] = fileparts(filename);

if ~exist('faces', 'var') | ~exist('vertices', 'var') | ~exist('sph_verts', 'var') | ~exist('fvec', 'var')
    disp('One or more of faces, vertices, spherical vertices, or SPHARM descriptor are missing');
    return;
end

switch deblank(char(confs.CPoint))
    case 'x'
        blue = 1;
    case 'y'
        blue = 2;
    case 'z'
        blue = 3;
end

switch deblank(char(confs.NPole))
    case 'x'
        yellow = 1;
    case 'y'
        yellow = 2;
    case 'z'
        yellow = 3;
end

blueyellow = [blue yellow];
degree = confs.MaxSPHARMDegree;

% rotate the parameter space
disp('<< Rotate the parameter space >>');
[fvec, sph_verts, expts] = param_rotate(fvec, vertices, sph_verts, faces, degree, blueyellow);

dirName = [confs.OutDirectory '/alignParam'];
if ~exist(dirName,'dir')
    mkdir(dirName);
end

new_name = sprintf('%s/%sFOE_prm.mat',dirName,name(1:end-3));
if exist(new_name,'file')
    prompt = {'Enter new filename:'};
    dlg_title = 'New File Name';
    num_lines = 1;
    def = {new_name};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    new_name = answer{1};
end
save(new_name, 'vertices', 'sph_verts', 'faces', 'fvec','expts');

% rotate in the object space
svs = [0 0 1; 1 0 0; 0 0 -1]; % north pole, intersection of dateline and equator, south pole
Z = calculate_SPHARM_basis(svs, 1);

disp('<< Rotate the object space >>');
vs = real(Z(:,2:4)*fvec(2:4,:));
R = object_rotate_R(vs);
% R = object_rotate_R([ellipAxes(:,1)';ellipAxes(:,3)']);
fvec = fvec*R';
vertices = vertices*R';

new_name = sprintf('%s/%sFOE_reg.mat',confs.OutDirectory,name(1:end-3));
if exist(new_name,'file')
    prompt = {'Enter new filename:'};
    dlg_title = 'New File Name';
    num_lines = 1;
    def = {new_name};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    new_name = answer{1};
end
save(new_name, 'vertices', 'sph_verts', 'faces', 'fvec','expts');

return;

%
% Factoring a Rotation Matrix as Rz*Ry*Rx (counterclockwise when looking towards the origin)
%

function [thetaX, thetaY, thetaZ] = factor_rot_xyz(R)

thetaY = asin(R(1,3));
if (thetaY < pi/2)
    if (thetaY > -pi/2)
        thetaX = atan2(-R(2,3),R(3,3));
        thetaZ = atan2(-R(1,2),R(1,1));
    else
        disp('WARNING (Factor Rotation): thetaY = -pi/2, not a unique solution, set thetaZ = 0');
        thetaX = -atan2(R(2,1),R(2,2));
        thetaZ = 0;
    end
else
    disp('WARNING (Factor Rotation): thetaY = pi/2, not a unique solution, set thetaZ = 0');
    thetaX = atan2(R(2,1),R(2,2));
    thetaZ = 0;
end

% disp(sprintf('Factor rotation xyz: %0.2f %0.2f %0.2f',thetaX/pi,thetaY/pi,thetaZ/pi));

return;

%
% rotation matrix in object space
%

function R = object_rotate_R(vs)

% fix north pole
[PHI,THETA] = cart2sph(vs(1,1),vs(1,2),vs(1,3));
ind = find(PHI<0); PHI(ind) = PHI(ind)+2*pi;
THETA = pi/2-THETA;
alpha = -PHI; beta = -THETA;
R = rotate_mat(0, beta, 0)*rotate_mat(0, 0, alpha);
vs = vs*R';
% fix intersection;
[PHI,THETA] = cart2sph(vs(2,1),vs(2,2),vs(2,3));
gamma = -PHI;
R1 = rotate_mat(0, 0, gamma); R = R1*R;
vs = vs*R1';

return;

%
% Parameter space rotation using degree 1 ellipsoid
%

function [fvec, sph_verts, expts] = param_rotate(fvec, vertices, sph_verts, faces, degree, blueyellow)

% calculate matrix A
coeffs = fvec(2:4,:);
A(:,1) = (coeffs(1,:)-coeffs(3,:))';
A(:,2) = -(coeffs(1,:)+coeffs(3,:))'*i; % there is a typo in the paper, should be - here.
A(:,3) = sqrt(2)*coeffs(2,:)';
A = real(A*sqrt(3)/(2*sqrt(2*pi)));

% SVD to find rotation and scaling matrics
%   [U,S,V] = svd(X) produces a diagonal matrix S of the same dimension as X, with
%   nonnegative diagonal elements in decreasing order, and unitary matrices U and V so
%   that X = U*S*V'.
% need to rotate object space first.
[U,S,V] = svd(A);
% extremum and saddle points
expts = U*S;

% set up rotation matrix
R = V(:,[3 2 1])';
if (det(R)<0)
    disp('WARNING (Parameter Space): rotoinversion!! Change back to pure rotation');
    R(2,:) = R(2,:)*(-1);
end

[thetaX, thetaY, thetaZ] = factor_rot_xyz(R');

% Question: why multiply thetaY by -1?
R = rotate_mat(thetaX, -thetaY, thetaZ);
disp(sprintf('Rotation xyz: %0.2f %0.2f %0.2f',thetaX/pi,thetaY/pi,thetaZ/pi));

% the parameter space rotation
sph_verts = (R*sph_verts')';
% create new spharm descriptor (degree 1 is enough)
[fvec, d, Z_tmp1, name_temp] = create_SPHARM_des_LSF(vertices, [], sph_verts, 1, '', '');

% calculate the blue point
svs = [0 0 1; 1 0 0]; % yellow (north pole), blue (intersection)
Z = calculate_SPHARM_basis(svs, 1);

vs = real(Z(:,2:4)*fvec(2:4,:));
disp(sprintf('blue: (%f, %f, %f)',vs));
% blue point and yellow point should be on the positive side of
% x (blue=1), y (blue=2), or z (blue=3) axis in the object space
% i.e., vs(blue) and vs(yellow) should be >0
blue = blueyellow(1); yellow = blueyellow(2);
if vs(2,blue)<0 | vs(1,yellow)<0
    if vs(2,blue)<0
        Rfix = rotate_mat(0, 0, pi);
    else
        Rfix = eye(3);
    end
    if vs(1,yellow)<0
        Rfix =  rotate_mat(pi, 0, 0)*Rfix;
    end
	% the parameter space rotation
	sph_verts = (Rfix*sph_verts')';
	% create new spharm descriptor (degree 1 is enough)
    [fvec, d, Z_tmp1, name_temp] = create_SPHARM_des_LSF(vertices, [], sph_verts, 1, '', '');
end

[fvec, d, Z_tmp1, name_temp] = create_SPHARM_des_LSF(vertices, [], sph_verts, degree, '', '');

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus