text:
classdef ComputingInnerAndBoundaryCutMesh < handle

    properties (Access = public)
        innerCutMesh
        cutCoordComputer
    end

    properties (Access = private)
        subCellCases
        cutEdgesComputer
        coord
    end

    properties (Access = private)
        backgroundMesh
        levelSet
        boundaryConnec
    end

    methods (Access = public)

        function obj = ComputingInnerAndBoundaryCutMesh(cParams)
            obj.init(cParams)
        end

        function compute(obj)
            obj.computeInnerCutMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.levelSet = cParams.levelSet;
            obj.backgroundMesh = cParams.backgroundMesh;
            obj.boundaryConnec = cParams.boundaryConnec;
        end

        function computeSubCellCases(obj)
            s.connec     = obj.backgroundMesh.connec;
            s.levelSet   = obj.levelSet;
            obj.subCellCases = SubCellsCasesComputer(s);
            obj.subCellCases.compute();
        end

        function computeCutEdges(obj)
            obj.backgroundMesh.computeEdges();
            e = obj.backgroundMesh.edges;
            s.nodesInEdges = e.nodesInEdges;
            s.levelSet     = obj.levelSet;
            c = CutEdgesComputer(s);
            c.compute();
            obj.cutEdgesComputer = c;
        end

        function computeCutCoordinateComputer(obj)
            e = obj.backgroundMesh.edges;
            sC.coord            = obj.backgroundMesh.coord;
            sC.nodesInEdges     = e.nodesInEdges;
            sC.xCutEdgePoint    = obj.cutEdgesComputer.xCutEdgePoint;
            sC.isEdgeCut        = obj.cutEdgesComputer.isEdgeCut;
            cComputer = CutCoordinatesComputer(sC);
            cComputer.compute();
            obj.cutCoordComputer = cComputer;
        end

        function computeInnerCutMesh(obj)
            obj.computeSubCellCases();
            obj.computeCutEdges();
            obj.computeCutCoordinateComputer();
            obj.coord = obj.cutCoordComputer.coord;

            e = obj.backgroundMesh.edges;
            s.edgesInElem   = e.edgesInElem;
            s.isEdgeCut     = obj.cutEdgesComputer.isEdgeCut;
            isEdgeCut       = EdgeCutInElemComputer(s);
            isEdgeCutInElem = isEdgeCut.compute();

            e = obj.backgroundMesh.edges;
            nEdgesCutCase   = [2 3 4];
            nSubCellsByElem = [3 4 6];

            subCell = cell(length(nEdgesCutCase),1);
            cN = cell(length(nEdgesCutCase),1);

            nCutEdges = sum(isEdgeCutInElem,1);
            for icases = 1:length(nEdgesCutCase)
                t = nCutEdges == nEdgesCutCase(icases);
                isEdgeCutInElemCase = isEdgeCutInElem(:,t);

                s.isEdgeCutInElem = isEdgeCutInElemCase;
                all2Cut = AllEdges2CutEdgesComputer(s);

                cEp.all2Cut = all2Cut;
                cEp.allNodesinElemParams.finalNodeNumber = size(obj.backgroundMesh.coord,1);
                cEp.allNodesinElemParams.connec = obj.backgroundMesh.connec(t,:);
                cEp.allNodesInElemCoordParams.localNodeByEdgeByElem = e.localNodeByEdgeByElem(t,:,:);
                cEp.edgesInElem = e.edgesInElem(t,:);
                cEp.nEdgeByElem = e.nEdgeByElem;
                cEp.isEdgeCut = obj.cutEdgesComputer.isEdgeCut;
                cEp.allNodesInElemCoordParams.xCutEdgePoint = obj.cutEdgesComputer.xCutEdgePoint;
                cEp.isEdgeCutInElem = isEdgeCutInElemCase;
                cE = CutPointsInElemComputer(cEp);
                cE.compute();

                if sum(t) ~= 0

                    cN{icases} = cE;

                    caseInfo = obj.subCellCases.caseInfo{icases};

                    nodes = obj.backgroundMesh.connec;
                    cutCells(:,1) = 1:size(nodes,1);

                    sS.bestSubCellCaseSelector.coord = obj.coord;
                    sA.subMeshConnecParams           = sS;
                    sA.xAllNodesInElem               = cE.xAllNodesInElem;
                    sA.allNodesInElem                = cE.allNodesInElem;
                    sA.subCellCases                  = caseInfo.subCellCases(t,:);

                    sI.allSubCellsConnecParams = sA;
                    sI.isSubCellInterior = caseInfo.isSubCellsInterior(:,t);
                    sI.cutElems = cutCells;

                    sI.nSubCellsByElem = nSubCellsByElem(icases);

                    subCell{icases} = InteriorSubCellsConnecComputer(sI);
                end

            end

            connecT = [];
            xCoordsIso = [];
            cellC = [];

            xCoordsIsoBoundary = [];
            cellContainingSubCellBoundary = [];

            for icase = 1:3
                subC = subCell{icase};
                scN  = cN{icases};
                if ~isempty(subC)
                    connecT = cat(1,connecT,subC.connec);
                    xCoordsIso = cat(3,xCoordsIso,subC.xCoordsIso);
                    cellC = cat(1,cellC,subC.cellContainingSubcell);
                end

                if ~isempty(scN)
                    xCoordsIsoBoundary = cat(3,xCoordsIsoBoundary,scN.xCutInElem);
                    % cellContainingSubCellBoundary = cat(1,cellContainingSubCellBoundary,obj.cutCells);
                end

            end

            connec                = connecT;
            % xCoordsIso            = xCoordsIso;
            cellContainingSubcell = cellC;

            sM.connec = connec;
            sM.coord  = obj.coord;

            m = Mesh(sM);
            s.mesh                  = m;
            s.xCoordsIso            = xCoordsIso;
            s.cellContainingSubcell = cellContainingSubcell;
            inCutMesh = InnerCutMesh(s);

            obj.innerCutMesh = inCutMesh;

        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [A,b,C,f,S,P,RESULTS,YrA] = run_CNMF_patches(data,K,patches,tau,p,options)

% Run the constrained NMF algorithm on a large dataset by operating on
% spatially overlapping patches in parallel and then merging the results.
% The inputs is memory mapped, allowing for large datasets to be processed
% with reduced memory requirements. Processing in patches also allows the
% identification of weaker neurons without the need of normalization.
% The components are also classified by retaining only the components that
% significantly overlap with the active pixels. The active pixels are
% determined by classifying the PSD of the data.

% INPUTS:

% data   :   .mat file containing
%   data.Y      (the data matrix in the original dimensions)
%   data.Yr     (the data matrix reshaped in 2d format)
%   data.sizY   (dimensions of the original dataset)
%   data.nY     (minimum value of dataset)
% OR the original dataset in 3d/4d format in which case a memory mapped
%   file is created

% K      :   number of components to be found in each patch
% patches:   cell array containing the start and end points of each patch
% tau    :   half-size of each cell for initializing the components
% p      :   order of autoregressive progress
% options:   struct for algorithm parameters

% OUTPUTS:

% A      :   Matrix of spatial components
% b      :   Spatial background
% C      :   Matrix of temporal components
% f      :   Temporal background
% P      :   Struct for neuron parameters
% RESULTS:   Results of the CNMF algorithm on individual patches
% YrA    :   Residual signal at the level of each component

% Author: Eftychios A. Pnevmatikakis, Simons Foundation, 2016

memmaped = isobject(data);
if memmaped
    sizY = data.sizY;
else    % create a memory mapped object named data_file.mat
    Y = data;
    clear data;
    sizY = size(Y);
    Yr = reshape(Y,prod(sizY(1:end-1)),[]);
    nY = min(Yr(:));
    %Yr = Yr - nY;
    save('data_file.mat','Yr','Y','nY','sizY','-v7.3');
    data = matfile('data_file.mat','Writable',true);
end

defoptions = CNMFSetParms;
if nargin < 6 || isempty(options)
    options = defoptions;
end

if nargin < 5 || isempty(p)
    p = 0;
end

if nargin < 4 || isempty(tau)
    tau = 5;
end

if nargin < 3 || isempty(patches)
    patches = construct_patches(sizY(1:end-1),[50,50]);
end

if nargin < 2 || isempty(K)
    K = 10;
end

if ~isfield(options,'cl_thr') || isempty(options.cl_thr)
    cl_thr = 0.8;
else
    cl_thr = options.cl_thr;
end

RESULTS(length(patches)) = struct();

%parfor_progress(length(patches)); %monitor parfor progress (requires parfor_progress from mathworks file exchange)
parfor i = 1:length(patches)
    if length(sizY) == 3
        Y = data.Y(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),:);
        [d1,d2,T] = size(Y);
        d3 = 1;
    else
        Y = data.Y(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6),:);
        [d1,d2,d3,T] = size(Y);
    end
    d = d1*d2*d3;
    options_temp = options;
    options_temp.d1 = d1; options_temp.d2 = d2; options_temp.d3 = d3;
    options_temp.nb = 1;
    [P,Y] = preprocess_data(double(Y),p);
    [Ain,Cin,bin,fin,center] = initialize_components(Y,K,tau,options_temp);  % initialize
    Yr = reshape(Y,d,T);
    %clear Y;
    options_temp.use_parallel = 0; % turn off parallel updating for spatial components
    [A,b,Cin] = update_spatial_components(Yr,Cin,fin,Ain,P,options_temp);
    P.p = 0;
    options_temp.temporal_parallel = 0;
    [C,f,P,S] = update_temporal_components(Yr,A,b,Cin,fin,P,options_temp);
    [Am,Cm,K_m,merged_ROIs,P,Sm] = merge_components(Yr,A,b,C,f,P,S,options_temp);
    [A2,b2,Cm] = update_spatial_components(Yr,Cm,f,Am,P,options_temp);
    P.p = p;
    [C2,f2,P2,S2] = update_temporal_components(Yr,A2,b2,Cm,f,P,options_temp);
    RESULTS(i).A = A2;
    RESULTS(i).C = C2;
    RESULTS(i).b = b2;
    RESULTS(i).f = f2;
    RESULTS(i).S = S2;
    RESULTS(i).P = P2;
    fprintf(['Finished processing patch # ',num2str(i),' out of ',num2str(length(patches)), '.\n']);
    %parfor_progress;
end
%parfor_progress(0);

%% combine results into one structure
fprintf('Combining results from different patches...');
d = prod(sizY(1:end-1));
A = sparse(d,length(patches)*K);
P.sn = zeros(sizY(1:end-1));
P.active_pixels = zeros(sizY(1:end-1));
IND = zeros(sizY(1:end-1));
P.b = {};
P.c1 = {};
P.gn = {};
P.neuron_sn = {};

if length(sizY) == 3
    P.psdx = zeros(patches{end}(2),patches{end}(4),size(RESULTS(1).P.psdx,2));
else
    P.psdx = zeros(patches{end}(2),patches{end}(4),patches{end}(6),size(RESULTS(1).P.psdx,2));
end

cnt = 0;
B = sparse(prod(sizY(1:end-1)),length(patches));
MASK = zeros(sizY(1:end-1));
F = zeros(length(patches),sizY(end));
for i = 1:length(patches)
    for k = 1:K
        if k <= size(RESULTS(i).A,2)
            cnt = cnt + 1;
            Atemp = zeros(sizY(1:end-1));
            if length(sizY) == 3
                Atemp(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = reshape(RESULTS(i).A(:,k),patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
            else
                Atemp(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = reshape(RESULTS(i).A(:,k),patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1,patches{i}(6)-patches{i}(5)+1);
            end
            A(:,cnt) = sparse(Atemp(:));
        end
    end
    if length(sizY) == 3
        b_temp = sparse(sizY(1),sizY(2));
        b_temp(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = reshape(RESULTS(i).b,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
        MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) + 1;
        P.sn(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = reshape(RESULTS(i).P.sn,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
        P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) + ...
            reshape(RESULTS(i).P.active_pixels,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1);
        IND(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) = IND(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4)) + 1;
        P.psdx(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),:) = reshape(RESULTS(i).P.psdx,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1,[]);
    else
        b_temp = sparse(sizY(1),sizY(2),sizY(3));
        b_temp(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = reshape(RESULTS(i).b,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1,patches{i}(6)-patches{i}(5)+1);
        MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = MASK(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) + 1;
        P.sn(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = reshape(RESULTS(i).P.sn,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1,patches{i}(6)-patches{i}(5)+1);
        P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = P.active_pixels(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) + ...
            reshape(RESULTS(i).P.active_pixels,patches{i}(2)-patches{i}(1)+1,patches{i}(4)-patches{i}(3)+1,patches{i}(6)-patches{i}(5)+1);
        IND(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) = IND(patches{i}(1):patches{i}(2),patches{i}(3):patches{i}(4),patches{i}(5):patches{i}(6)) + 1;
    end
    P.b = [P.b;RESULTS(i).P.b];
    P.c1 = [P.c1;RESULTS(i).P.c1];
    P.gn = [P.gn;RESULTS(i).P.gn];
    P.neuron_sn = [P.neuron_sn;RESULTS(i).P.neuron_sn];
    B(:,i) = b_temp(:);
    F(i,:) = RESULTS(i).f;
end
A(:,cnt+1:end) = [];
A = spdiags(1./MASK(:),0,prod(sizY(1:end-1)),prod(sizY(1:end-1)))*A;
B = spdiags(1./MASK(:),0,prod(sizY(1:end-1)),prod(sizY(1:end-1)))*B;
C = cell2mat({RESULTS(:).C}');
S = cell2mat({RESULTS(:).S}');
ff = find(sum(A,1)==0);
A(:,ff) = [];
C(ff,:) = [];
S(ff,:) = [];
fprintf(' done. \n');
%% estimate active pixels
fprintf('Classifying pixels...')
if length(sizY) == 3
    X = P.psdx(:,:,1:min(size(P.psdx,2),500));
else
    X = P.psdx(:,:,:,1:min(size(P.psdx,2),500));
end
X = reshape(X,[],size(X,ndims(X)));
X = bsxfun(@minus,X,mean(X,2));     % center
X = spdiags(std(X,[],2)+1e-5,0,size(X,1),size(X,1))\X;
[L,Cx] = kmeans_pp(X',2);
[~,ind] = min(sum(Cx(max(1,end-49):end,:),1));
P.active_pixels = (L==ind);
P.centroids = Cx;
fprintf(' done. \n');
%% merge results
fprintf('Merging overlaping components...')
Am = A;
Cm = C;
Pm = P;
Sm = S;
Km = 0;
Kn = size(A,2);

while Km < Kn
    Kn = size(Am,2);
    [Am,Cm,~,~,Pm,Sm] = merge_components([],Am,[],Cm,[],Pm,Sm,options);
    Km = size(Am,2);
end
fprintf(' done. \n');
%% classify components
%ff = classify_components(Am,Pm,options);
ff = true(size(Am,2),1);
A = Am(:,ff);
C = Cm(ff,:);

%% compute spatial and temporal background using a rank-1 fit

fin = mean(F);
for iter = 1:10
    bin = max(B*(F*fin')/norm(fin)^2,0);
    fin = max((bin'*B)*F/norm(bin)^2,0);
end
%% update spatial components
fprintf('Updating spatial components...');
options.d1 = sizY(1);
options.d2 = sizY(2);
if length(sizY) == 4; options.d3 = sizY(3); end
[A,b,C] = update_spatial_components(data,C,fin,A,Pm,options);
fprintf(' done. \n');
%% update temporal components
fprintf('Updating temporal components... ')
Pm.p = p;
options.temporal_iter = 2;
[C,f,P,S,YrA] = update_temporal_components(data,A,b,C,fin,Pm,options);
fprintf(' done. \n');

author: zhoupc


text:
function [ point, seed ] = p15_sample ( m, n, seed )

%*****************************************************************************80
%
%% P15_SAMPLE samples points from the region in problem 15.
%
%  Discussion:
%
%    With a little bit of work, we can guarantee that we don't have to
%    use a rejection method.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real POINT(M,N), the coordinates of the points.
%
%    Output, integer SEED, a seed for the random number generator.
%
  x1 = -8.0;
  x2 =  2.0;
  x3 = -2.0;
  x4 =  8.0;

  y1 = -1.0;
  y2 =  0.0;
  y3 =  0.0;
  y4 =  1.0;
%
%  Generate a batch of points in [0,1]x[0,1].
%
  [ point, seed ] = r8mat_uniform_01 ( m, n, seed );
%
%  Generate a batch of N probabilities.
%
  [ prob, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Map half the points into [X1,X2] x [Y1,Y2].
%
  i1 = find ( prob(1:n) < 0.5 );

  point(1,i1) = x1 + point(1,i1) * ( x2 - x1 );
  point(2,i1) = y1 + point(2,i1) * ( y2 - y1 );
%
%  Map the other points into [X3,X4] x [Y3,Y4].
%
  i2 = find ( 0.5 <= prob(1:n) );

  point(1,i2) = x3 + point(1,i2) * ( x4 - x3 );
  point(2,i2) = y3 + point(2,i2) * ( y4 - y3 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function[varargout]=vtranspose(varargin)
%VTRANSPOSE  Transpose multiple input arguments simultaneously.
%
%   [Y1,Y2, ... YN]=VTRANSPOSE(X1,X2, ... XN) is equivalent to
%
%       Y1=CONJ(X1'); Y2=CONJ(X2'); ... YN=CONJ(XN');
%
%   for matrix arguments X1,X2, ... XN.  This is the matrix transpose.
%   Recall X1' is the conjugate transpose of X1 in Matlab's notation.
%
%   VTRANSPOSE(X1,X2,...XN);  with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014 J.M. Lilly --- type 'help jlab_license' for details

for i=1:length(varargin)
  %varargout{i}=permute(varargin{i},[2 1]);
  varargout{i}=conj(varargin{i}');
end

eval(to_overwrite(nargin))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%% OPT_BATCH_LOCAL uses the BATCH command to run the OPT code locally.
%
%  Discussion:
%
%    The OPT code is a function, so first we must write a script
%    called OPT_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'OPT_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run OPT_SCRIPT locally.\n' );
%
%  BATCH sends the script for execution.
%
  my_job = batch ( 'opt_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'opt_fun' }, ...
    'matlabpool', 4 );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  LOAD makes the script's workspace available:
%
%  z_star = optimal control.
%  PAR = problem parameters.
%  exit = positive if calculation was successful.
%
  load ( my_job );
%
%  On successful conclusion, display the plot.
%
  if ( 0 <= exit && isfield ( PAR, 'plot' ) )
    feval ( PAR.plot, z_star, PAR )
  end

%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'OPT_BATCH_LOCAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function [label, center, bCon, sumD, D] = litekmeans(X, k, varargin)

if nargin < 2
    error('litekmeans:TooFewInputs','At least two input arguments required.');
end

[n, p] = size(X);

pnames = {   'distance' 'start'   'maxiter'  'replicates' 'onlinephase' 'clustermaxiter'};
dflts =  {'sqeuclidean' 'sample'       []        []        'off'              []        };
[eid,errmsg,distance,start,maxit,reps,online,clustermaxit] = getargs(pnames, dflts, varargin{:});
if ~isempty(eid)
    error(sprintf('litekmeans:%s',eid),errmsg);
end

if ischar(distance)
    distNames = {'sqeuclidean','cosine'};
    j = strcmpi(distance, distNames);
    j = find(j);
    if length(j) > 1
        error('litekmeans:AmbiguousDistance', ...
            'Ambiguous ''Distance'' parameter value:  %s.', distance);
    elseif isempty(j)
        error('litekmeans:UnknownDistance', ...
            'Unknown ''Distance'' parameter value:  %s.', distance);
    end
    distance = distNames{j};
else
    error('litekmeans:InvalidDistance', ...
        'The ''Distance'' parameter value must be a string.');
end

center = [];
if ischar(start)
    startNames = {'sample','cluster'};
    j = find(strncmpi(start,startNames,length(start)));
    if length(j) > 1
        error(message('litekmeans:AmbiguousStart', start));
    elseif isempty(j)
        error(message('litekmeans:UnknownStart', start));
    elseif isempty(k)
        error('litekmeans:MissingK', ...
            'You must specify the number of clusters, K.');
    end
    if j == 2
        if floor(.1*n) < 5*k
            j = 1;
        end
    end
    start = startNames{j};
elseif isnumeric(start)
    if size(start,2) == p
        center = start;
    elseif (size(start,2) == 1 || size(start,1) == 1)
        center = X(start,:);
    else
        error('litekmeans:MisshapedStart', ...
            'The ''Start'' matrix must have the same number of columns as X.');
    end
    if isempty(k)
        k = size(center,1);
    elseif (k ~= size(center,1))
        error('litekmeans:MisshapedStart', ...
            'The ''Start'' matrix must have K rows.');
    end
    start = 'numeric';
else
    error('litekmeans:InvalidStart', ...
        'The ''Start'' parameter value must be a string or a numeric matrix or array.');
end

% The maximum iteration number is default 100
if isempty(maxit)
    maxit = 100;
end

% The maximum iteration number for preliminary clustering phase on random
% 10% subsamples is default 10
if isempty(clustermaxit)
    clustermaxit = 10;
end

% Assume one replicate
if isempty(reps) || ~isempty(center)
    reps = 1;
end

if ~(isscalar(k) && isnumeric(k) && isreal(k) && k > 0 && (round(k)==k))
    error('litekmeans:InvalidK', ...
        'X must be a positive integer value.');
elseif n < k
    error('litekmeans:TooManyClusters', ...
        'X must have more rows than the number of clusters.');
end

bestlabel = [];
sumD = zeros(1,k);
bCon = false;

for t=1:reps
    switch start
        case 'sample'
            center = X(randsample(n,k),:);
        case 'cluster'
            Xsubset = X(randsample(n,floor(.1*n)),:);
            [dump, center] = litekmeans(Xsubset, k, varargin{:}, 'start','sample', 'replicates',1 ,'MaxIter',clustermaxit);
        case 'numeric'
    end

    last = 0;label=1;
    it=0;

    switch distance
        case 'sqeuclidean'
            while any(label ~= last) && it<maxit
                last = label;

                bb = full(sum(center.*center,2)');
                ab = full(X*center');
                D = bb(ones(1,n),:) - 2*ab;

                [val,label] = min(D,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    %disp([num2str(k-length(ll)),' clusters dropped at iter ',num2str(it)]);
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);

                    aa = sum(X.*X,2);
                    val = aa + val;
                    [dump,idx] = sort(val,1,'descend');
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)',0,k,k))'*X);    % compute center of each cluster
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if it>=maxit
                        aa = full(sum(X.*X,2));
                        bb = full(sum(center.*center,2));
                        ab = full(X*center');
                        D = bsxfun(@plus,aa,bb') - 2*ab;
                        D(D<0) = 0;
                    else
                        aa = full(sum(X.*X,2));
                        D = aa(:,ones(1,k)) + D;
                        D(D<0) = 0;
                    end
                    D = sqrt(D);
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if it>=maxit
                    aa = full(sum(X.*X,2));
                    bb = full(sum(center.*center,2));
                    ab = full(X*center');
                    D = bsxfun(@plus,aa,bb') - 2*ab;
                    D(D<0) = 0;
                else
                    aa = full(sum(X.*X,2));
                    D = aa(:,ones(1,k)) + D;
                    D(D<0) = 0;
                end
                D = sqrt(D);
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
        case 'cosine'
            while any(label ~= last) && it<maxit
                last = label;
                W=full(X*center');
                [val,label] = max(W,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);
                    [dump,idx] = sort(val);
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)',0,k,k))'*X);    % compute center of each cluster
                centernorm = sqrt(sum(center.^2, 2));
                center = center ./ centernorm(:,ones(1,p));
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if any(label ~= last)
                        W=full(X*center');
                    end
                    D = 1-W;
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if any(label ~= last)
                    W=full(X*center');
                end
                D = 1-W;
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
    end
end

label = bestlabel;
center = bestcenter;
if reps>1
    sumD = bestsumD;
    D = bestD;
elseif nargout > 3
    switch distance
        case 'sqeuclidean'
            if it>=maxit
                aa = full(sum(X.*X,2));
                bb = full(sum(center.*center,2));
                ab = full(X*center');
                D = bsxfun(@plus,aa,bb') - 2*ab;
                D(D<0) = 0;
            else
                aa = full(sum(X.*X,2));
                D = aa(:,ones(1,k)) + D;
                D(D<0) = 0;
            end
            D = sqrt(D);
        case 'cosine'
            if it>=maxit
                W=full(X*center');
            end
            D = 1-W;
    end
    for j = 1:k
        sumD(j) = sum(D(label==j,j));
    end
end

function [eid,emsg,varargout]=getargs(pnames,dflts,varargin)
%GETARGS Process parameter name/value pairs
%   [EID,EMSG,A,B,...]=GETARGS(PNAMES,DFLTS,'NAME1',VAL1,'NAME2',VAL2,...)
%   accepts a cell array PNAMES of valid parameter names, a cell array
%   DFLTS of default values for the parameters named in PNAMES, and
%   additional parameter name/value pairs.  Returns parameter values A,B,...
%   in the same order as the names in PNAMES.  Outputs corresponding to
%   entries in PNAMES that are not specified in the name/value pairs are
%   set to the corresponding value from DFLTS.  If nargout is equal to
%   length(PNAMES)+1, then unrecognized name/value pairs are an error.  If
%   nargout is equal to length(PNAMES)+2, then all unrecognized name/value
%   pairs are returned in a single cell array following any other outputs.
%
%   EID and EMSG are empty if the arguments are valid.  If an error occurs,
%   EMSG is the text of an error message and EID is the final component
%   of an error message id.  GETARGS does not actually throw any errors,
%   but rather returns EID and EMSG so that the caller may throw the error.
%   Outputs will be partially processed after an error occurs.
%
%   This utility can be used for processing name/value pair arguments.
%
%   Example:
%       pnames = {'color' 'linestyle', 'linewidth'}
%       dflts  = {    'r'         '_'          '1'}
%       varargin = {{'linew' 2 'nonesuch' [1 2 3] 'linestyle' ':'}
%       [eid,emsg,c,ls,lw] = statgetargs(pnames,dflts,varargin{:})    % error
%       [eid,emsg,c,ls,lw,ur] = statgetargs(pnames,dflts,varargin{:}) % ok

% We always create (nparams+2) outputs:
%    one each for emsg and eid
%    nparams varargs for values corresponding to names in pnames
% If they ask for one more (nargout == nparams+3), it's for unrecognized
% names/values

%   Original Copyright 1993-2008 The MathWorks, Inc.
%   Modified by Deng Cai (dengcai@gmail.com) 2011.11.27

% Initialize some variables
emsg = '';
eid = '';
nparams = length(pnames);
varargout = dflts;
unrecog = {};
nargs = length(varargin);

% Must have name/value pairs
if mod(nargs,2)~=0
    eid = 'WrongNumberArgs';
    emsg = 'Wrong number of arguments.';
else
    % Process name/value pairs
    for j=1:2:nargs
        pname = varargin{j};
        if ~ischar(pname)
            eid = 'BadParamName';
            emsg = 'Parameter name must be text.';
            break;
        end
        i = strcmpi(pname,pnames);
        i = find(i);
        if isempty(i)
            % if they've asked to get back unrecognized names/values, add this
            % one to the list
            if nargout > nparams+2
                unrecog((end+1):(end+2)) = {varargin{j} varargin{j+1}};
                % otherwise, it's an error
            else
                eid = 'BadParamName';
                emsg = sprintf('Invalid parameter name:  %s.',pname);
                break;
            end
        elseif length(i)>1
            eid = 'BadParamName';
            emsg = sprintf('Ambiguous parameter name:  %s.',pname);
            break;
        else
            varargout{i} = varargin{j+1};
        end
    end
end

varargout{nparams+1} = unrecog;

author: jwyang


text:
function [fh, plotImage] = plot(this, varargin)
% plots an MR image
%
%   Y  = MrImage
%   fh = Y.plot('ParameterName', ParameterValue)
%
% Note:
%   The first (specified) dimension will be plotted from left to right
%   The 2nd (specified) dimensions will be plotted from down to up
%   The 3rd (and higher) specified dimensions will be plotted in tiles or
%   figures (depending on options)
%
% IN
%   varargin    'ParameterName', 'ParameterValue'-pairs for the following
%               properties:
%               'plotType'          Type of plot that is created
%                                       'montage'   images are plotted as
%                                                   montages
%                                       'labeledMontage'
%                                                   as montage, but with
%                                                   labels (default)
%                                       'spm'       uses display functions
%                                                   from SPM (spm_display/
%                                                   spm_check_registration)
%                                                   to visualize 3D volumes
%                                                   with header information
%                                                   applied ("world space")
%                                                   Note: if multiple
%                                                   selected volumes are
%                                                   specified,
%                                                   spm_check_registration
%                                                   is used
%                                       'spmInteractive' /'spmi'
%                                                   same as SPM, but keeps
%                                                   temporary nifti files to
%                                                   allow clicking into spm
%                                                   figure
%                                       '3D'/'3d'/'ortho'
%                                                   See also tapas_uniqc_view3d
%                                                   Plots 3 orthogonal
%                                                   sections
%                                                   (with CrossHair) of
%                                                   3D image interactively
%               'linkOptions'       link another real-time plot to input
%                                   (e.g. mouse) on this one, using
%                                   positions on current plot as update to
%                                   plotDim-values for linked plot
%
%                                   shortcut string:
%                                   'timeseries_<plotDim>'/'ts_<plotDim>'
%                                       default: 'ts_4'/'ts_t'
%                                           plots time series of voxel that
%                                           mouse is currently pointing to
%                                       OR
%
%                                   struct of class MrLinkPlotOptions with fields/values
%                                   'plotType'  - 'timeseries' ...
%                                   'plotDim'   - dimension which should be
%                                                 plotted, default: 4
%                                   'fixedDimsPoint'
%                                               - cell(1,2*nFixedDims)
%                                                 dimension label/index array pairs that
%                                                 shall remain fixed/selected
%                                                 and are not on original plot
%                                                 default: {5 ,1, 6, 1,
%                                                 ..., nDims, 1}
%
%               'displayRange'      [1,2] vector for pixel value = black and
%                                                    pixel value = white
%               'overlay'           false (default) or true if an overlay
%                                   image is plotted
%               'signalPart'        for complex data, defines which signal
%                                   part shall be extracted for plotting
%                                       'all'       - take signal as is
%                                                     (default for
%                                                     real-valued data)
%                                       'abs'       - absolute value
%                                                     (default for complex
%                                                     data)
%                                       'phase'     - phase of signal
%                                       'real'      - real part of signal
%                                       'imag'      - imaginary part of
%                                                     signal
%               'plotMode'          transformation of data before plotting
%                                   'linear' (default), 'log'
%               'rotate90'          default: 0; 0,1,2,3; rotates image
%                                   by multiple of 90 degrees AFTER
%                                   flipping slice dimensions
%               'useSlider'         true or false
%                                   provides interactive slider for
%                                   slices/volumes
%               'colorMap'          string, any matlab colormap name
%                                   e.g. 'jet', 'gray'
%               'colorBar',         'on' or 'off' (default)
%                                   where applicable, determines whether
%                                   colorbar with displayRange shall be plotted
%                                   in figure;
%               'overlayImages'     (cell of) MrImages that will be
%                                   overlayed
%               'overlayMode'       'edge', 'mask', 'map'
%                                   'edge'  only edges of overlay are
%                                           displayed
%                                   'mask'  every non-zero voxel is
%                                           displayed (different colors for
%                                           different integer values, i.e.
%                                           clusters'
%                                   'map'   thresholded map in one colormap
%                                           is displayed (e.g. spmF/T-maps)
%                                           thresholds from
%                                           overlayThreshold
%               'overlayThreshold'  [minimumThreshold, maximumThreshold]
%                                   tresholds for overlayMode 'map'
%                                   default: [-Inf, Inf] = [minValue, maxValue]
%                                   everything below minValue will not be
%                                   displayed;
%                                   everything above maxValue
%                                   will have brightest color
%               'overlayAlpha'      transparency value of overlays
%                                   (0 = transparent; 1 = opaque; default: 0.1)
%               'edgeThreshold'     determines where edges will be drawn,
%                                   the higher, the less edges
%                                   Note: logarithmic scale, e.g. try 0.005
%                                   if 0.05 has too little edges
%
%               data selection      data selection uses MrImage.select/
%                                   MrImage.dimInfo.select
%               'dimLabel'          scalar or vector with array indices or
%                                   sampling points
%               'imagePlotDim'      1x3 vector of dimensions or cell array
%                                   of labels that constitue the image that
%                                   will be plotted, i.e. which three
%                                   dimensions define the volume that will
%                                   be plotted within one figure, e.g.
%                                   [1,2,3] (default)
%                                   {'x', 'y', 'z'}
%               'selectionType'     'index' (default) or 'label' selects
%                                   hoew 'dimlabel' vector is interpreted,
%                                   as array indices or sampling points
%
%               Orientation changes:
%               'sliceDimension'    (default: 3) determines which dimension
%                                   shall be plotted as a slice
%                                   can be entered as index or dimLabel
%               'rotate90'          default: 0; 0,1,2,3; rotates image
%                                   by multiple of 90 degrees AFTER
%                                   flipping slice dimensions
%               for montage plots:
%               'nRows'             default: NaN (automatic calculation)
%               'nCols'             default NaN (automatic calculation)
%               'FontSize'          font size of tile labels of montage
%               'plotTitle'         if true, title i.e. readible name of
%                                   image is put in plot
%               'plotLabels'        if true, slice labels are put into
%                                   montage
%
% OUT
%   fh          [nFigures,1] vector of figure handles
%
% EXAMPLE
%
%   Y.plot('z', [6:10])
%   Y.plot('displayRange', [0 1000])
%   Y.plot('useSlider', true, 'z', []);
%
%   See also MrLinkPlotOptions

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-05-21
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% check whether image object has data
if isempty(this.data)
    error('tapas:uniqc:MrImageEmptyDataMatrix', ...
        'Data matrix empty for MrImage-object %s', this.name);
end
%% set defaults
% default signal part and plot mode
if isreal(this)
    defaults.signalPart         = 'all';
else
    defaults.signalPart         = 'abs';
end
defaults.plotMode               = 'linear';

% data selection
defaults.selectionType          = 'index';

% plot appearance
defaults.plotType               = 'labeledMontage';

defaults.FigureSize             = [1600 900];
defaults.nRows                  = NaN; % automatically determined by image size
defaults.nCols                  = NaN;
defaults.FontSize               = 10;
defaults.plotTitle              = true;
defaults.plotLabels             = true;

defaults.rotate90               = 0;
defaults.sliceDimension         = 3;
defaults.displayRange           = [];
defaults.useSlider              = false;
defaults.colorMap               = 'gray';
defaults.colorBar               = 'off';
defaults.windowStyle            = 'docked';
defaults.imagePlotDim           = [1,2,3];

% overlay parameters
defaults.overlay                = false;
defaults.overlayImages          = {};
defaults.overlayMode            = 'mask';
defaults.overlayThreshold       = [];
defaults.overlayAlpha           = []; % depends on overlayMode
defaults.edgeThreshold          = [];

% linked plot options
defaults.linkOptions             = [];

% get arguments
[args, ~] = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

% check colorbar and overlays
doPlotColorBar = strcmpi(colorBar, 'on');
% plot overlays if overlay images are given, but not when plotType is spm -
% they are added via spm_check_registration
% make sure overlay images are cell
if ~iscell(overlayImages)
    overlayImages = {overlayImages};
end
doPlotOverlays = (overlay || ~isempty(overlayImages)) ...
    && ~(strcmpi(plotType, 'spm') || strcmpi(plotType, 'spmi'));

% convert imagePlotDim from label to index
if iscell(imagePlotDim)
    [~, imagePlotDim] = ismember(imagePlotDim, this.dimInfo.dimLabels);
end

doLinkPlot = ~isempty(linkOptions);
doMontage = ismember(lower(plotType), {'montage', 'labeledmontage'});

if doLinkPlot
    if ~isa(linkOptions, 'MrLinkPlotOptions')
        if ischar(linkOptions) % shortcut string to create options
            if contains(linkOptions, {'ts', 'timeseries'})
                % of the form timeseries_<iDim>, retrieve the second)
                splitString = regexp(linkOptions, '_', 'split');
                iDimLinkedPlot = str2num(splitString{2});
                linkOptions = MrLinkPlotOptions('ts', this.dimInfo, ...
                    imagePlotDim, iDimLinkedPlot);
            else
                error('tapas:uniqc:MrImage:InvalidLinkOptions', ...
                    'linkOptions must be a MrLinkOptions object or a shortcut string');
            end
        end
    end
end

%% select plot data as plotImage (dimension selection)

% check whether any input parameters specify which data to plot
plotDataSpecified = ismember(varargin(1:2:end), this.dimInfo.dimLabels);
% copy plot image for selection
plotImage = this.copyobj;

isPlotDataSpecified = any(plotDataSpecified);

% select plot data
if isPlotDataSpecified
    plotDataSpecified = repmat(plotDataSpecified, 2, 1);
    plotDataSpecified = reshape(plotDataSpecified, 1, []);
    stringSelection = varargin(plotDataSpecified);
else
    stringSelection = {};
    if ~useSlider % default: no slider used
        % 1 image with all samples of first three dimensions, for all further
        % dimensions only first sample is plotted
        if plotImage.dimInfo.nDims > 3
            nDimsSelect = plotImage.dimInfo.nDims - 3;
            dimLabelsSelect = plotImage.dimInfo.dimLabels;
            stringSelection(1:2:nDimsSelect*2) = dimLabelsSelect(4:end);
            stringSelection(2:2:nDimsSelect*2) = {1};
            [plotImage, selectionIndexArray] = plotImage.select('type', selectionType, ...
                stringSelection{:});
        else
            % create selectionIndexArray for labels later
            % simply select everything (by selecting everything in the
            % first dimension
            [~, selectionIndexArray] = plotImage.dimInfo.select(...
                plotImage.dimInfo.dimLabels{1}, 1:plotImage.dimInfo.nSamples(1));
        end
    else % use slider
        % 1 image with all samples of first FOUR dimensions, for all further
        % dimensions only first sample is plotted
        if plotImage.dimInfo.nDims > 4
            nDimsSelect = plotImage.dimInfo.nDims - 4;
            dimLabelsSelect = plotImage.dimInfo.dimLabels;
            stringSelection(1:2:nDimsSelect*2) = dimLabelsSelect(5:end);
            stringSelection(2:2:nDimsSelect*2) = {1};
        end

    end
end
[plotImage, selectionIndexArray] = plotImage.select('type', selectionType, ...
    stringSelection{:});

%% extract signal part, plot mode and display range

% signal part
switch signalPart
    case 'all'
        % do nothing, leave dataPlot as is
    case 'abs'
        plotImage = abs(plotImage);
    case {'angle', 'phase'}
        plotImage = angle(plotImage) + pi;
    case 'real'
        plotImage = real(plotImage);
    case 'imag'
        plotImage = imag(plotImage);
end

% linear or logarithmic plot
switch plotMode
    case 'linear' % nothing happens
    case 'log'
        plotImage = log(abs(plotImage));
end

% display range
if isempty(displayRange)
    if islogical(plotImage.data) % for logical arrays (masks)
        displayRange = [0 1];
    else
        displayRange = [min(plotImage), ...
            prctile(plotImage,98)];
        % check whether valid display range has been specified
        % if values equal or second range larger then first, specify by
        % hand
        if diff(displayRange) <= 0
            displayRange(2) = displayRange(1) + 1;
        end
    end
end

% Manipulate orientation for plot
if ischar(sliceDimension) % convert dimLabel to index
    sliceDimension = plotImage.dimInfo.get_dim_index(sliceDimension);
end
nDims = plotImage.dimInfo.nDims;
switch sliceDimension
    case 1
        permuteArray = [3 2 1 4 5:nDims];
        plotImage = permute(plotImage, permuteArray(1:nDims));
        selectionIndexArray = selectionIndexArray(permuteArray(1:nDims));
    case 2
        permuteArray = [1 3 2 4 5:nDims];
        plotImage = permute(plotImage, permuteArray(1:nDims));
        selectionIndexArray = selectionIndexArray(permuteArray(1:nDims));
    case 3
        %   as is...
    otherwise
        plotImage = permute(plotImage, [1 2 sliceDimension]);
        selectionIndexArray = selectionIndexArray([3 2 sliceDimension]);
end

if rotate90
    plotImage = rot90(plotImage, rotate90);
end

%% extract data for overlay image
if doPlotOverlays

    % check background image is 3D image
    nDimsPlotImage = sum(plotImage.dimInfo.nSamples > 1);
    is3dBackground = nDimsPlotImage < 4;
    if ~is3dBackground
        error('tapas:uniqc:MrImage:plot:UnderlayNot3D', ...
            ['Background image is not 3D but has ', ...
            num2str(nDimsPlotImage), ' dimensions.']);
    end
    % extract data from background image
    % extract plot data and sort
    plotData = squeeze(plotImage.data);
    backgroundNSamples = size(plotData);
    % set default Alpha depending on define mode'
    if isempty(overlayAlpha)
        switch overlayMode
            case {'mask', 'edge'}
                overlayAlpha = 1;
            case 'blend'
                overlayAlpha = 0.2;
            case 'map'
                overlayAlpha = 0.7;
            otherwise
                overlayAlpha = 0.1;
        end
    end

    % settings
    nColorsPerMap   = 256;
    overlayImages   = reshape(overlayImages, [], 1);

    % loop over overlays and extract data
    nOverlays       = numel(overlayImages);
    dataOverlays    = cell(nOverlays,1);

    for iOverlay = 1:nOverlays
        thisOverlay = overlayImages{iOverlay};

        %% for map: overlayThreshold image only,
        %  for mask: binarize
        %  for edge: binarize, then compute edge

        switch overlayMode
            case {'map', 'maps'}
                thisOverlay.threshold(overlayThreshold);
            case {'mask', 'masks'}
                thisOverlay.threshold(0, 'exclude');
            case {'edge', 'edges'}
                thisOverlay.threshold(0, 'exclude');
                % for cluster mask with values 1, 2, ...nClusters,
                % leave values of edge same as cluster values
                thisOverlay = edge(thisOverlay,'sobel', edgeThreshold);
        end

        if any(plotDataSpecified)
            stringSelection = varargin(plotDataSpecified);
            [plotOverlay, ~, ~] = thisOverlay.select('type', selectionType, ...
                stringSelection{:});
        else
            plotOverlay = thisOverlay.copyobj;
        end

        switch sliceDimension
            case 1
                permuteArray = [3 2 1 4];
                plotOverlay = permute(plotOverlay, permuteArray(1:nDims));
            case 2
                permuteArray = [1 3 2 4];
                plotOverlay = permute(plotOverlay, permuteArray(1:nDims));
            case 3
                %   as is...
            otherwise
                plotOverlay = permute(plotOverlay, [1 2 sliceDimension]);
        end

        % apply rotation
        if rotate90
            plotOverlay = rot90(plotOverlay, rotate90);
        end
        % extract plot data and sort
        dataOverlays{iOverlay} = squeeze(plotOverlay.data);

        % check that background and overlay image have same dimension
        overlayNSamples = size(dataOverlays{iOverlay});
        equalDimBackgroundOverlay = ...
            numel(overlayNSamples) == numel(backgroundNSamples) && ...
            all(overlayNSamples == backgroundNSamples);
        if any(~equalDimBackgroundOverlay)
            error('tapas:uniqc:MrImage:plot:NumberOfSamplesUnderlayOverlayMismatch', ...
                ['Different number of samples for background (', ...
                num2str(backgroundNSamples), ') and overlay image (', ...
                num2str(size(dataOverlays{iOverlay})), ').']);
        end
    end

    % Define color maps for different cases:
    %   map: hot
    %   mask/edge: one color per mask image, faded colors for different
    %   clusters within same mask

    functionHandleColorMaps = {
        @hot
        @cool
        @spring
        @summer
        @winter
        @jet
        @hsv
        };

    overlayColorMap = cell(nOverlays,1);
    switch overlayMode
        case {'mask', 'edge', 'masks', 'edges'}
            baseColors = hsv(nOverlays);

            % determine unique color values and make color map
            % a shaded version of the base color
            for iOverlay = 1:nOverlays
                indColorsOverlay = unique(dataOverlays{iOverlay});
                nColorsOverlay = max(2, round(...
                    max(indColorsOverlay) - min(indColorsOverlay)));
                overlayColorMap{iOverlay} = tapas_uniqc_get_brightened_color(...
                    baseColors(iOverlay,:), 1:nColorsOverlay - 1, ...
                    nColorsOverlay -1, 0.7);

                % add for transparency
                overlayColorMap{iOverlay} = [0,0,0; ...
                    overlayColorMap{iOverlay}];
            end

        case {'map', 'maps'}
            for iOverlay = 1:nOverlays
                overlayColorMap{iOverlay} = ...
                    functionHandleColorMaps{iOverlay}(nColorsPerMap);
            end

    end

    % Assemble RGB-image for montage by adding overlays with transparency as
    % RGB in right colormap
    rangeOverlays   = cell(nOverlays, 1);
    rangeImage      = cell(nOverlays, 1);

    for iOverlay = 1:nOverlays
        [plotData, rangeOverlays{iOverlay}, rangeImage{iOverlay}] = ...
            tapas_uniqc_add_overlay(plotData, dataOverlays{iOverlay}, ...
            overlayColorMap{iOverlay}, ...
            overlayThreshold, ...
            overlayAlpha, ...
            displayRange);
    end
end
%% plot

% slider view
if useSlider
    % extract plot data
    plotData = squeeze(plotImage.data);
    % make sure only data with 4 or less dimension is used
    nDimsPlotImage = sum(plotImage.dimInfo.nSamples > 1);
    is4dor3dPlotImage = (nDimsPlotImage == 3 || nDimsPlotImage == 4);
    if ~is4dor3dPlotImage
        error('tapas:uniqc:MrImage:plot:SliderImageTooManyDimensions', ...
            ['Selected plot image is not 3D or 4D but has ', ...
            num2str(nDimsPlotImage), ' dimensions.']);
    end
    nSlices = plotImage.dimInfo.nSamples(3);
    tapas_uniqc_slider4d(plotData, @(Y, iDynSli, fh, yMin, yMax) ...
        tapas_uniqc_plot_abs_image(Y, iDynSli, fh, yMin, yMax, colorMap, colorBar), ...
        nSlices, displayRange(1), displayRange(2), this.name);

else % different plot types: montage, 3D, spm
    switch lower(plotType)
        case {'montage', 'labeledmontage'} % this is the default setting
            % make labels
            if strcmpi(plotType, 'labeledMontage') && plotImage.dimInfo.nDims >= 3
                stringLabels = cellfun(@(x,y) sprintf('%3.1f [%d]',x,y), ...
                    num2cell(plotImage.dimInfo.samplingPoints{imagePlotDim(3)}),...
                    num2cell(selectionIndexArray{imagePlotDim(3)}),...
                    'UniformOutput', false);
            else
                stringLabels = [];
            end

            % which dims need their own figure, i.e. are not in the image?
            dimsWithFig = setdiff(1:plotImage.dimInfo.nDims, imagePlotDim);
            if isempty(dimsWithFig), dimsWithFig = 4; end % for 3D data
            % how many additional dims are given
            nDimsWithFig = length(dimsWithFig);
            % extract plot data and sort
            if ~doPlotOverlays
                plotData = permute(plotImage.data, [imagePlotDim, dimsWithFig]);
                % number of samples in imagePlotDim
                nSamplesImagePlotDim = plotImage.dimInfo.nSamples(imagePlotDim(1:min(plotImage.dimInfo.nDims,3)));
                % reshape plot data to 4D matrix
                if plotImage.dimInfo.nDims > 3
                    % number of samples in dimsWithFig
                    nSamplesDimsWithFig = plotImage.dimInfo.nSamples(dimsWithFig);
                    plotData = reshape(plotData, ...
                        nSamplesImagePlotDim(1), nSamplesImagePlotDim(2), nSamplesImagePlotDim(3), []);
                else
                    % number of samples in dimsWithFig
                    nSamplesDimsWithFig = 1;
                end
                % total number of figures
                nFigures = size(plotData, 4);
            else
                nDimsWithFig = 1;
                nFigures = 1;
                nSamplesDimsWithFig = 1;
            end

            % now plot
            for n = 1:nFigures
                % make title string
                titleString = [];
                % sampling positions for titleString
                samplingPosArray = cell(1,nDimsWithFig);
                % convert index to subscript values for titleString
                [samplingPosArray{:}] = ind2sub(nSamplesDimsWithFig, n);
                % loop over nDimsWithFig
                for nTitle = 1:nDimsWithFig % number of labels in the title
                    % pos of label in dimInfo.dimLabel
                    labelPos = dimsWithFig(nTitle);
                    % build title string from label and corresponding sampling position

                    if labelPos <= plotImage.dimInfo.nDims % 3D and smaller, have no label!
                        titleString = [titleString, ...
                            plotImage.dimInfo.dimLabels{labelPos}, ...
                            num2str(plotImage.dimInfo.samplingPoints{labelPos}(samplingPosArray{nTitle}), ...
                            '%4.0f') ' ']; %#ok<AGROW>
                    end
                end

                % add info to figure title, if only one slice
                if numel(stringLabels) == 1
                    titleString = plotImage.dimInfo.index2label(1,3);
                    titleString = titleString{1}{1};
                end

                titleString = tapas_uniqc_str2label([plotImage.name, ' ', titleString]);
                % open figure
                fh(n,1) = figure('Name', titleString, 'Position', ...
                    [1 1 FigureSize(1), FigureSize(2)], 'WindowStyle', windowStyle);
                % montage
                if doPlotOverlays
                    thisPlotData = plotData;
                else
                    thisPlotData = permute(plotData(:,:,:,n), [1, 2, 4, 3]);
                end
                if plotLabels
                    [~, montageSize] = tapas_uniqc_labeled_montage(thisPlotData, ...
                        'DisplayRange', displayRange, ...
                        'LabelsIndices', stringLabels, ...
                        'Size', [nRows nCols], ...
                        'FontSize', FontSize);
                else
                    montage(thisPlotData, ...
                        'DisplayRange', displayRange, ...
                        'ThumbnailSize', [], ...
                        'Size', [nRows nCols]);
                end

                resolutions = abs(plotImage.dimInfo.resolutions);
                samplingWidths = abs(plotImage.dimInfo.samplingWidths);
                % use sampling width if resolution if lost
                posNanResolutions = find(isnan(resolutions));
                if ~isempty(posNanResolutions)
                    for idxRes = posNanResolutions
                        % check sampling widths
                        if isnan(samplingWidths)
                            resolutions(idxRes) = 1;
                        else
                            resolutions(idxRes) = samplingWidths(idxRes);
                        end
                    end
                end
                resolutions((end+1):3) = 1;
                resolutions(4:end) = [];
                set(gca, 'DataAspectRatio', resolutions);

                % Display title, colorbar, colormap, if specified
                if plotTitle
                    title(titleString);
                end

                if doPlotColorBar
                    colorbar;
                end
                colormap(gca, colorMap);
                drawnow;
            end
            %

        case {'3d', 'ortho'}

            % check plot image is 3D image
            nDimsPlotImage = sum(plotImage.dimInfo.nSamples > 1);
            is3dPlotImage = nDimsPlotImage == 3;
            if ~is3dPlotImage
                error('tapas:uniqc:MrImage:plot:3DPlotImageNot3D', ...
                    ['Selected plot image is not 3D but has ', ...
                    num2str(nDimsPlotImage), ' dimensions.']);
            end
            % get voxel size ratio
            nonSingleDims = plotImage.dimInfo.nSamples ~=1;
            voxelSizeRatio = abs(plotImage.dimInfo.resolutions);
            voxelSizeRatio = abs(voxelSizeRatio(nonSingleDims));
            % call tapas_uniqc_view3d on plotImage data
            tapas_uniqc_view3d(squeeze(plotImage.data), voxelSizeRatio);
            if doPlotOverlays
                disp('Overlay function for plotType 3d not yet implemented.');
            end
        case {'spm', 'spminteractive', 'spmi'}
            % calls spm_image-function (for single volume) or
            % spm_check_registration (multiple volumes)

            fileNameNifti = plotImage.write_temporary_nifti_for_spm();
            if iscell(fileNameNifti) && numel(fileNameNifti) == 1
                fileNameNifti = fileNameNifti{1};
            elseif ischar(fileNameNifti)
            else
                error('tapas:uniqc:MrImage:plot:SPMPlotTooManyDimensions', ...
                    'SPM plots not implemented for 5+dimensional data yet');
            end

            % select Volumes
            fileNameVolArray = strvcat(tapas_uniqc_get_vol_filenames(fileNameNifti));

            % check if additional (overlay) images have been specified
            doPlotAdditionalImages = ~isempty(overlayImages);
            if doPlotAdditionalImages
                nAddImages = numel(overlayImages);
                for iAddImages = 1:nAddImages
                    fileNameAdditionalNiftis = overlayImages{iAddImages}.write_temporary_nifti_for_spm();
                    if iscell(fileNameAdditionalNiftis) && numel(fileNameAdditionalNiftis) == 1
                        fileNameAdditionalNiftis = fileNameAdditionalNiftis{1};
                    elseif ischar(fileNameAdditionalNiftis) && ~isempty(fileNameAdditionalNiftis)
                    else
                        error('tapas:uniqc:MrImage:plot:SPMPlotTooManyDimensions', ...
                            'High dimensional plotting with SPM not implemented yet');
                    end
                    volArrayFileNameNiftiAddImages{iAddImages} = strvcat(tapas_uniqc_get_vol_filenames(fileNameAdditionalNiftis));
                end

                fileNameVolArray = strvcat(fileNameVolArray, ...
                    volArrayFileNameNiftiAddImages{:});
            else
                fileNameAdditionalNiftis = '';
            end

            % display image
            nImages = size(fileNameVolArray, 1);
            if nImages == 1
                % use display option if only one image selected
                spm_image('Display', fileNameVolArray);

            else
                nMaxImages = 24;
                % check reg all if less than 25 (SPM only supports up to 24
                % volumes)
                % check if filenames are unique first
                spm_check_registration(fileNameVolArray(1:min(nImages,nMaxImages), :));
                if nImages > nMaxImages
                    warning('Only first 24 volumes are displayed');
                end
            end

            switch lower(plotType)
                case {'spminteractive', 'spmi'}
                    input('Press Enter to leave interactive mode');
            end

            % clean up temporary nifti files
            tapas_uniqc_delete_with_hdr(fileNameNifti);
            [~,~] = rmdir(fileparts(fileNameNifti));
            tapas_uniqc_delete_with_hdr(fileNameAdditionalNiftis);
            [~,~] = rmdir(fileparts(fileNameAdditionalNiftis));

    end % plotType
end % use Slider

if doLinkPlot
    ha = gca;
    hi = findobj(ha.Children,'Type','Image');
    hf = gcf;
    stringTitle = sprintf('Linked timeseries Plot %s', this.name);
    hFigLinePlot = figure('Name', stringTitle);
    hAxLinePlot = axes;

    iZ = find(cellfun(@(x) strcmpi(x, 'z'), stringSelection));
    if ~isempty(iZ)
        idxSlicePlotted = stringSelection{iZ+1};
    else
        idxSlicePlotted = 1;
    end

    if doMontage
        dimInfoSelection = plotImage.dimInfo;
        % conversion of coordinates follows from image size and number of
        % slices put into montage rows/columns
        linkOptions.convertMousePosToSelection = ...
            @(x) tapas_uniqc_convert_montage_position_to_selection(x, montageSize, ...
            dimInfoSelection, selectionIndexArray);
    else
        % single slice plot
        linkOptions.convertMousePosToSelection = @(x) [x(2) x(1) idxSlicePlotted];
    end

    hCallback = @(x,y) tapas_uniqc_lineplot_callback(x, y, this, hAxLinePlot, ...
        linkOptions.convertMousePosToSelection);
    ha.ButtonDownFcn = hCallback;
    hi.ButtonDownFcn = hCallback;
    hf.WindowButtonMotionFcn  = hCallback;

end
end

author: translationalneuromodeling

repo: tapas


text:
function [nP,dP]=paderm(tau,r,m)
c(1)=1;
for i=2:r+m+1, c(i)=-c(i-1)*tau/(i-1); end
w=-c(r+2:m+r+1)';
vv=[c(r+1:-1:1)'; zeros(m-1-r,1)];
W=rot90(hankel(c(m+r:-1:r+1),vv));
V=rot90(hankel(c(r:-1:1)));
x=[1 (W\w)'];
dP=x(m+1:-1:1)/x(m+1);
y=[1 x(2:r+1)*V'+c(2:r+1)];
nP=y(r+1:-1:1)/x(m+1);

author: Sable

repo: mcbench-benchmarks


text:
function mappedX = gda(X, Y, no_dims, varargin)
%GDA Perform Generalized Discriminant Analysis
%
%	mappedX = gda(X, Y, no_dims)
%	mappedX = gda(X, Y, no_dims, kernel)
%	mappedX = gda(X, Y, no_dims, kernel, param1)
%	mappedX = gda(X, Y, no_dims, kernel, param1, param2)
%
% Perform Generalized Discriminant Analysis. GDA or Kernel LDA is the
% nonlinear generalization of LDA by means of the kernel trick. X is the
% data on which to perform GDA, Y are the corresponding labels.
% The value of kernel determines the used kernel. Possible values are 'linear',
% 'gauss', 'poly', 'subsets', or 'princ_angles' (default = 'gauss'). For
% more info on setting the parameters of the kernel function, type HELP
% GRAM.
% The function returns the locations of the embedded trainingdata in
% mappedX.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Process inputs
    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    kernel = 'gauss';
    param1 = 1;
	param2 = 0;
    if length(varargin) > 0 & strcmp(class(varargin{1}), 'char'), kernel = varargin{1}; end
	if length(varargin) > 1 & strcmp(class(varargin{2}), 'double'), param1 = varargin{2}; end
	if length(varargin) > 2 & strcmp(class(varargin{3}), 'double'), param2 = varargin{3}; end

	% Make sure labels are nice
	[foo, bar, Y] = unique(Y, 'rows');

	% Get dimensions
	[n, dim] = size(X);
	nclass = max(Y);

	% Sort data according to labels
	[foo, ind] = sort(Y);
	Y = Y(ind);
	X = X(ind,:);

	% Compute kernel matrix
    disp('Computing kernel matrix...');
    K = gram(X, X, kernel, param1, param2);

    % Compute centering matrix
    ell = size(X, 1);
    D = sum(K) / ell;
    E = sum(D) / ell;
    J = ones(ell, 1) * D;
    K = K - J - J' + E * ones(ell, ell);

    % Perform eigenvector decomposition of kernel matrix (Kc = P * gamma * P')
    disp('Performing eigendecomposition of kernel matrix...');
    K(isnan(K)) = 0;
    K(isinf(K)) = 0;
    [P, gamma] = eig(K);

	if size(P, 2) < n
		error('Singularities in kernel matrix prevent solution.');
    end

	% Sort eigenvalues and vectors in descending order
	[gamma, ind] = sort(diag(gamma), 'descend');
	P = P(:,ind);

	% Remove eigenvectors with relatively small value
	minEigv = max(gamma) / 1e5;
	ind = find(gamma > minEigv);
	P = P(:,ind);
	gamma = gamma(ind);
	rankK = length(ind);

	% Recompute kernel matrix
    K = P * diag(gamma) * P';

	% Construct diagonal block matrix W
	W = [];
	for i=1:nclass
		num_data_class = length(find(Y == i));
		W = blkdiag(W, ones(num_data_class) / num_data_class);
	end

	% Determine target dimensionality of data
	old_no_dims = no_dims;
    no_dims = min([no_dims rankK nclass]);
    if old_no_dims > no_dims
        warning(['Target dimensionality reduced to ' num2str(no_dims) '.']);
    end

	% Perform eigendecomposition of matrix (P' * W * P)
	disp('Performing GDA eigenanalysis...');
    [Beta, lambda] = eig(P' * W * P);
	lambda = diag(lambda);

	% Sort eigenvalues and eigenvectors in descending order
	[lambda, ind] = sort(lambda, 'descend');
	Beta = Beta(:,ind(1:no_dims));

	% Compute final embedding mappedX
	mappedX = P * diag(1 ./ gamma) * Beta;

	% Normalize embedding
	for i=1:no_dims
		mappedX(:,i) = mappedX(:,i) / sqrt(mappedX(:,i)' * K * mappedX(:,i));
	end

repo: cluster


text:
function [wmap, wconv] = mri_field_map_pcg_ls(yik, etime, varargin)
% This function is a wrapper for fmap_est_pcg_ls.m to make it callable in a
% similar fashion as mri_field_map_reg.m

%|	yik	[(N) nset ncoil]	complex images at different echo times
%|				(these should be undistorted)
%|				(N) can be (nx,ny) or (nx,ny,nz) etc.
%|	etime	[nset]		vector of echo times [usually in seconds]
%|
%| options
%|	see fmap_est_pcg_ls.m for options
%|
%| out
%|	wmap	[(N)]	regularized estimate of field map [rad/sec]
%|	wconv	""	conventional estimate from first two scans
%|
%| Written 2021-06-03, Melissa Haskell, University of Michigan

if nargin >= 1 && streq(yik, 'test')
	wmap = mri_field_map_pcg_ls_test(yik, varargin{:});
	if ~nargout, clear wmap, end
return
end

%% call function
func_name = 'fmap_est_pcg_ls';
[wmap, wconv] = call_fmap(yik, etime, func_name, varargin{:});

end

% mri_field_map_pcg_ls_test()
% built-in test/example, copied from mri_field_map_reg_test
function wmap = mri_field_map_pcg_ls_test(~, varargin)

% simulate multiple images based on a given mag & field map.
printm 'simulate noisy data'
etime = [0 2 10] * 1e-3; % echo times
nset = length(etime);
R2 = 20; % 1/sec decay
SNR = 20; % dB

wtrue = 2*pi * ir_get_data(fullfile('mri','2001-phase-data','fieldmap128.fld')); % "true" fieldmap
mag = ir_get_data(fullfile('mri','2001-phase-data','mag128.fld')); % "true" magnitude
[nx, ny] = size(mag); % 128

im plc 2 4
im(1, wtrue/(2*pi), 'true field map', [-40, 128]), cbar('Hz')
im(5, mag, 'true mag'), cbar

image_power = 10*log10(sum(sum(mag.^2))/(nx*ny));
noise_power = image_power - SNR;
noise_std = sqrt(10^(noise_power/10));
noise_std = noise_std / 2; % because complex

yik = zeros(nx,ny,nset);
for kk=1:nset
	yik(:,:,kk) = mag ...
		.* exp(1i * wtrue * (etime(kk) - etime(1))) ...
		.* exp(-R2 * (etime(kk) - etime(1)));
end
rng(0)
yik = yik + noise_std * (randn(size(yik)) + 1i * randn(size(yik)));

%yik(1:10,1:10,:) = 0; % stress test all 0 regions

printm 'estimate field map'
[wmap, wconv] = mri_field_map_pcg_ls(yik, etime, 'l2b', -6, varargin{:});

mask = mag > 0.05 * max(mag(:));
im(8, mask), cbar
titlef('Mask for RMSE')

im(2, wconv / (2*pi), 'Conventional field map', [-40,128]), cbar('Hz')
err = (wconv - wmap) / (2*pi);
im(6, err, 'error', [-25 25]), cbar('Hz')
xlabelf('RMSE = %.1f Hz', sqrt(mean(err(mask).^2)))

im(3, wmap / (2*pi), 'Regularized field map', [-40,128]), cbar('Hz')
err = (wtrue - wmap) / (2*pi);
im(7, err, 'error', [-25 25]), cbar('Hz')
xlabelf('RMSE = %.1f Hz', sqrt(mean(err(mask).^2)))
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [kJ] = erg2kJ(erg)
% Convert energy or work from ergs to kilojoules.
% Chad A. Greene 2012
kJ = erg*1e-10;


text:
% changeaxisunits(factorx,factory)
% changeaxisunits(factorx,factory,hax)
% Changes the units of the x and y axes.
% Suppose the original units of the x axis is frames and the original unit
% of the y-axis is pixels, and we want to change to seconds and
% centimeters, respectively. If a frame is .05 seconds and a pixel is .025
% cm, then we would call
% changeaxisunits(.05,.025)
% currently, changes all lines, patches, and text.
%
function changeaxisunits(factorx,factory,hax)

if nargin < 3,
  hax = gca;
end

chil = findobj(hax,'-property','xdata','-and','-property','ydata','-and','-not','type','hggroup');
%oktypes = {'line','patch','text','hggroup'};
for i = 1:length(chil),
  if chil(i) == hax,
    continue;
  end
  x = get(chil(i),'xdata');
  y = get(chil(i),'ydata');
  set(chil(i),'xdata',x*factorx);
  set(chil(i),'ydata',y*factory);
end
chil = findobj(hax,'-property','position','-and','units','data');
%oktypes = {'line','patch','text','hggroup'};
for i = 1:length(chil),
  if chil(i) == hax,
    continue;
  end
  pos = get(chil(i),'position');
  pos(1) = pos(1)*factorx;
  pos(2) = pos(2)*factory;
  set(chil(i),'position',pos);
end
ax = axis;
ax(1:2) = ax(1:2)*factorx;
ax(3:4) = ax(3:4)*factory;
axis(ax);

if strcmpi(get(hax,'xtickmode'),'manual'),
  xtick = get(hax,'xtick');
  xtick = xtick*factorx;
  set(hax,'xtick',xtick);
end

if strcmpi(get(hax,'ytickmode'),'manual'),
  ytick = get(hax,'ytick');
  ytick = ytick*factory;
  set(hax,'ytick',ytick);
end

author: kristinbranson

repo: JAABA


text:
classdef Folder < Datasets.CDDataset
    properties (Constant)
        dirT1 = 't1';
        dirT2 = 't2';
        dirGT = 'gt';
        fmt = {'.png', '.bmp', '.tif'}  % Reserved
    end
    methods
        function obj = Folder(data_path)
            obj@Datasets.CDDataset(data_path)
        end
        function initFileSys(obj)
            % Note that in this initial version, Folder will process
            % ALL non-dir files with any extension name in the subdirs
            %
            % Note that the correspondance is deduced from the alphabetical
            % order of the files inside each subdir, which means that you
            % have to ensure when files in a folder are sorted
            % alphabetically, then the first readable file in t1 subdir
            % corresponds to the first one in t2 subdir, and also the first
            % one in gt subdir, etc.
            listingT1 = dir(fullfile(obj.dataPath, obj.dirT1));
            listingT2 = dir(fullfile(obj.dataPath, obj.dirT2));
            listingGT = dir(fullfile(obj.dataPath, obj.dirGT));
            obj.t1List = obj.getFileList__(listingT1);
            obj.t2List = obj.getFileList__(listingT2);
            obj.refList = obj.getFileList__(listingGT);
        end
    end
    methods(Access=private)
        function paths = getFileList__(obj, listing)
            % Sort in ascending order to ensure the alignment of files from
            % different subfolders
            [names, idxs] = sort({listing.name});
            folders = {listing.folder};
            folders = folders(idxs);
            flags = {listing.isdir};
            flags = flags(idxs);
            paths = {};
            for ii = 1:size(names, 2)
                if ~flags{ii}
                    paths = [paths, fullfile(folders{ii}, names{ii})];
                end
            end
        end
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function test_glomin_all ( )

%*****************************************************************************80
%
%% TEST_GLOMIN_ALL tests the Brent global minimizer on all test functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_GLOMIN_ALL\n' );
  fprintf ( 1, '  Test the Brent GLOMIN routine, which seeks\n' );
  fprintf ( 1, '  a global minimizer of a function F(X)\n' );
  fprintf ( 1, '  in an interval [A,B],\n' );
  fprintf ( 1, '  given some upper bound M for F".\n' );

  machep = eps;
  e = sqrt ( machep );
  t = sqrt ( machep );

  a = 7.0;
  b = 9.0;
  c = ( a + b ) / 2.0;
  m = 0.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_01, ...
    'h_01(x) = 2 - x' );

  a = 7.0;
  b = 9.0;
  c = ( a + b ) / 2.0;
  m = 100.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_01, ...
    'h_01(x) = 2 - x' );

  a = -1.0;
  b = +2.0;
  c = ( a + b ) / 2.0;
  m = 2.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_02, ...
    'h_02(x) = x * x' );

  a = -1.0;
  b = +2.0;
  c = ( a + b ) / 2.0;
  m = 2.1;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_02, ...
    'h_02(x) = x * x' );

  a = -0.5;
  b =  +2.0;
  c = ( a + b ) / 2.0;
  m = 14.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_03, ...
    'h_03(x) = x^3 + x^2' );

  a = -0.5;
  b =  +2.0;
  c = ( a + b ) / 2.0;
  m = 28.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_03, ...
    'h_03(x) = x^3 + x^2' );

  a = -10.0;
  b = +10.0;
  c = ( a + b ) / 2.0;
  m = 72.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_04, ...
    'h_04(x) = ( x + sin(x) ) * exp(-x*x)' );

  a = -10.0;
  b = +10.0;
  c = ( a + b ) / 2.0;
  m = 72.0;

  test_glomin_one ( a, b, c, m, machep, e, t, @h_05, ...
    'h_05(x) = ( x - sin(x) ) * exp(-x*x)' );

  return
end


text:
function [SPE_limit,T2_limit] = comtupeLimit(model)
% DESCRIPTION
% Compute the control limit of T2 statistic and SPE statistic
%
%       [SPE_limit,T2_limit] = comtupeLimit(model)
%
% INPUT
%   model       KPCA model
%
% OUTPUT
%   SPE_limit   the control limit of SPE statistic
%   T2_limit    the control limit of T2 statistic
%
% Created on 18th April 2019, by Kepeng Qiu.
%-------------------------------------------------------------%

% Compute Hotelling's T2 statistic
% T2 = diag(model.mappedX/diag(model.lambda(1:model.dims))*model.mappedX');

% Compute the squared prediction error (SPE)
temp1 = model.K_c*model.V_s;
temp2 = model.K_c*model.V_s(:,1:model.dims);
SPE = sum(temp1.^2,2)-sum(temp2.^2 ,2);

% Compute the T2 limit (the Chi-square Distribution)
k = model.dims*(model.L-1)/(model.L-model.dims);
T2_limit = k *finv(model.beta,model.dims,model.L-model.dims);

% Compute the SPE limit (the Chi-square Distribution)
a = mean(SPE);
b = var(SPE);
g = b/2/a;
h = 2*a^2/b;
SPE_limit = g*chi2inv(model.beta,h);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
clear;
clc
close all;

%% 已知
w = [2 -3 5]';
dt = 0.01; %姿态更新周期: 0.01s = 100Hz
dur = 100; %积分时长 单位s

N = dur / dt; %积分次数
eul = zeros(N, 3);
fprintf("已知b系角速度(陀螺仪输出为):%.3fdeg/s %.3fdeg/s %.3fdeg/s\n", w(1), w(2), w(3));

Cb2n = eye(3);
theta = deg2rad(w)*dt; %定轴转动下: 等效旋转矢量 = 角增量 = 角速度*dt

for i = 1:N
    C_m2m_1 = ch_rv2m(theta);
    Cb2n = Cb2n * C_m2m_1;

    %单位阵正交化
    Cb2n = ch_mnormlz(Cb2n);

    %记录每一步的欧拉角
    eul(i,:) = rad2deg(ch_m2eul(Cb2n));
end

plot(eul)
legend("PITCH(deg)", "ROLL(deg)", "YAW(deg)");

tmp =eul(end,:);
fprintf("最终欧拉角: pich:%.4f° roll:%.4f° yaw:%.4f°\n", tmp(1), tmp(2), tmp(3));

author: yandld


text:
function [y] = spm_MNpdf (m, C, x)
% Evaluate a Multivariate Gaussian PDF
% FORMAT [y] = spm_MNpdf (m, C, x)
%
% m     [d x 1] mean
% C     [d x d] covar
% x     [n x d] points at which to evaluate
%
% y     [n x 1] density at n points
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_MNpdf.m 3995 2010-07-13 17:19:49Z karl $

ic     = inv(C);
[n, d] = size(x);

x    = x - ones(n, 1)*m(:)';
fact = sum(((x*ic).*x), 2);

y    = exp(-0.5*fact);
y    = y./sqrt((2*pi)^d*det(C));

author: spm


text:
function [lp,lj,tr] = findline(pa,pb,pp,varargin)
%FINDLINE "point-on-line" queries in d-dimensional space.
%   [LP,LI] = FINDLINE(PA,PB,PI) finds the set of d-dimensi-
%   onal line-segments that intersect with a given spatial
%   query. Lines are specified as a set of endpoints [PA,PB]
%   where both PA and PB are NL-by-ND arrays of coordinates,
%   where ND is the number of dimensions.
%
%   A set of intersecting lines is returned for each query
%   point in PI, such that the II-th point is associated
%   with the lines LI(LP(II,1):LP(II,2)). Unenclosed points
%   have LP(II,1) == 0.
%
%   [LP,LI,TR] = FINDLINE(PA,PB,PI) additionally returns the
%   supporting aabb-tree used internally to compute the que-
%   ry. If the underlying collection [PA,PB] is static, the
%   tree TR may be passed to subsequent calls, via
%   [LP,LI,TR] = FINDLINE(PA,PB,PI,TR). This syntax may lead
%   to improved performance, especially when the number of
%   lines is large w.r.t. the number of query points. Note
%   that in such cases the distribution of underlying lines
%   is NOT permitted to change between calls, or erroneous
%   results may be returned. Additional parameters used to
%   govern the creation of the underlying aabb-tree may be
%   passed via [...] = FINDLINE(...,TR,OP). See MAKETREE for
%   additional information.
%
%   See also MAKETREE, FINDTRIA, FINDBALL, LINELINE

% Please see the following for additional information:
%
%   Darren Engwirda, "Locally-optimal Delaunay-refinement &
%   optimisation-based mesh generation". Ph.D. Thesis, Scho-
%   ol of Mathematics and Statistics, Univ. of Sydney, 2014:
%   http://hdl.handle.net/2123/13148

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 02/08/2017
%-----------------------------------------------------------

    lp = []; lj = []; tr = []; op = [];

%---------------------------------------------- basic checks
    if (nargin < +3 || nargin > +5)
        error('findline:incorrectNumInputs', ...
            'Incorrect number of inputs.');
    end

    if (nargin >= +4), tr = varargin{1}; end
    if (nargin >= +5), op = varargin{2}; end

%------------------------------ quick return on empty inputs
    if (isempty(pa)), return; end
    if (isempty(pb)), return; end

%---------------------------------------------- basic checks
    if (~isnumeric(pa) || ~isnumeric(pb) || ...
        ~isnumeric(pp) )
        error('findline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end

    if (ndims(pa) ~= +2 || size(pa,2) < +2 || ...
        ndims(pb) ~= +2 || size(pb,2) < +2 )
        error('findline:incorrectDimensions', ...
            'Incorrect input dimensions.');
    end
    if (ndims(pp) ~= +2 || ...
            size(pa,1) ~= size(pb,1) || ...
            size(pa,2) ~= size(pp,2) || ...
            size(pb,2) ~= size(pp,2) )
        error('findline:incorrectDimensions', ...
            'Incorrect input dimensions.');
    end

    if (~isempty(tr) && ~isstruct(tr) )
        error('findline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end
    if (~isempty(op) && ~isstruct(op) )
        error('findline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end

    if (isempty(tr))
%------------------------------ compute aabb-tree for d-line
        nd = size(pp,2) ;
        nl = size(pa,1) ;

        ab = zeros(nl,nd*2);  % compute aabb-tree
        for ax = +1:nd
        ab(:,ax+nd*0) = min(pa(:,ax), ...
                            pb(:,ax)) ;
        ab(:,ax+nd*1) = max(pa(:,ax), ...
                            pb(:,ax)) ;
        end
        tr = maketree(ab,op) ;
    end

%------------------------------ compute tree-to-vert mapping
    tm = mapvert (tr,pp);

%------------------------------ compute intersection rel-tol
    p0 = min([pa; pb],[],+1) ;
    p1 = max([pa; pb],[],+1) ;

    zt = max(p1-p0) * eps^.8 ;

%------------------------------ compute vert-to-line queries
   [li,ip,lj] = queryset( ...
       tr,tm,@linekern,pp,pa,pb,zt) ;

%------------------------------ re-index onto full obj. list
    lp = zeros(size(pp,1),2) ;
    lp( :,2) = -1 ;

    if (isempty(li)), return ; end

    lp(li,:) = ip ;

end

function [ip,il] = linekern(pk,lk,pp,pa,pb,zt)
%LINEKERN d-dim. node//line intersection kernel routine.

        mp = length(pk);
        ml = length(lk);

    %-------------------------- push line/vert onto n*m tile
        pk = pk.' ;

        pk = pk(ones(ml,1),:);
        pk = pk(:);
        lk = lk(:,ones(1,mp));
        lk = lk(:);

    %-------------------------- compute O(n*m) intersections
        mm = (pa(lk,:)+pb(lk,:)) * +.5 ;
        DD = (pb(lk,:)-pa(lk,:)) * +.5 ;

        mp = mm-pp(pk,:) ;

        tt =-sum(mp.*DD,2)./ ...
             sum(DD.*DD,2) ;
        tt = max(min(tt,+1.),-1.);

        nd = size(pp,2);

        qq = mm + repmat(tt,1,nd) .* DD;

        on = ...
        sum((pp(pk,:)-qq).^2,2) <= zt^2;

        ip = pk(on) ;
        il = lk(on) ;

end

author: CHLNDDEV


text:
function subs = tt_ind2sub(siz,idx)
%TT_IND2SUB Multiple subscripts from linear indices.
%
%   SUBS = TT_IND2SUB(SIZ,INDS) returns that subscripts equivalent
%   to the linear indices in INDS for a tensor of size SIZ.
%
%   See also TT_SUB2IND, IND2SUB.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if isempty(idx)
    subs = [];
    return;
end

k = [1 cumprod(siz(1:end-1))];
n = length(siz);
idx = idx - 1;
for i = n : -1 : 1
    subs(:,i) = floor(idx / k(i)) + 1;
    idx = rem(idx,k(i));
end


text:
function [M_R] = vox2ras_rsolveAA(Vc_C, inPlaneRotation, varargin)
%%
%% NAME
%%
%%     vox2ras_rsolveAA.m (vox2ras_r{otation}solveA{uto}A{lign})
%%
%% AUTHOR
%%
%%	Rudolph Pienaar
%%
%% VERSION
%%
%% 	$Id: vox2ras_rsolveAA.m,v 1.5 2011/03/02 00:04:13 nicks Exp $
%%
%% SYNOPSIS
%%
%%     [M_R] = vox2ras_rsolveAA(Vc_C, inPlaneRotation, ch_orientation)
%%
%% ARGUMENTS
%%
%%      Vc_C		in      column vector defining a direction cosine for
%%					a volume
%%	inPlaneRotation	in	scalar float defining the in plane rotation as
%%					read from the data meas.asc file
%%	ch_orientation	in	optional character string definining the plane of
%%					the main slab orientation
%%	M_R		out	vox2ras rotational candidate
%%
%% DESCRIPTION
%%
%%	"vox2ras_rsolveAA" attempts to find a candidate vox2ras rotation matrix using
%%	techniques described in autoaligncorrect.cpp.
%%
%%	The optional character string, ch_orientation, defaults to 's' (sagittal)
%%	and defines the orientation of the main slab. This string is one of:
%%
%%		's'	sagittal
%%		'c'	coronal
%%		't'	transverse
%%
%% PRECONDITIONS
%%
%%	o The vector C is read from a Siemens meas.asc file such that
%%		ci	= sSliceArray.asSlice[0].sNormal.dSag
%%		cj	= sSliceArray.asSlice[0].sNormal.dCor
%%		ck	= sSliceArray.asSlice[0].sNormal.dTra
%%
%% POSTCONDITIONS
%%
%%	o All returned matrices are 4x4.
%%	o Only the rotations of the vox2ras matrix are determined by this function.
%%		The center of k-space is not determined.
%%
%% SEE ALSO
%%
%%	vox2ras_rsolve	- determine the rotational component of a vox2ras matrix
%%				using Siemens reference orientations indirectly
%%	vox2ras_ksolve	- determine the k-space col in RAS of a vox2ras matrix
%%	vox2ras_dfmeas	- main function: determines the vox2ras matrix from a
%%			  Siemens meas.asc file.
%%
%% HISTORY
%%
%% 02 June 2004
%% o Initial design and coding.
%%

%
% vox2ras_rsolveAA.m
%
% Original Author: Rudolph Pienaar
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.5 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

M_R		= zeros(4, 4);
Vc_Cn		= Vc_C./(norm(Vc_C));
%Vc_Cn		= Vc_C;
ch_orientation	= 's';
if length(varargin)
	ch_orientation	= varargin{1}(1);
end

%% phase reference vector -
%%	adapted from Andre van der Kouwe's "autoaligncorrect.cpp"
Vc_P	= zeros(3, 1);
switch ch_orientation
    case 't'
    	Vc_P(1)	= 0;
	Vc_P(2)	=  Vc_Cn(3)*sqrt(1/(Vc_Cn(2)*Vc_Cn(2)+Vc_Cn(3)*Vc_Cn(3)));
	Vc_P(3)	= -Vc_Cn(2)*sqrt(1/(Vc_Cn(2)*Vc_Cn(2)+Vc_Cn(3)*Vc_Cn(3)));
    case 'c'
	Vc_P(1)	=  Vc_Cn(2)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
	Vc_P(2)	= -Vc_Cn(1)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
    	Vc_P(3)	= 0;
    case 's'
	Vc_P(1)	= -Vc_Cn(2)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
	Vc_P(2)	=  Vc_Cn(1)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
    	Vc_P(3)	= 0;
    otherwise
        fprintf(1, 'Unknown orientation parameter passed. Returning with dummy M_R');
	return;
end

%% The readout reference vector is the cross product of Vc_Cn and Vc_P
Vc_R	= cross(Vc_Cn, Vc_P);

M_R(1:3, 1)	= Vc_P;
M_R(1:3, 2)	= Vc_R;
M_R(1:3, 3)	= Vc_C;

%% The above calculated rotation matrices define an (x,y) plane
%%	given by the first two column vectors. These reference
%%	vectors need to be rotated by inPlaneRotation to
%%	reach the final rotational vox2ras.

theta_f		= inPlaneRotation;

M3_Mu	= [	 cos(theta_f)	 sin(theta_f)	0
		-sin(theta_f)	 cos(theta_f)	0
		 	0		0	1];
M3_R	= M_R(1:3,1:3)	* M3_Mu;
M_R(1:3,1:3)	= M3_R;

%% The MGH vox2ras matrix inverts the Readout column
M_R		= M_R * [ 1  0  0  0
			  0 -1  0  0
			  0  0  1  0
			  0  0  0  1];

%% All done!

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% DEMO_smoothCurve.m
%
% This script demonstrates how to use smoothCurve to create a piecewise
% smooth curve that is twice differentiable
%

% The number of knots is the number of segments + 1
nKnot = 5;

% P = function [value; slope; curvature] at knot points
P = [...
    2*(1-2*rand(1,nKnot));  %Fucntion value
    1*randn(1,nKnot);   %Function slope
    5*randn(1,nKnot)];  %Function curvature

% T is the location of each knot point:
T = cumsum(0.8+0.2*rand(1,nKnot));

% Evaluate the piecewise function:
t = linspace(T(1),T(end),250);
[x,dx,ddx] = pwPoly5(T,P,t);

% Plot the results:
figure(87); clf;
subplot(3,1,1); hold on;
plot(t,x,'k-');
plot(T,P(1,:),'bo','MarkerSize',10)
ylabel('x')
title('Test Function')
subplot(3,1,2); hold on;
plot(t,dx,'k-');
plot(T,P(2,:),'bo','MarkerSize',10)
ylabel('dx')
subplot(3,1,3); hold on;
plot(t,ddx,'k-');
plot(T,P(3,:),'bo','MarkerSize',10)
ylabel('ddx')
xlabel('t')

author: MatthewPeterKelly

repo: dscTutorials


text:
function pdf = english_sentence_length_pdf ( x )

%*****************************************************************************80
%
%% ENGLISH_SENTENCE_LENGTH_PDF evaluates the English Sentence Length PDF.
%
%  Discussion:
%
%    PDF(A,B;X) = B(X) if 1 <= X <= A
%                = 0    otherwise
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Henry Kucera, Winthrop Francis,
%    Computational Analysis of Present-Day American English,
%    Brown University Press, 1967.
%
%  Parameters:
%
%    Input, integer X, the word length whose probability is desired.
%
%    Output, real PDF, the value of the PDF.
%
  word_length_max = 79;

  pdf_vec = [ ...
    0.00806, ...
    0.01370, ...
    0.01862, ...
    0.02547, ...
    0.03043, ...
    0.03189, ...
    0.03516, ...
    0.03545, ...
    0.03286, ...
    0.03533, ...
    0.03562, ...
    0.03788, ...
    0.03669, ...
    0.03751, ...
    0.03518, ...
    0.03541, ...
    0.03434, ...
    0.03305, ...
    0.03329, ...
    0.03103, ...
    0.02867, ...
    0.02724, ...
    0.02647, ...
    0.02526, ...
    0.02086, ...
    0.02178, ...
    0.02128, ...
    0.01801, ...
    0.01690, ...
    0.01556, ...
    0.01512, ...
    0.01326, ...
    0.01277, ...
    0.01062, ...
    0.01051, ...
    0.00901, ...
    0.00838, ...
    0.00764, ...
    0.00683, ...
    0.00589, ...
    0.00624, ...
    0.00488, ...
    0.00477, ...
    0.00406, ...
    0.00390, ...
    0.00350, ...
    0.00318, ...
    0.00241, ...
    0.00224, ...
    0.00220, ...
    0.00262, ...
    0.00207, ...
    0.00174, ...
    0.00174, ...
    0.00128, ...
    0.00121, ...
    0.00103, ...
    0.00117, ...
    0.00124, ...
    0.00082, ...
    0.00088, ...
    0.00061, ...
    0.00061, ...
    0.00075, ...
    0.00063, ...
    0.00056, ...
    0.00052, ...
    0.00057, ...
    0.00031, ...
    0.00029, ...
    0.00021, ...
    0.00017, ...
    0.00021, ...
    0.00034, ...
    0.00031, ...
    0.00011, ...
    0.00011, ...
    0.00008, ...
    0.00006 ];
  pdf_sum = 0.99768;

  if ( 1 <= x & x <= word_length_max )
    pdf = pdf_vec(x) / pdf_sum;
  else
    pdf = 0.0;
  end

  return
end