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

function test_1(testCase)
    b = [ 1 2 3; 3 4 3];
    bb = spx.dict.MatrixOperator(b);
    [m, n] = size(bb);
    verifyEqual(testCase, [m, n], size(b));
    bbb = double(bb);
    verifyEqual(testCase, bbb, b);
    v = [1 2 3]';
    verifyEqual(testCase, b*v, bb*v);
    v = [1 2]';
    verifyEqual(testCase, b'*v, bb'*v);
    c = bb';
    verifyEqual(testCase, b', c.A);
    c = bb.';
    verifyEqual(testCase, b.', c.A);
    b2 = bb.columns_operator([1 2]);
    verifyEqual(testCase, double(b2), b(:, [1 2]));
end

function test_dct_basis(testCase)
    N = 4;
    Dict = spx.dict.DCTBasis(N);
    A = double(Dict);

    tolerance = 1e-10;
    verifyEqual(testCase, A, dctmtx(N)', 'AbsTol', tolerance);

    tolerance = 1e-3;
    x = [1 1 1 1]';
    y = Dict * x;
    verifyEqual(testCase, y, [1.9239
   -0.3827
    0.3827
    0.0761], 'AbsTol', tolerance);

    x = [1 1 1 1]';
    y = Dict.apply_transpose(x);
    verifyEqual(testCase, y, [2.0000
         0
         0
         0], 'AbsTol', tolerance);
    y = Dict.adjoint(x);
    verifyEqual(testCase, y, [2.0000
         0
         0
         0], 'AbsTol', tolerance);

end

function test_dft_basis(testCase)
    N = 4;
    Dict = spx.dict.DFTBasis(N);
    A = double(Dict);

    tolerance = 1e-10;
    verifyEqual(testCase, A, dftmtx(N)' / sqrt(N), 'AbsTol', tolerance);

    A = Dict.apply(eye(4));
    verifyEqual(testCase, A, dftmtx(N)' / sqrt(N), 'AbsTol', tolerance);

    A = Dict.adjoint(eye(4));
    verifyEqual(testCase, A, dftmtx(N) / sqrt(N), 'AbsTol', tolerance);

    tolerance = 1e-3;
    x = [1 1 1 1]';
    y = Dict * x;
    verifyEqual(testCase, y, [2
     0
     0
     0], 'AbsTol', tolerance);

    x = [1 1 1 1]';
    y = Dict.apply_transpose(x);
    verifyEqual(testCase, y, [2.0000
         0
         0
         0], 'AbsTol', tolerance);
    y = Dict.adjoint(x);
    verifyEqual(testCase, y, [2.0000
         0
         0
         0], 'AbsTol', tolerance);

end

function test_partial_dct_1(testCase)
    rng default;
    M = 2;
    N = 4;
    p = randperm(N);
    % select some rows randomly
    row_pics = sort(p(1:M));
    % make sure that the first row is there
    row_pics(1) = 1;
    col_perm = randperm(N);
    Dict = spx.dict.PartialDCT(row_pics, col_perm);
    x = ones(N, 1);
    y  = Dict.apply(x);
    z = Dict.adjoint(y);
    tolerance = 1e-3;
    verifyEqual(testCase, y, [2.0000
         0], 'AbsTol', tolerance);
    verifyEqual(testCase, z, [1.0000
    1.0000
    1.0000
    1.0000], 'AbsTol', tolerance);

    x(1) = x(1) + 2;
    y  = Dict.apply(x);
    z = Dict.adjoint(y);
    verifyEqual(testCase, y, [3.0000
    1.0000], 'AbsTol', tolerance);
    verifyEqual(testCase, z, [2.0000
    2.0000
    1.0000
    1.0000], 'AbsTol', tolerance);

end

function test_partial_dct_orthogonal(testCase)
    rng default;
    M = 32;
    N = 64;
    p = randperm(N);
    % select some rows randomly
    row_pics = sort(p(1:M));
    % make sure that the first row is there
    row_pics(1) = 1;
    col_perm = randperm(N);
    Dict = spx.dict.PartialDCT(row_pics, col_perm);
    y = eye(M);
    z = Dict.adjoint(y);
    u = Dict.apply(z);
    tolerance = 1e-8;
    verifyEqual(testCase, y, u, 'AbsTol', tolerance);
end

function test_partial_dft_1(testCase)
    rng default;
    M = 8;
    N = 16;
    p = randperm(N);
    % select some rows randomly
    row_pics = sort(p(1:M));
    % make sure that the first row is there
    row_pics(1) = 1;
    col_perm = randperm(N);
    Dict = spx.dict.PartialDFT(row_pics, col_perm);
    x = ones(N, 1);
    y  = Dict.apply(x);
    z = Dict.adjoint(y);
    tolerance = 1e-3;
    verifyEqual(testCase, y, [4
     0
     0
     0
     0
     0
     0
     0], 'AbsTol', tolerance);
    verifyEqual(testCase, z, [1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1], 'AbsTol', tolerance);

    x(1) = x(1) + 2;
    y  = Dict.apply(x);
    z = Dict.adjoint(y);
    verifyEqual(testCase, y, [4.5000 + 0.0000i
   0.5000 + 0.0000i
   0.0000 - 0.5000i
  -0.5000 + 0.0000i
   0.0000 + 0.5000i
  -0.5000 + 0.0000i
   0.0000 - 0.5000i
   0.0000 + 0.5000i], 'AbsTol', tolerance);
   verifyEqual(testCase, z, [2.0000 + 0.0000i
   0.9482 + 0.1250i
   0.9482 - 0.1250i
   1.3018 - 0.1250i
   1.0000 - 0.3536i
   1.3018 + 0.1250i
   1.0000 + 0.0000i
   1.0000 + 0.3536i
   1.3018 - 0.1250i
   1.0000 - 0.3536i
   1.3018 + 0.1250i
   0.9482 - 0.1250i
   0.9482 + 0.1250i
   1.0000 + 0.3536i
   1.0000 + 0.0000i
   1.0000 + 0.0000i], 'AbsTol', tolerance);
end

function test_partial_dft_orthogonal(testCase)
    rng default;
    M = 32;
    N = 64;
    p = randperm(N);
    % select some rows randomly
    row_pics = sort(p(1:M));
    % make sure that the first row is there
    row_pics(1) = 1;
    col_perm = randperm(N);
    Dict = spx.dict.PartialDFT(row_pics, col_perm);
    y = eye(M);
    z = Dict.adjoint(y);
    u = Dict.apply(z);
    tolerance = 1e-8;
    verifyEqual(testCase, y, u, 'AbsTol', tolerance);
end

author: indigits

repo: sparse-plex


text:
function [mv1, mv2, mv3] = specificIMCDetails(img, varargin)
%SPECIFICIMCDETAILS Ohser's Integral of Mean Curvature
%
%   this version is just for debugging.
%   It allows to extract contribution for each directions.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2005.
%

%   HISTORY

img = img~=0;

delta = [1 1 1];

% square faces of the unit cell
kr = [...
    1 2 3 4;...
    1 2 5 6;...
    1 3 5 7;...
    1 2 7 8;...
    3 5 4 6;...
    1 6 3 8;...
    2 4 5 7;...
    2 3 6 7;...
    1 5 3 8];

% triangle faces of the unit cell
kt = [...
    1 7 6;...
    2 5 8;...
    4 7 6;...
    3 5 8;...
    2 3 8;...
    4 1 7;...
    2 3 5;...
    4 1 6];

% unit surface for each cell
c1 = 0.045778;
c2 = 0.036981;
c3 = 0.035196;
c = [c1 c1 c1   c2 c2 c2   c2 c2 c2   c3 c3 c3 c3];

% distances between pixels
d1 = delta(1);
d2 = delta(2);
d3 = delta(3);
d12  = sqrt(delta(1)*delta(1) + delta(2)*delta(2));
d13  = sqrt(delta(1)*delta(1) + delta(3)*delta(3));
d23  = sqrt(delta(2)*delta(2) + delta(3)*delta(3));
s = (d12 + d23 + d13)/2;
a123 = 2*sqrt(s*(s-d12)*(s-d13)*(s-d23));

a = [d1*d2 d1*d3 d2*d3   d3*d12 d3*d12 d2*d13 d2*d13 d1*d23 d1*d23  a123 a123 a123 a123];

% compute gray-tone histogram of the image
h = grayHist(img);

mv1 = 0;
mv2 = 0;
mv3 = 0;

% for each type of configuration
for l=1:256

    v = l-1;
    b(1) = bitand(v,1)~=0;
    b(2) = bitand(v,2)~=0;
    b(3) = bitand(v,4)~=0;
    b(4) = bitand(v,8)~=0;
    b(5) = bitand(v,16)~=0;
    b(6) = bitand(v,32)~=0;
    b(7) = bitand(v,64)~=0;
    b(8) = bitand(v,128)~=0;

    % for each square face
    for nu=1:3

        b1 = b(kr(nu, 1));
        b2 = b(kr(nu, 2));
        b3 = b(kr(nu, 3));
        b4 = b(kr(nu, 4));

        s = sum([b1 b2 b3 b4]);
        if s==1
            mv1 = mv1 + h(l)*c(nu)/4/a(nu);
        elseif s==3
            mv1 = mv1 - h(l)*c(nu)/4/a(nu);
        end
    end

    % for each square diagonal face
    for nu = 4:9

        b1 = b(kr(nu, 1));
        b2 = b(kr(nu, 2));
        b3 = b(kr(nu, 3));
        b4 = b(kr(nu, 4));

        s = sum([b1 b2 b3 b4]);
        if s==1
            mv2 = mv2 + h(l)*c(nu)/4/a(nu);
        elseif s==3
            mv2 = mv2 - h(l)*c(nu)/4/a(nu);
        end
    end

    % for each triangular face
    for nu=10:13

        b1 = b(kt(nu-9, 1));
        b2 = b(kt(nu-9, 2));
        b3 = b(kt(nu-9, 3));
        b4 = b(kt(nu-5, 1));
        b5 = b(kt(nu-5, 2));
        b6 = b(kt(nu-5, 3));

        s1 = sum([b1 b2 b3]);
        s2 = sum([b4 b5 b6]);
        if s1==1
            mv3 = mv3 + h(l)*c(nu)/3/a(nu);
        end
        if s2==2
            mv3 = mv3 - h(l)*c(nu)/3/a(nu);
        end
    end
end

%imc = 4*pi*mv/sum(h(:));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function r=rot4x(theta)
    % r = rot4x(theta)
    %
    % rotx produces a 4x4 rotation matrix representing
    % a rotation by theta radians about the x axis.
    %
    %	Argument definitions:
    %
    %	theta = rotation angle in radians
    c = cos(theta);
    s = sin(theta);
    r = [1  0  0 0;
        0  c -s 0;
        0  s  c 0;
        0  0  0 1];

author: CelsoReyes


text:
classdef DesignImagVariable < handle

   properties (Access = public)
      value
      valueOld
      valueOldOld
   end

   properties (GetAccess = public, SetAccess = private)
      xLength
   end

   properties (Access = private)

   end

   methods (Access = public)

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

       function update(obj)
           obj.valueOldOld = obj.valueOld;
           obj.valueOld = obj.value;
       end

   end

   methods (Access = private)

       function init(obj,cParams)
          obj.xLength = cParams.xLength;
          obj.valueOldOld = zeros(obj.xLength,1);
          obj.valueOld    = zeros(obj.xLength,1);
          obj.value       = zeros(obj.xLength,1);
       end

   end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function CrowdDis = CrowdingDistance(PopObj,FrontNo)
%CrowdingDistance - Calculate the crowding distances of solutions front
%by front.
%
%   CD = CrowdingDistance(F) calculates the crowding distances of solutions
%   according to their objective values in F.
%
%   CD = CrowdingDistance(F,FrontNo) calculates the crowding distances of
%   solutions in each non-dominated front, where FrontNo is the front
%   numbers of solutions.
%
%   Example:
%       CrowdDis = CrowdingDistance(PopObj,FrontNo)

%------------------------------- Reference --------------------------------
% S. Kukkonen and K. Deb, Improved pruning of non-dominated solutions based
% on crowding distance for bi-objective optimization problems, Proceedings
% of the IEEE Congress on Evolutionary Computation, 2006, 1179-1186.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,M] = size(PopObj);
    if nargin < 2
        FrontNo = ones(1,N);
    end
    CrowdDis = zeros(1,N);
    Fronts   = setdiff(unique(FrontNo),inf);
    for f = 1 : length(Fronts)
        Front = find(FrontNo==Fronts(f));
        Fmax  = max(PopObj(Front,:),[],1);
        Fmin  = min(PopObj(Front,:),[],1);
        for i = 1 : M
            [~,Rank] = sortrows(PopObj(Front,i));
            CrowdDis(Front(Rank(1)))   = inf;
            CrowdDis(Front(Rank(end))) = inf;
            for j = 2 : length(Front)-1
                CrowdDis(Front(Rank(j))) = CrowdDis(Front(Rank(j)))+(PopObj(Front(Rank(j+1)),i)-PopObj(Front(Rank(j-1)),i))/(Fmax(i)-Fmin(i));
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function metabolites = formula2mets(formula)
% Takes a rxn formula and produces a list of metabolites
%
% USAGE:
%
%    metabolites = formula2mets(formula)
%
% INPUT:
%    formula:       rxn formula
%
% OUTPUT:
%    metabolites:   list of metabolites
%
% NOTE:
%    uses `parseRxnFormula.m` (Cobra Toolbox)
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is
mets = parseRxnFormula(formula);
metabolites = regexprep(mets, '\[.+\]', '');

author: opencobra

repo: cobratoolbox


text:
function x=findQuadraticFitVertex(x1,x2,param1,param2,param3)
%%FINDQUADRATICFIRVERTEX Find the minimum of a parabola that fits given
%           scalar data when presented with one of two parameterizations.
%           the first parameterization uses two function values and one
%           derivaitve value. The second parameterization uses just two
%           derivative values.
%
%INPUTS: x1, x2 The two points at which the values are taken.
% param1,param2,param3 These thee values can be given in two ways, which
%               determine their meaning. The first is:
%               f1 The value of the function at x1.
%               f2 The value of the function at x2.
%               g1 The derivative of the function at x1.
%               The other parameterization is
%               g1 The gradient at x1.
%               g2 The gradient at x2.
%               param3 is either omitted or an empty matrix is passed when
%               using the second formulation.
%
%OUTPUTS: x The interpolated critical point (minimum or maximum) of the
%           fitted parabola.
%
%For the first formulation of the problem, we are given equations with
%unknown coefficients (a,b,c):
% f1==a*x1^2+b*x1+c
% f2==a*x2^2+b*x2+c
% gl==b+2*a*x1
%Solving for (a,b,c) one gets
%a=(-f1+f2+g1*(x1-x2))/(x1-x2)^2
%b=(2*f1*x1-x1*(2*f2+g1*x1)+g1*x2^2)/(x1-x2)^2
%c=(f2*x1^2+g1*x1*(x1-x2)*x2+f1*x2*(-2*x1+x2))/(x1-x2)^2
%Substitutign into the expression for the derivative and setting the
%derivative equal to zero, one gets a value of x of
%x=x1+(g1*(x2-x1)^2)/(2*(f1-f2+g1*(x2-x1)))
%
%For the second formulation of the problem, we are given equations with
%unknown coefficients (a,b,c):
% f1==a*x1^2+b*x1+c
% g1==b+2*a*x1
% g2==b+2*a*x2
%Solving for (a,b,c) one gets
% a=(g1-g2)/(2*x1-2*x2)
% b=(g2*x1-g1*x2)/(x1-x2)
% c=f1+(x1*(-(g1+g2)*x1+2*g1*x2))/(2*(x1-x2))
%Substituting into the expression for the derivative and setting the
%derivative equal to zero, one gets a value of x of
% x=x2+g2*(x1-x2)/(g2-g1);
%which does not require the value f1.
%
%EXAMPLE:
%Here, we use the function f=3*(x-12)*(x-14)=3*x^2-78*x+504;
% x1=-10;
% x2=24;
% f1=3*x1^2-78*x1+504;
% f2=3*x2^2-78*x2+504;
% g1=2*3*x1-78;
% g2=2*3*x2-78;
% xSol1=findQuadraticFitVertex(x1,x2,f1,f2,g1)
% xSol2=findQuadraticFitVertex(x1,x2,g1,g2)
%One will find that both xSol1 and xSol2 =13, which is the minimum point of
%the parabola.
%
%February 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin>4&&~isempty(param3))
    %Given the first type of parameterization.
    f1=param1;
    f2=param2;
    g1=param3;

    diff=x2-x1;
    x=x1+(g1*diff^2)/(2*(f1-f2+g1*diff));
else
    %Given the second type of parameterization.
    g1=param1;
    g2=param2;

    x=x2+g2*(x1-x2)/(g2-g1);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [f] = fact(n)
%FACT  Vectorized Factorial function
%
%usage: f = fact(n)
%
%tested under version 5.3.1
%
%     This function computes the factorial of
%     the elements of N.
%     N can be any size but must contain
%     Real, Non-Negative, Integers.
%
%     This routine is much more robust than
%     the built in FACTORIAL function.
%
%see also: Gamma, Prod, Factorial, Binomial

%Paul Godfrey
%pgodfrey@conexant.com
%8-23-00

[row,col]=size(n);
n=n(:);

f=NaN*n;

pp=find(imag(n)==0 & real(n)>=0 & round(n)==n);
%find integer values
nn=n(pp);
ff=zeros(length(pp),1);

s=1;
p=[];
if ~isempty(nn)
   p=find(nn==0);
end
if ~isempty(p)
   ff(p)=s;
end

%upper limit here depends upon realmax
if ~isempty(nn)
for k=1:170
    s=s*k;
%empty=scalar warning
    p=find(nn==k);
    if ~isempty(p)
        ff(p)=s;
    end
end
end

p=find(nn>170);
if ~isempty(p)
    ff(p)=Inf;
end

f(pp)=ff;
f=reshape(f,row,col);

return

author: Sable

repo: mcbench-benchmarks


text:
function m = cvecrep(v,c)

% CVECREP  Column vector replicate
%
%   M = cvecrep(V, C) Replicates a Nx1 dimensional column vector V, C times to generate a
%   NxC dimensional matrix M.
%
%   See also
%   RVECREP, REPMAT
%   Copyright (c) Oregon Health & Science University (2006)
%
%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
%   academic use only (see included license file) and can be obtained from
%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
%   software should contact rebel@csee.ogi.edu for commercial licensing information.
%
%   See LICENSE (which should be part of the main toolkit distribution) for more
%   detail.

%=============================================================================================

if isempty(v)

  m = zeros(0,c);

else

  m = v(:,ones(c,1));

end

repo: ekfmonoslam


text:
% The COBRAToolbox: testReadSBML.m
%
% Purpose:
%     - reads all the sbml files in this folder and
%       checks if all the parameters are correctly written
%     - loads certain xml files, runs and FBA, and compares the
%       solution to pre-calculated values from FBAs run with .mat files
%
% Authors:
%     - Partial original file: Joseph Kang 04/07/09
%     - CI integration: Laurent Heirendt
%
% Note:
%     - The solver libraries must be included separately

global CBTDIR

%Check the requirements (a LP solver is necessary)
solvers = prepareTest('needsLP',true,'requireOneSolverOf',{'gurobi','ibm_cplex','glpk','mosek','quadMinos'});

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testReadSBML'));
cd(fileDir);

%Check if the model contains the same information

% initialize the test
cd([CBTDIR, filesep, 'test', filesep, 'models']);

% load the models
modelArr = { 'Abiotrophia_defectiva_ATCC_49176.xml', 'STM_v1.0.xml', 'iIT341.xml'};

%loop through the models
for i = 1:length(modelArr)
    %reading the models takes quite a bit of time, so only do it once for
    %all solvers.
    % output a line before launching the test for model i
    fprintf('   Testing %s ...\n', modelArr{i});

    % load the model (actually supply the full filename of the path
    % where the model is found)
    model = getDistributedModel(modelArr{i});
    for k = 1:length(solvers.LP)
        fprintf(' -- Running testReadSBML using the solver interface: %s ... ', solvers.LP{k});
        solverOK = changeCobraSolver(solvers.LP{k}, 'LP', 0);

        fprintf('   Testing loaded model ... \n');

        % set the tolerance
        tol = 1e-6;

        % define the maximum objective values calculated from pre-converted .mat files
        modelFBAf_max = [0.149475406282249; 0.477833660760744; 0.692812693473487];

        % define the minimum objective values
        modelFBAf_min = [0.0; 0.0; 0.0];

        % solve the maximisation problem
        FBA = optimizeCbModel(model, 'max');

        % test the maximisation solution
        assert(FBA.stat == 1);
        assert(abs(FBA.f - modelFBAf_max(i)) < tol);
        assert(norm(model.S * FBA.x) < tol);

        % solve the minimisation problem
        FBA = optimizeCbModel(model, 'min');

        % test the minimisation solution
        assert(FBA.stat == 1);
        assert(abs(FBA.f - modelFBAf_min(i)) < tol);
        assert(norm(model.S * FBA.x) < tol);

        % print a line for success of loop i
        fprintf(' Done.\n');

        % test that gene rules are generated correctly
        % needs testing on model with large number of genes, i.e.,
        % 'Abiotrophia_defectiva_ATCC_49176.xml'
        if strcmp(modelArr{i}, 'Abiotrophia_defectiva_ATCC_49176.xml')
            % test that rules are correctly generated, i.e. no gene partially matched
            % no indication of x(1)23 instead of x(123).
            assert(~any(~cellfun(@isempty, regexp(model.rules, '\(\d+\)\d+')))) %incorrect
        end
    end
end

% test reading COBRA models with symbols in objective reactions and multiple objective reactions
for jTest = 1:2
    if jTest == 1
        % test objective reactions with symbols
        fprintf('   Testing readSBML for models with symbols in objective reactions ...\n');
        model = createModel({'EX_a(e)'; 'EX_b(e)'}, {'Test A'; 'Test B'}, {'a[e] <=>'; 'b[e] <=>'});
        model.c = [1; 0];

    elseif jTest == 2
        % test more than one objective reactions with >1 objective reactions
        fprintf('   Testing readSBML for models with >1 objective reactions ...\n');
        model = createModel({'EX_a'; 'EX_b'}, {'Test A'; 'Test B'}, {'a[e] <=>'; 'b[e] <=>'});
        model.c = [1; -2];
    end

    model.lb = model.lb(:);
    model.ub = model.ub(:);
    model.comps = {'e'};
    model.compNames = {'ExtraCellular'};
    % add the fields outputted by readSBML
    [model.modelVersion.SBML_level, model.modelVersion.SBML_version, model.modelVersion.fbc_version] = deal(3, 1, 2);
    model.metCharges = zeros(numel(model.mets), 1);
    model.metFormulas = {'C'; 'C'};
    model.osense = -1;
    model.description = 'test_sbml_obj.xml';
    model.genes = {'gene1'; 'gene2'};
    model.rules = {''; ''};
    model.rxnGeneMat = zeros(2, 2);

    model = convertOldStyleModel(model);

    % write the model
    writeCbModel(model, 'sbml', 'test_sbml_obj');

    % read in the model
    model2 = readCbModel('test_sbml_obj.xml');
    %We are creating a few default fields in readCbModel.
    DefaultFields = {'S','b','csense','lb','ub','c','osense','rxns','mets','genes','rules','subSystems'};

    model2 = convertOldStyleModel(model2);

    modelFields = fieldnames(model);
    modelFields = setdiff(modelFields,'rxnGeneMat'); %rxnGeneMat did not contain any information and is not created by readCbModel.
    for i = 1:numel(modelFields)
        if iscell(model.(modelFields{i}))
            if ~all(cellfun(@isempty, model.(modelFields{i})))  || any(ismember(modelFields{i},DefaultFields))
                assert(isequal(model.(modelFields{i}),model2.(modelFields{i})));
            else
                assert(~isfield(model2, modelFields{i}));
            end
        else
            if isnumeric(model.(modelFields{i}))
                if ~all(isnan(model.(modelFields{i})(:))) || any(ismember(modelFields{i},DefaultFields))
                    assert(isequal(model.(modelFields{i})(~isnan(model.(modelFields{i}))),model2.(modelFields{i})(~isnan(model.(modelFields{i})))));
                    assert(isequal(~isnan(model.(modelFields{i})),~isnan(model2.(modelFields{i}))));
                else
                    assert(~isfield(model2, modelFields{i}));
                end
            else
                %Not numeric, i.e. another field.
                assert(isequal(model.(modelFields{i}),model2.(modelFields{i})));
            end
        end
    end
    fprintf(' Done.\n\n');
end
delete('test_sbml_obj.xml')

% change the directory
cd(currentDir)


text:
function v = planeIntersect(N1,N2,alpha1,alpha2)
% compute the intersection of two planes

d = dot(N1,N2);
c1 = (alpha1 - alpha2 .* d)./(1-d.^2);
c2 = (alpha2 - alpha1 .* d)./(1-d.^2);

% the line intersecting both planes is
% a + lambda b
a = c1 .* N1 + c2 .* N2;
b = cross(N1,N2);

% determine lambda such that norm(a + lambda b) = 1
p = dot(a,b) ./ norm(b).^2;
q = (norm(a).^2 - 1)./norm(b).^2;

lambda1 = - p + sqrt(p.^2 - q);
lambda2 = - p - sqrt(p.^2 - q);

% compute the points on the sphere
v = [a(:) + lambda1(:) .* b(:), a(:) + lambda2(:) .* b(:)];

end

author: mtex-toolbox

repo: mtex

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


text:
function [zCart,RCart]=ruv2CartStdRefracCubature(zRUVBiased,SR,useHalfRange,zTx,zRx,M,Ns,xi,w,ce,rE,spherCent,xMax)
%%RUV@CARTSTDREFRACCUBATURE Use cubature integration to approximate the
%              moments of measurements converted from refraction-corrupted
%              bistatic r-u-v coordinates into Cartesian coordinates using
%              a standard exponential atmospheric model. For a two-way
%              monostatic conversion, set zTx=[0;0;0]; to make the
%              transmitter and receiver collocated.
%
%INPUTS: z A 3XnumMeas matrix of numMeas vectors to convert. Each has
%          elements [r;u;v], where r is the bistatic range from the
%          transmitter to the target to the receiver, and u and v are
%          direction cosines.
%       SR The 3X3XnumMeas lower-triangular square roots of the measurement
%          covariance matrices for the measurements. If all of the matrices
%          are the same, then this can just be a single 3X3 matrix.
% useHalfRange A boolean value specifying whether the bistatic range value
%          should be divided by two. This normally comes up when operating
%          in monostatic mode, so that the range reported is a one-way
%          range. The default if this parameter is not provided (or an
%          empty matrix is provided) is false.
%      zTx The 3X1 [x;y;z] location vector of the transmitter in global
%          Cartesian coordinates.  If this parameter is omitted or an empty
%          matrix is passed, then the receiver is placed at the origin.
%      zRx The 3X1 [x;y;z] location vector of the receiver in global
%          Cartesian coordinates. If this parameter is omitted or an empty
%          matrix is passed, then the receiver is placed at the origin.
%        M A 3X3 rotation matrix to go from the alignment of the global
%          coordinate system to the local alignment fo the receiver. The z
%          vector of the local coordinate system of the receiver is the
%          pointing direction of the receiver. If this matrix is omitted,
%          then the identity matrix is used.
%       Ns The atmospheric refractivity reduced to the reference sphere.
%          Note that the refractivity is (n-1)*1e6, where n is the index
%          of refraction. The function reduceStdRefrac2Spher can be used
%          to reduce a refractivity to the surface of a reference
%          ellipsoid. This function does not allow different
%          refractivities to be used as the transmitter and receiver. If
%          this parameter is omitted or an empty matrix is passed, a
%          default value of 313 is used.
%       xi A 3 X numCubaturePoints matrix of cubature points for the
%          numeric integration. If this and the next parameter are omitted
%          or empty matrices are passed, then fifthOrderCubPoints is used
%          to generate cubature points.
%        w A numCubaturePoints X 1 vector of the weights associated with
%          the cubature points.
%       ce The optional decay constant of the exponential model. The
%          refractivity N at height h is N=Ns*exp(-ce*(h-h0)) where h0 is
%          the reference height (in this function, the height of the
%          reference ellipsoid surface is used). ce is related to the
%          change in refractivity at an elevation of 1km based on the
%          refractivity at sea level as
%          ce=log(Ns/(Ns+DeltaN))/1000;%Units of inverse meters.
%          where the change in refractivity for a change in elevation of
%          1km is DeltaN=-multConst*exp(expConst*Ns); In [1], standard
%          values for the two constants are expConst=0.005577; and
%          multConst=7.32; If ce is omitted or an empty matrix is passed,
%          the value based on the standard model is used.
% rE,spherCent The radius of the Earth to use for the spherical Earth
%           approximation used in the model and also the offset between the
%           global model and the local spherical model. It is assumed that
%           zC,zTx,and zRx are all given in the global model and will need
%           to be transformed to the local model to the used. If rE is
%           omitted or an empty matrix is passed, then the default of
%           [rE,spherCent]=osculatingSpher4LatLon(Cart2Ellipse(zRx)) is
%           used. The defaults here mean that a WGS-84 reference ellipsoid
%           is approximated by the local osculating sphere.
%     xMax This function traces the ray to a maximum displacement in the
%          local tangent plane prior (or vertically for nearly vertical
%          directions) to attempting to find a particular range. This is an
%          optional parameter specifying the maximum distance to trace the
%          ray. The default if this parameter is omitted or an empty matrix
%          is passed is 1000e3 (1000km).
%
%OUTPUTS: zCart The approximate means of the PDF of the measurements
%               in global Cartesian [x;y;z] coordinates. This is a
%               3XnumMeas matrix.
%         RCart The approximate 3X3XnumMeas covariance matrices of the
%               PDF of the Cartesian converted measurements. This is a
%               3X3XnumMeas hypermatrix.
%
%The basic cubature conversion approach is detailed in [1] and [2]. The
%standard exponential measurement model is from [3] and the model is
%discussed in more detail in the comments to ruv2CartStdRefrac.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using 3D monostatic and bistatic
%    measurements in refractive environments," IEEE Aerospace and
%    Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 54-75, Aug.
%    2014.
%[2] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[3] B. R. Bean and G. D. Thayer, CRPL Exponential Reference Atmosphere.
%    Washington, D.C.: U. S. Department of Commerce, National Bureau of
%    Standards, Oct. 1959. [Online]. Available:
%    http://digicoll.manoa.hawaii.edu/techreports/PDF/NBS4.pdf
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numMeas=size(zRUVBiased,2);

if(nargin<3||isempty(useHalfRange))
    useHalfRange=false;
end

if(nargin<4||isempty(zTx))
   zTx=zeros(3,1);
end

if(nargin<5||isempty(zRx))
   zRx=zeros(3,1);
end

if(nargin<6||isempty(M))
   M=eye(3,3);
end

if(nargin<7||isempty(Ns))
   Ns=313;
end

if(nargin<8||isempty(xi))
    [xi,w]=fifthOrderCubPoints(3);
end

if(nargin<10||isempty(ce))
    expConst=0.005577;
    multConst=7.32;

    %The change in refractivity at an elevation of 1km based on the
    %refractivity on the surface of the Earth.
    DeltaN=-multConst*exp(expConst*Ns);
    ce=log(Ns/(Ns+DeltaN))/1000;%Units of inverse meters.
end

if(nargin<11||isempty(rE))
    %Use the radius of the Earth that is the radius of the osculating
    %sphere at the location of the observed. This will be the radius used
    %in the local spherical Earth approximation for computing atmospheric
    %refraction. This uses the WGS-84 reference ellipsoid.
    [rE,spherCent]=osculatingSpher4LatLon(Cart2Ellipse(zRx));
end

if(nargin<13||isempty(xMax))
    xMax=1000e3;%1000 kilometer assumed maximum x displacement.
end

if(size(SR,3)==1)
    SR=repmat(SR,[1,1,numMeas]);
end

zCart=zeros(3,numMeas);
RCart=zeros(3,3,numMeas);
for curMeas=1:numMeas
    %Transform the cubature points to match the given Gaussian.
    cubPoints=transformCubPoints(xi,zRUVBiased(:,curMeas),SR(:,:,curMeas));

    %Convert all of the points into Cartesian space
    cubPoints=ruv2CartStdRefrac(cubPoints,useHalfRange,zTx,zRx,M,Ns,ce,rE,spherCent,xMax);

    %Extract the first two moments of the transformed points.
    [zCart(:,curMeas),RCart(:,:,curMeas)]=calcMixtureMoments(cubPoints,w);
end
end

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


text:
function c = wavecdf97(x, nlevel)
%WAVECDF97: Multi-level discrete 2-D wavelet transform
%with the Cohen-Daubechies-Feauveau (CDF) 9/7 wavelet.
%
% c = wavecdf97(x, nlevel) does the follows according to the value of
%   nlevel:
%   nlevel > 0:   decomposes 2-dimension matrix x up to nlevel level;
%   nlevel < 0:   does the inverse transform to nlevel level;
%   nlevel = 0:   sets c equal to x;
%   omitted:      does the same as nlevel=5.
%
% The boundary handling method is symmetric extension.
%
% x may be of any size; it need not have size divisible by 2^L.
%   For example, if x has length 9, one stage of decomposition
%   produces a lowpass subband of length 5 and a highpass subband
%   of length 4.  Transforms of any length have perfect
%   reconstruction (exact inversion).
%   NOTE: the 5 lines above are quoted directly form [3].
%
% If nlevel is so large that the approximation coefficients become
%   a 1-D array, any further decomposition will be performed as for 1-D
%   decomposition until the approximation coefficients be a scale number.
%
% Lifting algorithm is not used here; we use subband filters directly.
%   Lifting algorithm and spline 5/3 wavelets and other jpeg2000 related
%   codes will be available soon.
%
% Example:
%   Y = wavecdf97(X, 5);    % Decompose image X up to 5 level
%   R = wavecdf97(Y, -5);   % Reconstruct from Y
%
% You can test wavecdf97.m with the following lines:
%   % get a 2-D uint8 image
%   x=imread('E:\study\jpeg2000\images\lena.tif');
%   % decompose
%   y=wavecdf97(x,2);
%   % show decomposed result
%   figure;imshow(mat2gray(y));
%   % reconstruct without change of anything
%   ix=wavecdf97(y,-2);
%   % show and compare the original and reconstructed images
%   figure;subplot(1,2,1);imshow(x);subplot(1,2,2);imshow(uint8(ix));
%   % look at the MSE difference
%   sum(sum((double(x)-ix).^2))/numel(x)
%
% Reference:
%   [1] D.S.Taubman et al., JPEC2000 Image Compression: F. S. & P.,
%       Chinese Edition, formula 10.6-10.9 in section 10.3.1
%       and formula 10.13 in section 10.4.1.
%   [2] R.C.Gonzalez et al., Digital Image Processing Using MATLAB,
%       Chinese Edition, function wavefast in section 7.2.2.
%   [3] Pascal Getreuer, waveletcdf97.m from Matlab file Exchange website
%   [4] Matlab files: biorwavf.m, wavdec2.m, wawrec2.m, etc.
%
% Contact information:
%   Email/MSN messenger:  wangthth@hotmail.com
%
% Tianhui Wang at Beijing, China,   July, 2006
%                  Last Revision:   Aug 5, 2006

%---------------------- input arguments checking  ----------------------%
error(nargchk(1,2,nargin));
if nargin == 1
    nlevel = 5; % default level
end
% check x
if ~isreal(x) || ~isnumeric(x) || (ndims(x) > 2)
    error('WAVELIFT:InArgErr', ['The first argument must' ...
        ' be a real, numeric 2-D or 1-D matrix.']);
end
if isinteger(x)
    x = double(x);
end
% check nlevel
if ~isreal(nlevel) || ~isnumeric(nlevel) || round(nlevel)~=nlevel
    error('WAVELIFT:InArgErr', ['The 2nd argument shall be ' ...
        'a real and numeric integer.']);
end
%---------------- forming low-pass and high-pass filters ---------------%
% CDF 9/7 filters: decomposition low-pass lp and high-pass hp
%                  reconstruction low-pass lpr and high-pass hpr
% The filter coefficients have several forms.
% What D.S.Taubman et al. suggest in [1] are used here:
lp = [.026748757411 -.016864118443 -.078223266529 .266864118443];
lp = [lp .602949018236 fliplr(lp)];
hp = [.045635881557 -.028771763114 -.295635881557];
hp = [hp .557543526229 fliplr(hp)];
lpr = hp .* [-1 1 -1 1 -1 1 -1] * 2;
hpr = lp .* [1 -1 1 -1 1 -1 1 -1 1] * 2;
% Matlab 'bior4.4' use the varied version (see Matlab's biorwavf.m):
%  lp=lp*sqrt(2);hp=hp*(-sqrt(2));lpr=lpr*(1/sqrt(2));hpr=hpr*(-1/sqrt(2));
% P.Getreuer's waveletcdf97.m [3] alters the Taubman's version as follows:
%  lp=lp*sqrt(2);hp=hp*sqrt(2);lpr=lpr*(1/sqrt(2));hpr=hpr*(1/sqrt(2));
% while R.C.Gonzalez et al in [2] alter the Taubman's version as follows:
%  lp=lp;hp=hp*(-2);lpr=lpr;hpr=hpr*(-1/2);
%----------------  remain unchanged when nlevel = 0  -------------------%
if nlevel == 0
    c = x;
%--------------------  decomposition,  if nlevel < 0  ------------------%
elseif nlevel > 0
    c = zeros(size(x));
    x = double(x);
    for i = 1:nlevel
        % [ll, hl; lh, hh]: 1-level FWT for x
        temp = symconv2(x, hp, 'col');    % high filtering
        temp = temp(2:2:end, :);          % down sampling
        hh = symconv2(temp, hp, 'row');   % high filtering
        hh = hh(:, 2:2:end);              % down sampling
        lh = symconv2(temp, lp, 'row');   % low filtering
        lh = lh(:, 1:2:end);              % down sampling

        temp = symconv2(x, lp, 'col');    % low filtering
        temp = temp(1:2:end, :);          % down sampling
        hl = symconv2(temp, hp, 'row');   % high filtering
        hl = hl(:, 2:2:end);              % down sampling
        ll = symconv2(temp, lp, 'row');   % low filtering
        ll = ll(:, 1:2:end);              % down sampling
        % update coefficient matrix
        c(1:size(x,1), 1:size(x,2)) = [ll, hl; lh, hh];
        % replace x with ll for next level FWT
        x = ll;
        % give a warning if nlevel is too large
        if size(x,1)<=1 && size(x,2)<=1 && i~=nlevel
            warning('WAVECDF97:DegradeInput', ['Only decompose to ' ...
                num2str(i) '-level instead of ' num2str(nlevel) ...
                ', \nas the approximation coefficients at ' num2str(i) ...
                '-level has row or/and column of length 1.']);
            break
        end
    end
%--------------------  reconstruction,  if nlevel < 0  -----------------%
else
    sx = size(x);
    % find reconstruction level
    nl = -nlevel;
    while sx(1)/2^nl<=1/2 && sx(2)/2^nl<=1/2,  nl = nl-1;  end
    if nl ~= -nlevel
        warning('WAVECDF97:DegradeInput', ['Only reconstruct to ' ...
            num2str(nl) '-level instead of ' num2str(-nlevel) ...
            ', \nas the approximation coefficients at ' num2str(nl) ...
            '-level has row or/and column of length 1.']);
    end
    % nl-level reconstruction
    for i = 1:nl
        % find the target ll hl lh hh blocks
        sLL = ceil(sx/2^(nl-i+1));
        sConstructed = ceil(sx/2^(nl-i));
        sHH = sConstructed - sLL;
        lrow = sConstructed(1); lcol = sConstructed(2);

        ll = x(1:sLL(1), 1:sLL(2));
        hl = x(1:sLL(1), sLL(2)+1:sLL(2)+sHH(2));
        lh = x(sLL(1)+1:sLL(1)+sHH(1), 1:sLL(2));
        hh = x(sLL(1)+1:sLL(1)+sHH(1), sLL(2)+1:sLL(2)+sHH(2));

        % upsample rows and low filter
        temp = zeros(sLL(1), lcol); temp(:, 1:2:end) = ll;
        ll = symconv2(temp, lpr, 'row');
        % upsample rows and high filter
        temp = zeros(sLL(1), lcol); temp(:, 2:2:end) = hl;
        hl = symconv2(temp, hpr, 'row');
        % upsample columns and low filter
        temp = zeros(lrow, lcol); temp(1:2:end, :) = ll + hl;
        l = symconv2(temp, lpr, 'col');

        % upsample rows and high filter
        temp = zeros(sHH(1), lcol); temp(:, 1:2:end) = lh;
        lh = symconv2(temp, lpr, 'row');
        % upsample rows and high filter
        temp = zeros(sHH(1), lcol); temp(:, 2:2:end) = hh;
        hh = symconv2(temp, hpr, 'row');
        % upsample rows and high filter
        temp = zeros(lrow, lcol); temp(2:2:end, :) = lh + hh;
        h = symconv2(temp, hpr, 'col');

        % update x with the new ll, ie. l+h
        x(1:lrow, 1:lcol) = l + h;
    end
    % output
    c = x;
end
%------------------------- internal function  --------------------------%
%       2-dimension convolution with edges symmetrically extended       %
%-----------------------------------------------------------------------%
function y = symconv2(x, h, direction)
% symmetrically extended convolution(see section 6.5.2 in [1]):
%    x[n], E<=n<=F-1, is extended to x~[n] = x[n], 0<=n<=N-1;
%                                  x~[E-i] = x~[E+i], for all integer i
%                                x~[F-1-i] = x~[F-1+i], for all integer i
%    For odd-length h[n], to convolve x[n] and h[n], we just need extend x
%    by (length(h)-1)/2  for both left and right edges.
% The symmetric extension handled here is not the same as in Matlab
%  wavelets toolbox nor in [2]. The last two use the following method:
%    x[n], E<=n<=F-1, is extended to x~[n] = x[n], 0<=n<=N-1;
%                                  x~[E-i] = x~[E+i-1], for all integer i
%                                x~[F-1-i] = x~[F+i], for all integer i

l = length(h); s = size(x);
lext = (l-1)/2; % length of h is odd
h = h(:)'; % make sure h is row vector
y = x;
if strcmp(direction, 'row') % convolving along rows
    if ~isempty(x) && s(2) > 1 % unit length array skip convolution stage
        for i = 1: lext
            x = [x(:, 2*i), x, x(:, s(2)-1)]; % symmetric extension
        end
        x = conv2(x, h);
        y = x(:, l:s(2)+l-1);
    end
elseif strcmp(direction, 'col') % convolving along columns
    if ~isempty(x) && s(1) > 1 % unit length array skip convolution stage
        for i = 1: lext
            x = [x(2*i, :); x; x(s(1)-1, :)]; % symmetric extension
        end
        x = conv2(x, h');
        y = x(l:s(1)+l-1, :);
    end
end
% EOF


text:
function cl = orthviews(image_obj, varargin)
% Orthviews display (SPM) for CANlab object
%
% :Usage:
% ::
%
%    cl = orthviews(image_object)
%
%    % OR
%
%    cl = orthviews(image_object, handle_number of existing orthviews)
%
% :Features:
% - Uses cluster_orthviews.m
% - Will take inputs to cluster_orthviews
% - e.g., orthviews(pstat, 'unique', 'solid');
%
% Output is clusters structure (see also region.m)
%
% Pass in 'largest_region' to center the orthviews on the largest region in the  image
%
%
% :Example:
% ::
%
%    % T-test, Construct a stats_image object, threshold and display:
%    statsimg = ttest(fmridat, .001, 'unc');
%
%    % Re-threshold and display:
%    statsimg = threshold(statsimg, .000001, 'unc');
%    orthviews(statsimg);
%
%    statsimg = threshold(statsimg, .01, 'fdr');
%    orthviews(statsimg);
%
%    % Create an orthviews and view at multiple thresholds in different panes:
%    overlay = which('SPM8_colin27T1_seg.img');
%    spm_check_registration(repmat(overlay, n, 1));
%    statsimg = ttest(fmridat);
%    statsimg = threshold(statsimg, .001, 'unc');
%    orthviews(statsimg, 'handle', 1);
%
%    statsimg = threshold(statsimg, .000001, 'unc');
%    orthviews(statsimg, 'handle', 2);
%
% :See also: statistic_image.multi_threshold

input_handle = [];
cl = [];
doreg=0;
overlay = which('keuken_2014_enhanced_for_underlay.img');
%dounique = 0; uniquestr = 'nounique';

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'han', 'handle', 'input_handle'}, input_handle = varargin{i+1};
            case 'largest_region', doreg = 1;
            case 'overlay', overlay = varargin{i + 1}; varargin{i + 1} = [];
            %case 'unique', dounique = 1; uniquestr = 'unique';

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if isempty(input_handle)
    % re-initialize new orthviews

    n = size(image_obj.dat, 2);

    spm_check_registration(repmat(overlay, n, 1));

    handle_indices = 1:n;
else

    % use existing
    handle_indices = input_handle;

end

image_indx = 1;

if isempty(image_obj.sig), image_obj.sig = true(size(image_obj.dat)); end

% replace missing voxels if necessary
image_obj.sig = zeroinsert(image_obj.removed_voxels, image_obj.sig);
image_obj.dat = zeroinsert(image_obj.removed_voxels, image_obj.dat);

% do not do this, because it may replace deliberately removed imgs
% image_obj = replace_empty(image_obj);

for i = handle_indices

    if size(image_obj.sig, 2) < image_indx
        disp('No sig vector; displaying all voxels')
    end

    % use .sig field if we have it, otherwise use dat
    if size(image_obj.sig, 2) >= image_indx && any(image_obj.sig(:, image_indx))
        cl{i} = iimg_indx2clusters(image_obj.dat(:, image_indx) .* double(image_obj.sig(:, image_indx)), image_obj.volInfo);
    else
        cl{i} = iimg_indx2clusters(image_obj.dat(:, image_indx), image_obj.volInfo);
        %disp('No sig vector; displaying all voxels')
    end

    if ~isempty(cl{i})
        if ~isempty(image_obj.image_names) && ischar(image_obj.image_names)

            if size(image_obj.image_names, 1) < image_indx
                image_obj.image_names = char(image_obj.image_names, 'NO NAME');
            end

            fprintf('Displaying image %3.0f, %3.0f voxels: %s\n', i, sum(cat(1, cl{i}.numVox)), image_obj.image_names(image_indx, :));
        end

        cluster_orthviews(cl{i}, 'add', 'handle', i, varargin{:});

    elseif size(image_obj.image_names, 1) >= image_indx
        fprintf('Image %3.0f empty: %s\n', i, image_obj.image_names(image_indx, :));
    else
        fprintf('Image %3.0f empty\n', i);
    end

    image_indx = image_indx + 1;
end

%if ~dounique
    spm_orthviews_change_colormap([0 0 1], [1 1 0], [0 1 1], [.5 .5 .5], [1 .5 0]);
%end

% set crosshairs
if doreg
    [~,wh] = max(cat(1,cl{1}.numVox));
    spm_orthviews('Reposition', cl{1}(wh).mm_center);
else
    spm_orthviews('Reposition', [0 0 0]);
end
drawnow

% Set colormap
spm_orthviews_hotcool_colormap(image_obj.dat(:), 0);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [zValues,outputStatistics] = ...
      findEmbeddings(projections,trainingData,trainingEmbedding,parameters)
%findEmbeddings finds the optimal embedding of a data set into a previously
%found t-SNE embedding
%
%   Input variables:
%
%       projections -> N x (pcaModes x numPeriods) array of projection values
%       trainingData -> Nt x (pcaModes x numPeriods) array of wavelet
%                       amplitudes containing Nt data points
%       trainingEmbedding -> Nt x 2 array of embeddings
%       parameters -> struct containing non-default choices for parameters
%
%
%   Output variables:
%
%       zValues -> N x 2 array of embedding results
%       outputStatistics -> struct containing embedding outputs
%
%
% (C) Gordon J. Berman, 2014
%     Princeton University

    if nargin < 4
        parameters = [];
    end
    parameters = setRunParameters(parameters);

    setup_parpool(parameters.numProcessors)

    d = length(projections(1,:));
    numModes = parameters.pcaModes;
    numPeriods = parameters.numPeriods;

    if d == numModes*numPeriods

        data = projections;
        data(:) = bsxfun(@rdivide,data,sum(data,2));

        minT = 1 ./ parameters.maxF;
        maxT = 1 ./ parameters.minF;
        Ts = minT.*2.^((0:numPeriods-1).*log(maxT/minT)/(log(2)*(numPeriods-1)));
        f = fliplr(1./Ts);

    else

        fprintf(1,'Finding Wavelets\n');
        [data,f] = findWavelets(projections,numModes,parameters);
        data(:) = bsxfun(@rdivide,data,sum(data,2));

    end

    fprintf(1,'Finding Embeddings\n');
    [zValues,zCosts,zGuesses,inConvHull,meanMax,exitFlags] = ...
        findTDistributedProjections_fmin(data,trainingData,...
                                    trainingEmbedding,parameters);

    outputStatistics.zCosts = zCosts;
    outputStatistics.f = f;
    outputStatistics.numModes = numModes;
    outputStatistics.zGuesses = zGuesses;
    outputStatistics.inConvHull = inConvHull;
    outputStatistics.meanMax = meanMax;
    outputStatistics.exitFlags = exitFlags;

    if parameters.numProcessors > 1  && parameters.closeMatPool
        close_parpool
    end

author: gordonberman

repo: MotionMapper


text:
function  [M,mv,alpha,unA] = ...
           select_taylor_degree(A,b,m_max,p_max,prec,shift,bal,force_estm)
%SELECT_TAYLOR_DEGREE   Select degree of Taylor approximation.
%   [M,MV,alpha,unA] = SELECT_TAYLOR_DEGREE(A,m_max,p_max) forms a matrix M
%   for use in determining the truncated Taylor series degree in EXPMV
%   and EXPMV_TSPAN, based on parameters m_max and p_max.
%   MV is the number of matrix-vector products with A or A^* computed.

%   Reference: A. H. Al-Mohy and N. J. Higham, Computing the action of
%   the matrix exponential, with an application to exponential
%   integrators. MIMS EPrint 2010.30, The University of Manchester, 2010.

%   Awad H. Al-Mohy and Nicholas J. Higham, October 26, 2010.

if nargin < 8, force_estm = false; end
if nargin < 4 || isempty(p_max), p_max = 8; end
if nargin < 3 || isempty(m_max), m_max = 55; end

if p_max < 2 || m_max > 60 || m_max + 1 < p_max*(p_max - 1)
    error('>>> Invalid p_max or m_max.')
end
n = length(A);
if nargin < 7 || isempty(bal), bal = false; end
if bal
    [D B] = balance(A);
    if norm(B,1) < norm(A,1), A = B; end
end
if nargin < 5 || isempty(prec), prec = class(A); end
switch prec
    case 'double'
        load theta_taylor
    case 'single'
        load theta_taylor_single
    case 'half'
        load theta_taylor_half
end
if shift
    mu = trace(A)/n;
    A = A-mu*speye(n);
end
mv = 0;
if ~force_estm, normA = norm(A,1); end

if ~force_estm && normA <= 4*theta(m_max)*p_max*(p_max + 3)/(m_max*size(b,2));
% if true
    % Base choice of m on normA, not the alpha_p.
    unA = 1;
    c = normA;
    alpha = c*ones(p_max-1,1);
else
    unA = 0;
    eta = zeros(p_max,1); alpha = zeros(p_max-1,1);
    for p = 1:p_max
        [c,k] = normAm(A,p+1);
        c = c^(1/(p+1));
        mv = mv + k;
        eta(p) = c;
    end
    for p = 1:p_max-1
        alpha(p) = max(eta(p),eta(p+1));
    end
end
M = zeros(m_max,p_max-1);
for p = 2:p_max
    for m = p*(p-1)-1 : m_max
        M(m,p-1) = alpha(p-1)/theta(m);
    end
end


text:
%gohome   positions the command prompt at top level directory
%   IFISS scriptfile: DJS; 24 May 2006.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage (see readme.m)
%
cd ( '/panfs/panasas1/users/jburkardt/public_html/m_src/ifiss3.3' );

author: johannesgerer

repo: jburkardt-m


text:
function data = loadubjson(fname,varargin)
%
% data=loadubjson(fname,opt)
%    or
% data=loadubjson(fname,'param1',value1,'param2',value2,...)
%
% parse a JSON (JavaScript Object Notation) file or string
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% created on 2013/08/01
%
% $Id$
%
% input:
%      fname: input file name, if fname contains "{}" or "[]", fname
%             will be interpreted as a UBJSON string
%      opt: a struct to store parsing options, opt can be replaced by
%           a list of ('param',value) pairs - the param string is equivalent
%           to a field in opt. opt can have the following
%           fields (first in [.|.] is the default)
%
%           opt.SimplifyCell [0|1]: if set to 1, loadubjson will call cell2mat
%                         for each element of the JSON data, and group
%                         arrays based on the cell2mat rules.
%           opt.IntEndian [B|L]: specify the endianness of the integer fields
%                         in the UBJSON input data. B - Big-Endian format for
%                         integers (as required in the UBJSON specification);
%                         L - input integer fields are in Little-Endian order.
%           opt.NameIsString [0|1]: for UBJSON Specification Draft 8 or
%                         earlier versions (JSONLab 1.0 final or earlier),
%                         the "name" tag is treated as a string. To load
%                         these UBJSON data, you need to manually set this
%                         flag to 1.
%
% output:
%      dat: a cell array, where {...} blocks are converted into cell arrays,
%           and [...] are converted to arrays
%
% examples:
%      obj=struct('string','value','array',[1 2 3]);
%      ubjdata=saveubjson('obj',obj);
%      dat=loadubjson(ubjdata)
%      dat=loadubjson(['examples' filesep 'example1.ubj'])
%      dat=loadubjson(['examples' filesep 'example1.ubj'],'SimplifyCell',1)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

global pos inStr len  esc index_esc len_esc isoct arraytoken fileendian systemendian

if(regexp(fname,'[\{\}\]\[]','once'))
   string=fname;
elseif(exist(fname,'file'))
   fid = fopen(fname,'rb');
   string = fread(fid,inf,'uint8=>char')';
   fclose(fid);
else
   error('input file does not exist');
end

pos = 1; len = length(string); inStr = string;
isoct=exist('OCTAVE_VERSION','builtin');
arraytoken=find(inStr=='[' | inStr==']' | inStr=='"');
jstr=regexprep(inStr,'\\\\','  ');
escquote=regexp(jstr,'\\"');
arraytoken=sort([arraytoken escquote]);

% String delimiters and escape chars identified to improve speed:
esc = find(inStr=='"' | inStr=='\' ); % comparable to: regexp(inStr, '["\\]');
index_esc = 1; len_esc = length(esc);

opt=varargin2struct(varargin{:});
fileendian=upper(jsonopt('IntEndian','B',opt));
[os,maxelem,systemendian]=computer;

jsoncount=1;
while pos <= len
    switch(next_char)
        case '{'
            data{jsoncount} = parse_object(opt);
        case '['
            data{jsoncount} = parse_array(opt);
        otherwise
            error_pos('Outer level structure must be an object or an array');
    end
    jsoncount=jsoncount+1;
end % while

jsoncount=length(data);
if(jsoncount==1 && iscell(data))
    data=data{1};
end

%%-------------------------------------------------------------------------
function object = parse_object(varargin)
    parse_char('{');
    object = [];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1); % TODO
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        count=double(parse_number());
    end
    if next_char ~= '}'
        num=0;
        while 1
            if(jsonopt('NameIsString',0,varargin{:}))
                str = parseStr(varargin{:});
            else
                str = parse_name(varargin{:});
            end
            if isempty(str)
                error_pos('Name of value at position %d cannot be empty');
            end
            %parse_char(':');
            val = parse_value(varargin{:});
            num=num+1;
            object.(valid_field(str))=val;
            if next_char == '}' || (count>=0 && num>=count)
                break;
            end
            %parse_char(',');
        end
    end
    if(count==-1)
        parse_char('}');
    end
    if(isstruct(object))
        object=struct2jdata(object);
    end

%%-------------------------------------------------------------------------
function [cid,len]=elem_info(type)
id=strfind('iUIlLdD',type);
dataclass={'int8','uint8','int16','int32','int64','single','double'};
bytelen=[1,1,2,4,8,4,8];
if(id>0)
    cid=dataclass{id};
    len=bytelen(id);
else
    error_pos('unsupported type at position %d');
end
%%-------------------------------------------------------------------------

function [data, adv]=parse_block(type,count,varargin)
global pos inStr isoct fileendian systemendian
[cid,len]=elem_info(type);
datastr=inStr(pos:pos+len*count-1);
if(isoct)
    newdata=int8(datastr);
else
    newdata=uint8(datastr);
end
id=strfind('iUIlLdD',type);
if(fileendian~=systemendian)
    newdata=swapbytes(typecast(newdata,cid));
end
data=typecast(newdata,cid);
adv=double(len*count);

%%-------------------------------------------------------------------------

function object = parse_array(varargin) % JSON array is written in row-major order
global pos inStr
    parse_char('[');
    object = cell(0, 1);
    dim=[];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1);
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        if(next_char=='[')
            dim=parse_array(varargin{:});
            count=prod(double(dim));
        else
            count=double(parse_number());
        end
    end
    if(~isempty(type))
        if(count>=0)
            [object, adv]=parse_block(type,count,varargin{:});
            if(~isempty(dim))
                object=reshape(object,dim);
            end
            pos=pos+adv;
            return;
        else
            endpos=matching_bracket(inStr,pos);
            [cid,len]=elem_info(type);
            count=(endpos-pos)/len;
            [object, adv]=parse_block(type,count,varargin{:});
            pos=pos+adv;
            parse_char(']');
            return;
        end
    end
    if next_char ~= ']'
         while 1
            val = parse_value(varargin{:});
            object{end+1} = val;
            if next_char == ']'
                break;
            end
            %parse_char(',');
         end
    end
    if(jsonopt('SimplifyCell',0,varargin{:})==1)
      try
        oldobj=object;
        object=cell2mat(object')';
        if(iscell(oldobj) && isstruct(object) && numel(object)>1 && jsonopt('SimplifyCellArray',1,varargin{:})==0)
            object=oldobj;
        elseif(size(object,1)>1 && ismatrix(object))
            object=object';
        end
      catch
      end
    end
    if(count==-1)
        parse_char(']');
    end

%%-------------------------------------------------------------------------

function parse_char(c)
    global pos inStr len
    skip_whitespace;
    if pos > len || inStr(pos) ~= c
        error_pos(sprintf('Expected %c at position %%d', c));
    else
        pos = pos + 1;
        skip_whitespace;
    end

%%-------------------------------------------------------------------------

function c = next_char
    global pos inStr len
    skip_whitespace;
    if pos > len
        c = [];
    else
        c = inStr(pos);
    end

%%-------------------------------------------------------------------------

function skip_whitespace
    global pos inStr len
    while pos <= len && isspace(inStr(pos))
        pos = pos + 1;
    end

%%-------------------------------------------------------------------------
function str = parse_name(varargin)
    global pos inStr
    bytelen=double(parse_number());
    if(length(inStr)>=pos+bytelen-1)
        str=inStr(pos:pos+bytelen-1);
        pos=pos+bytelen;
    else
        error_pos('End of file while expecting end of name');
    end
%%-------------------------------------------------------------------------

function str = parseStr(varargin)
    global pos inStr
 % len, ns = length(inStr), keyboard
    type=inStr(pos);
    if type ~= 'S' && type ~= 'C' && type ~= 'H'
        error_pos('String starting with S expected at position %d');
    else
        pos = pos + 1;
    end
    if(type == 'C')
        str=inStr(pos);
        pos=pos+1;
        return;
    end
    bytelen=double(parse_number());
    if(length(inStr)>=pos+bytelen-1)
        str=inStr(pos:pos+bytelen-1);
        pos=pos+bytelen;
    else
        error_pos('End of file while expecting end of inStr');
    end

%%-------------------------------------------------------------------------

function num = parse_number(varargin)
    global pos inStr isoct fileendian systemendian
    id=strfind('iUIlLdD',inStr(pos));
    if(isempty(id))
        error_pos('expecting a number at position %d');
    end
    type={'int8','uint8','int16','int32','int64','single','double'};
    bytelen=[1,1,2,4,8,4,8];
    datastr=inStr(pos+1:pos+bytelen(id));
    if(isoct)
        newdata=int8(datastr);
    else
        newdata=uint8(datastr);
    end
    if(fileendian~=systemendian)
        newdata=swapbytes(typecast(newdata,type{id}));
    end
    num=typecast(newdata,type{id});
    pos = pos + bytelen(id)+1;

%%-------------------------------------------------------------------------

function val = parse_value(varargin)
    global pos inStr

    switch(inStr(pos))
        case {'S','C','H'}
            val = parseStr(varargin{:});
            return;
        case '['
            val = parse_array(varargin{:});
            return;
        case '{'
            val = parse_object(varargin{:});
            return;
        case {'i','U','I','l','L','d','D'}
            val = parse_number(varargin{:});
            return;
        case 'T'
            val = true;
            pos = pos + 1;
            return;
        case 'F'
            val = false;
            pos = pos + 1;
            return;
        case {'Z','N'}
            val = [];
            pos = pos + 1;
            return;
    end
    error_pos('Value expected at position %d');
%%-------------------------------------------------------------------------

function error_pos(msg)
    global pos inStr len
    poShow = max(min([pos-15 pos-1 pos pos+20],len),1);
    if poShow(3) == poShow(2)
        poShow(3:4) = poShow(2)+[0 -1];  % display nothing after
    end
    msg = [sprintf(msg, pos) ': ' ...
    inStr(poShow(1):poShow(2)) '<error>' inStr(poShow(3):poShow(4)) ];
    error( ['JSONparser:invalidFormat: ' msg] );

%%-------------------------------------------------------------------------

function str = valid_field(str)
global isoct
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',char(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' )))
        return;
    end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos))
            return;
        end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',str0(pos(i)))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
    %str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';

%%-------------------------------------------------------------------------
function endpos = matching_quote(str,pos)
len=length(str);
while(pos<len)
    if(str(pos)=='"')
        if(~(pos>1 && str(pos-1)=='\'))
            endpos=pos;
            return;
        end
    end
    pos=pos+1;
end
error('unmatched quotation mark');
%%-------------------------------------------------------------------------
function [endpos, e1l, e1r, maxlevel] = matching_bracket(str,pos)
global arraytoken
level=1;
maxlevel=level;
endpos=0;
bpos=arraytoken(arraytoken>=pos);
tokens=str(bpos);
len=length(tokens);
pos=1;
e1l=[];
e1r=[];
while(pos<=len)
    c=tokens(pos);
    if(c==']')
        level=level-1;
        if(isempty(e1r))
            e1r=bpos(pos);
        end
        if(level==0)
            endpos=bpos(pos);
            return
        end
    end
    if(c=='[')
        if(isempty(e1l))
            e1l=bpos(pos);
        end
        level=level+1;
        maxlevel=max(maxlevel,level);
    end
    if(c=='"')
        pos=matching_quote(tokens,pos+1);
    end
    pos=pos+1;
end
if(endpos==0)
    error('unmatched "]"');
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [densfield] = denserfocalv2(rho,theta,radius)
% Determine a density field in a stereonet type plot
%
% [densfield] = denserfocalv2(rho,theta,radius)
%
%input in polar coordinates
%rho: the distance of the points
%theta: angle of the points
%radius: radius of the countercircle, kind of grid size
%output is a matrix cartesian coordinates (x,y,density)

%get the number of events
totalev=length(rho);
rhotor=rho;

%first do the middle circle (densR=0)

%find the values lower than radius and count
indi=find(rhotor(:,1)<=radius);
counting=length(indi);

if counting>0
    densfield(1,1)=0;
    densfield(1,2)=0;
    densfield(1,3)=counting/totalev;
    rhotor(indi,1)=NaN;
else
    densfield(1,3) = NaN;
    densfield(1,1) = 0;
    densfield(1,2) = 0;
end

%set densR to start value radius
densR=radius;

%set the counters for the result matrix
j=2;

%loop for the distance
while densR<=1+radius

    %calculate stepwidth for the angle
        dalpha=2*asin(radius/(2*densR));

    %set angle to 0
        densalpha=0;

    %second loop for the angle
        while densalpha<=2*pi
            %calculate the distance between the middle of the circle and
            %the points
            distery=(rhotor.^2+densR^2-2.*rhotor.*densR.*cos(abs(densalpha-theta))).^0.5;

            %find the values lower than radius and count
            indi=find(distery(:,1)<=radius);
            counting=length(indi);

            %write values if counting>0
            if counting>0
                  densfield(j,3) = counting/totalev;
                  densfield(j,1) = densR * cos(densalpha);
                  densfield(j,2) = densR * sin(densalpha);
                  rhotor(indi,:)=NaN;
            else
                  densfield(j,3) = NaN;
                  densfield(j,1) = densR * cos(densalpha);
                  densfield(j,2) = densR * sin(densalpha);

            end

                   %increase j
                  j=j+1;
            %increase densalpha
            densalpha=densalpha+dalpha;

        end

      %increase densR
      densR=densR+radius;
end

end


text:
% Reconstruction of 2D Cartesian Pulseq data
% provides an example on how data reordering can be detected from the MR
% sequence with almost no additional prior knowledge
%
% it loads Matlab .mat files with the rawdata in the format
%     adclen x channels x readouts
% if Matlab .mat file not available it attempt to load Siemens .dat (which needs mapVBVD in the path)
% but first it seeks the accompanying .seq file with the same name to interpret
%     the data

%% Load the latest file from the specified directory
path='../IceNIH_RawSend/'; % directory to be scanned for data files
%path='/data/Dropbox/ismrm2021pulseq_liveDemo/dataLive/Vienna_7T_Siemens'; % directory to be scanned for data files

pattern='*.seq';
D=dir([path filesep pattern]);
[~,I]=sort([D(:).datenum]);
seq_file_path=[path filesep D(I(end-0)).name]; % use end-1 to reconstruct the second-last data set, etc...
                                                % or replace I(end-0) with I(1) to process the first dataset, I(2) for the second, etc...
%seq_file_path='../interpreters/siemens/data_example/gre_example.seq'

% keep basic filename without the extension
[p,n,e] = fileparts(seq_file_path);
basic_file_path=fullfile(p,n);

% try loading Matlab data
data_file_path=[basic_file_path '.mat'];
if isfile(data_file_path)
    fprintf(['loading `' data_file_path '´ ...\n']);
    data_unsorted = load(data_file_path);
    if isstruct(data_unsorted)
        fn=fieldnames(data_unsorted);
        assert(length(fn)==1); % we only expect a single variable
        data_unsorted=data_unsorted.(fn{1});
    end
else
    % revert to Siemens .dat file
    data_file_path=[basic_file_path '.dat'];
    fprintf(['loading `' data_file_path '´ ...\n']);
    twix_obj = mapVBVD(data_file_path);
    if iscell(twix_obj)
        data_unsorted = twix_obj{end}.image.unsorted();
    else
        data_unsorted = twix_obj.image.unsorted();
    end
end
[adc_len,channels,readouts]=size(data_unsorted);

%% Load sequence from file
fprintf(['loading `' seq_file_path '´ ...\n']);
seq = mr.Sequence();              % Create a new sequence object
seq.read(seq_file_path,'detectRFuse');
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();
figure; plot(ktraj(1,:),ktraj(2,:),'b',...
             ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % a 2D plot
axis('equal'); title('2D kx/ky k-space trajectory');

figure; plot(ktraj(1,:),ktraj(3,:),'b',...
             ktraj_adc(1,:),ktraj_adc(3,:),'r.'); % another 2D plot
axis('equal'); title('2D kx/kz k-space trajectory');

%% Analyze the trajectory data (ktraj_adc)
fprintf('analyzing the k-space trajectory ...\n');
k_extent=max(abs(ktraj_adc),[],2);
k_scale=max(k_extent);
k_threshold=k_scale/5000;

% detect unused dimensions and delete them
if any(k_extent<k_threshold)
    ktraj_adc(k_extent<k_threshold,:)=[]; % delete rows
    k_extent(k_extent<k_threshold)=[];
end

% detect dK, k-space reordering and repetitions (or slices, etc)
kt_sorted=sort(ktraj_adc,2);
dk_all=kt_sorted(:,2:end)-kt_sorted(:,1:(end-1));
dk_all(dk_all<k_threshold)=NaN;
dk_min=min(dk_all,[],2);
dk_max=max(dk_all,[],2);
dk_all(dk_all-dk_min(:,ones(1,size(dk_all,2)))>k_threshold)=NaN;
dk_all_cnt=sum(isfinite(dk_all),2);
dk_all(~isfinite(dk_all))=0;
dk=sum(dk_all,2)./dk_all_cnt;
dk(~isfinite(dk))=0;
[~,k0_ind]=min(sum(ktraj_adc.^2,1));
kindex=round((ktraj_adc-ktraj_adc(:,k0_ind*ones(1,size(ktraj_adc,2))))./dk(:,ones(1,size(ktraj_adc,2))));
kindex(~isfinite(kindex))=0;
kindex_min=min(kindex,[],2);
kindex_mat=kindex-kindex_min(:,ones(1,size(ktraj_adc,2)))+1;
kindex_end=max(kindex_mat,[],2);
sampler=zeros(kindex_end');
repeat=zeros(1,size(ktraj_adc,2));
for i=1:size(kindex_mat,2)
    if (size(kindex_mat,1)==3)
        ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i),kindex_mat(3,i));
    else
        ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i));
    end
    repeat(i)=sampler(ind);
    sampler(ind)=repeat(i)+1;
end
if (max(repeat(:))>0)
    kindex=[kindex;(repeat+1)];
    kindex_mat=[kindex_mat;(repeat+1)];
    kindex_end=max(kindex_mat,[],2);
end
%figure; plot(kindex(1,:),kindex(2,:),'.-');

%% sort the k-space data into the data matrix
% the incoming data order is [kx coils acquisitions]
data_coils_last = permute(data_unsorted, [1, 3, 2]);
data_coils_last = reshape(data_coils_last, [adc_len*readouts, channels]);

data=zeros([kindex_end' channels]);
if (size(kindex,1)==3)
    for i=1:size(kindex,2)
        data(kindex_mat(1,i),kindex_mat(2,i),kindex_mat(3,i),:)=data_coils_last(i,:);
    end
else
    for i=1:size(kindex,2)
        data(kindex_mat(1,i),kindex_mat(2,i),:)=data_coils_last(i,:);
    end
end

if size(kindex,1)==3
    nImages=size(data,3);
else
    nImages=1;
    data=reshape(data, [size(data,1) size(data,2) 1 size(data,3)]); % we need a dummy images/slices dimension
end

%figure; imab(log(abs(data))); title('k-space data');

%% Reconstruct coil images

images = zeros(size(data));
%figure;

for ii = 1:channels
    images(:,:,:,ii) = fftshift(fft2(fftshift(data(:,:,:,ii)))); % 1.4.0. does not need inversion of the read direction
    %for ni = 1:nImages
        %tmp = abs(images(:,:,ni,ii));
        %tmp = tmp./max(tmp(:));
        %imwrite(tmp, ['img_coil_' num2str(ii) '.png'])
    %end
end

% Phase images (possibly channel-by-channel and echo-by-echo)
%figure;imab(angle(images));colormap('jet');
%figure;imab(abs(images));colormap('gray');

%% Image display with optional sum of squares combination
% figure;
% if nCoils>1
%     sos=abs(sum(images.*conj(images),ndims(images))).^(1/2);
%     sos=sos./max(sos(:));
%     imab(sos);
%     %imwrite(sos, ['img_combined.png']
% else
%     imab(abs(images));
% end
% colormap('gray');

%% 3D image now

images3D=fftshift(fft(fftshift(images,3),[],3),3);
if channels>1
    sos3D=abs(sum(images3D.*conj(images3D),ndims(images3D))).^(1/2);
    im3D=sos3D./max(sos3D(:));
else
    im3D=abs(images3D);
end

figure;imab(im3D);colormap('gray');title('all reconstructed image(s)');
figure;imab(im3D(:,:,end/2+1));colormap('gray'); title('central partition');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [beta_gibbs omega_gibbs F_gibbs L_gibbs phi_gibbs sigma_gibbs lambda_t_gibbs sigma_t_gibbs sbar]=tvbvar2gibbs(G,sigmahat,T,chi,psi,kappa,betahat,q,n,It,Bu,I_tau,I_om,H,Xbar,y,alpha0,yt,Xbart,upsilon0,f0,delta0,gamma,pick,pickf)

% preliminary elements for the algorithm
% compute the product G'*I_gamma*G (to speed up computations of deltabar)
GIG=G'*I_om*G;
% set tau as a large value
tau=10000;
% set omega as a large value
om=5;
% compute psibar
chibar=(chi+T)/2;
% compute alphabar
kappabar=T+kappa;
% compute alphabar
alphabar=T+alpha0;

% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices and cells
beta_gibbs=[];
omega_gibbs=[];
F_gibbs=[];
L_gibbs=[];
phi_gibbs=[];
sigma_gibbs=[];
lambda_t_gibbs={};
sigma_t_gibbs={};

% step 1: determine initial values for the algorithm

% initial value for B
B=kron(ones(T,1),betahat);
% initial value Omega
omega=diag(diag(betahat*betahat'));
% invert Omega
invomega=diag(1./diag(omega));
% initial value for f_2,...,f_n
% obtain the triangular factorisation of sigmahat
[Fhat Lambdahat]=bear.triangf(sigmahat);
% obtain the inverse of Fhat
[invFhat]=bear.invltod(Fhat,n);
% create the cell storing the different vectors of invF
Finv=cell(n,1);
% store the vectors
for ii=2:n
    Finv{ii,1}=invFhat(ii,1:ii-1);
end
% initial values for L_1,...,L_n
L=zeros(T,n);
% initial values for phi_1,...,phi_n
phi=ones(1,n);
% initiate invsigmabar
invsigmabar=sparse(kron(eye(T),inv(sigmahat)));

% step 2: determine the sbar values and Lambda
sbar=diag(Lambdahat);
Lambda=sparse(diag(sbar));

% step 3: recover the series of initial values for lambda_1,...,lambda_T and sigma_1,...,sigma_T
lambda_t=repmat(diag(sbar),1,1,T);
sigma_t=repmat(sigmahat,1,1,T);

hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler');  %create the progress bar

% run the Gibbs sampler
while count<=It
    % count

    hbar.iterate(1);   % update progress by one iteration

    % step 4: draw B
    invomegabar=H'*kron(I_tau,invomega)*H+Xbar'*invsigmabar*Xbar;
    % compute the choleski of invomegabar
    C=chol(bear.nspds(invomegabar),'Lower');
    % compute temporary value
    temp=Xbar'*invsigmabar*y;
    % smoothing phase: solve by back substitution
    temp1=C\temp;
    % smoothing phase: solve by forward substitution
    Bbar=C'\temp1;
    % simulation phase:
    B=Bbar+C'\randn(q*T,1);
    % reshape
    Beta=reshape(B,q,T);

    % step 5: draw omega from its posterior
    % compute the summ
    summ=(1/tau)*Beta(:,1)*Beta(:,1)';
    for ii=2:T
        summ=summ+(Beta(:,ii)-Beta(:,ii-1))*(Beta(:,ii)-Beta(:,ii-1))';
    end
    summ=diag(summ);
    % obtain Qbar
    psibar=summ+psi;
    % draw omega
    omega=diag(arrayfun(@bear.igrandn,kron(ones(q,1),chibar),psibar));
    % invert it for next iteration
    invomega=diag(1./diag(omega));

    % step 6: draw the series f_2,...,f_n from their conditional posteriors
    % recover first the residuals
    for jj=1:T
        epst(:,:,jj)=yt(:,:,jj)-Xbart{jj,1}*Beta(:,jj);
    end
    % then draw the vectors in turn
    for jj=2:n
        % first compute the summations required for upsilonbar and fbar
        summ1=zeros(jj-1,jj-1);
        summ2=zeros(jj-1,1);
        % run the summation
        for kk=1:T
            prodt=epst(1:jj-1,1,kk)*exp(-L(kk,jj));
            summ1=summ1+prodt*epst(1:jj-1,1,kk)';
            summ2=summ2+prodt*epst(jj,1,kk)';
        end
        summ1=(1/sbar(jj,1))*summ1;
        summ2=(-1/sbar(jj,1))*summ2;
        % then obtain the inverse of upsilon0
        invupsilon0=diag(1./diag(upsilon0{jj,1}));
        % obtain upsilonbar
        invupsilonbar=summ1+invupsilon0;
        C=chol(bear.nspd(invupsilonbar));
        invC=C\speye(jj-1);
        upsilonbar=full(invC*invC');
        % recover fbar
        fbar=upsilonbar*(summ2+invupsilon0*f0{jj,1});
        % finally draw f_i^(-1)
        Finv{jj,1}=fbar+chol(bear.nspd(upsilonbar),'lower')*randn(jj-1,1);
    end
    % recover the inverse of F
    invF=eye(n);
    for jj=2:n
        invF(jj,1:jj-1)=Finv{jj,1};
    end
    % eventually recover F
    F=bear.invltod(invF,n);
    % then update sigma
    sigma=F*Lambda*F';

    % step 7: draw the series phi_1,...,phi_n from their conditional posteriors
    % draw the parameters in turn
    for jj=1:n
        % estimate deltabar
        deltabar=L(:,jj)'*GIG*L(:,jj)+delta0;
        % draw the value phi_i
        phi(1,jj)=bear.igrandn(alphabar/2,deltabar/2);
    end

    % step 8: draw the series lambda_i,t from their conditional posteriors, i=1,...,n and t=1,...,T
    % consider variables in turn
    for jj=1:n
        % consider periods in turn
        for kk=1:T
            % a candidate value will be drawn from N(lambdabar,phibar)
            % the definitions of lambdabar and phibar varies with the period, thus define them first
            % if the period is the first period
            if kk==1
                lambdabar=(gamma*L(2,jj))/(1/om+gamma^2);
                phibar=phi(1,jj)/(1/om+gamma^2);
                % if the period is the final period
            elseif kk==T
                lambdabar=gamma*L(T-1,jj);
                phibar=phi(1,jj);
                % if the period is any period in-between
            else
                lambdabar=(gamma/(1+gamma^2))*(L(kk-1,jj)+L(kk+1,jj));
                phibar=phi(1,jj)/(1+gamma^2);
            end
            % now draw the candidate
            cand=lambdabar+phibar^0.5*randn;
            % compute the acceptance probability
            prob=bear.mhprob2(jj,cand,L(kk,jj),sbar(jj,1),epst(:,1,kk),Finv{jj,1});
            % draw a uniform random number
            draw=rand;
            % keep the candidate if the draw value is lower than the prob
            if draw<=prob
                L(kk,jj)=cand;
                % if not, just keep the former value
            end
        end
    end
    % then recover the series of matrices lambda_t and sigma_t
    for jj=1:T
        lambda_t(:,:,jj)=diag(sbar).*diag(exp(L(jj,:)));
        sigma_t(:,:,jj)=F*lambda_t(:,:,jj)*F';
    end

    % record phase
    % if the burn-in sample phase is not yet over
    if count<=Bu
        % simply add 1 to the iteration count
        count=count+1;
        % on the other hand, if the burn-in sample phase is over
    elseif count>Bu
        % adding one iteration to the count will depend on wether post-burn selection applies
        % if there is no post burn selection
        if pick==0
            % record the results
            beta_gibbs(:,count-Bu)=B;
            omega_gibbs(:,count-Bu)=diag(omega);
            F_gibbs(:,:,count-Bu)=F;
            L_gibbs(:,:,count-Bu)=L;
            phi_gibbs(count-Bu,:)=phi;
            sigma_gibbs(:,count-Bu)=sigma(:);
            for jj=1:T
                lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
            end
            % then add one to the count
            count=count+1;
            % if there is post burn selection, only one draw over 'fpick' draws will be retained
        elseif pick==1
            % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1, and do record the acceptance rate of the Metropolis-Hastings step
            if pickcount~=pickf
                pickcount=pickcount+1;
                % on the other hand, if the iteration does correspond to fpick
            elseif pickcount==pickf
                % do record the results
                beta_gibbs(:,count-Bu)=B;
                omega_gibbs(:,count-Bu)=diag(omega);
                F_gibbs(:,:,count-Bu)=F;
                L_gibbs(:,:,count-Bu)=L;
                phi_gibbs(count-Bu,:)=phi;
                sigma_gibbs(:,count-Bu)=sigma(:);
                for jj=1:T
                    lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                    sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
                end
                % then increase the regular count by 1 and re-initialise pickcount
                count=count+1;
                pickcount=1;
            end
        end
    end
end
close(hbar);   %close progress bar

% turn beta_gibbs into cell
beta_gibbs=mat2cell(beta_gibbs,repmat(q,T,1),It-Bu);

author: european-central-bank

repo: BEAR-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SpharmMatUtilImportObjs(confs, objs)

numSbj = length(objs);
if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end
if (~exist([confs.OutDirectory '/Logs'],'dir'))
    mkdir([confs.OutDirectory '/Logs']);
end

resampleFactor = confs.ResampleFactor;
%class(confs) - struct
%confs.vars -list of variable
%class(objs) - cell

h = waitbar(0,'Please wait...');
for i = 1:numSbj
    file = objs{i};
    [path, name, ext] = fileparts(file);
    diary(fullfile([confs.OutDirectory '/Logs'],[date, '_', num2str(round(cputime)), '_' name '_importObj.log']));

    vxsize = [1 1 1];
    switch upper(ext)
        case '.STL'
            [vertices,faces] = readSTL(file);
            origin = [0 0 0];
            new_name = [confs.OutDirectory '/' name '_obj.mat'];
            if exist(new_name,'file')
                prompt = {'Enter new filename:'};
                dlg_title = 'New File Name';
                num_lines = 1;
                def = {new_name};
                answer = inputdlg(prompt,dlg_title,num_lines,def);
                new_name = answer{1};
            end
            save(new_name, 'vertices', 'faces', 'origin');
        case '.HDR'
            [bim, origin] = NIfTI2bim(file, resampleFactor);
            new_name = [confs.OutDirectory '/' name '_bim.mat'];
            if exist(new_name,'file')
                prompt = {'Enter new filename:'};
                dlg_title = 'New File Name';
                num_lines = 1;
                def = {new_name};
                answer = inputdlg(prompt,dlg_title,num_lines,def);
                new_name = answer{1};
            end
            save(new_name, 'bim', 'origin', 'vxsize');
        case '.NII'
            [bim, origin] = NIfTI2bim(file, resampleFactor);
            new_name = [confs.OutDirectory '/' name '_bim.mat'];
            if exist(new_name,'file')
                prompt = {'Enter new filename:'};
                dlg_title = 'New File Name';
                num_lines = 1;
                def = {new_name};
                answer = inputdlg(prompt,dlg_title,num_lines,def);
                new_name = answer{1};
            end
            save(new_name, 'bim', 'origin', 'vxsize');
        case '.M'
            [vertices,faces] = readM(file);
            origin = [0 0 0];
            new_name = [confs.OutDirectory '/' name '_obj.mat'];
            if exist(new_name,'file')
                prompt = {'Enter new filename:'};
                dlg_title = 'New File Name';
                num_lines = 1;
                def = {new_name};
                answer = inputdlg(prompt,dlg_title,num_lines,def);
                new_name = answer{1};
            end
            save(new_name, 'vertices', 'faces', 'origin');
    end
    diary('off');
    waitbar(i/numSbj)
end
close(h);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout = transpose(varargin) %#ok<STOUT,*INUSD>
%TRANSPOSE   DELTAFUN objects are not transposable.

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

error('CHEBFUN:DELTAFUN:transpose:notPossible', ...
    'DELTAFUN objects are not transposable.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Test Problems for SEDUMI Interface
clc
clear

%Check SEDUMI is available
if(~checkSolver('sedumi',0))
    fprintf(2,'\n\nSEDUMI IS NOT AVAILABLE, NO SEDUMI TESTS CAN BE RUN\n\n');
    return;
end

files = {'sdp_truss1.mat','sdp_truss2.mat','sdp_truss3.mat','sdp_truss4.mat','sdp_hinf01.mat',...
         'sdp_hinf02.mat','sdp_hinf03.mat','sdp_hinf04.mat','sdp_arch0.mat','sdp_arch2.mat',...
         'sdp_arch4.mat','sdp_qap05.mat','sdp_qap06.mat','sdp_qap07.mat','sdp_qap08.mat',...
         'sdp_mcp100.mat','sdp_mcp124-1.mat','sdp_mcp124-2.mat','sdp_mcp124-3.mat','sdp_mcp124-4.mat',...
         'sdp_mcp250-1.mat','sdp_mcp250-2.mat','sdp_mcp250-3.mat','sdp_mcp250-4.mat'};

len = length(files);
cfval = zeros(len,1);
sfval = zeros(len,1);
optsC = optiset('solver','csdp','maxtime',15);
optsS = optiset('solver','sedumi','maxtime',15);

fprintf('SEDUMI SDP Checking\n');
%For each problem
for i = 1:len
    %Read In Problem
    prob = sdpRead(files{i});
    %Solve
    fprintf('Running Problem %s....',files{i});
    [~,cfval(i)] = solve(opti(prob,optsC));
    [~,sfval(i)] = solve(opti(prob,optsS));
    fprintf(' Done\n');
end

%Display fvals
display([cfval sfval])
if(norm(cfval-sfval) > 1e-1)
    error('Failed SeDuMi Checking');
else
    fprintf('Passed SeDuMi Checking\n');
end

% %% Complex Checking
% % Create a positive (semi)definite matrix:
% R=randn(5)+j*randn(5); R=R*R';
% % Make sure it's Hermitian:
% R=(R+R')/2;
%
% % Primal formulation
% clear At b c K
% c=vec(eye(5));
% At=vec(R);
% b=1;
% K.s=5;
% K.scomplex=1;
%
% [x,y,info]=sedumi(At,b,c,K);
%
% % Normed eigenvector:
% e=x(1:5)/norm(x(1:5));
% % Verify the eigenvalue
% [real(e'*R*e), max(eig(R))]


text:
function self = reduce(self, waveType, sourcelat, sourcelon, stationlat, stationlon, varargin)
    % s.reduce('waveType', 'surface', 'waveSpeed', 2000, 'f', 2.0, );
    % s.distance and waveSpeed assumed to be in metres (m)
    % (INPUT) s.data assumed to be in nm or Pa
    % (OUTPUT) s.data in cm^2 or Pa.m

    % REQUIRES EXTRA PARAMETERS FOR RSAM OBJECTS
       %reduced = struct('Q', Inf, 'sourcelat', NaN, 'sourcelon', NaN, 'distance', NaN, 'waveType', '', 'isReduced', false, 'f', NaN, 'waveSpeed', NaN, 'stationlat', NaN, 'stationlon', NaN);
       %use = true;
%   REDUCED:    a structure that is set is data are "reduced", i.e. corrected
%               for geometric spreading (and possibly attenuation)
%               Has 4 fields:
%                   REDUCED.Q = the value of Q used to reduce the data
%                   (Inf by default, which indicates no attenuation)
%                   REDUCED.SOURCELAT = the latitude used for reducing the data
%                   REDUCED.SOURCELON = the longitude used for reducing the data
%                   REDUCED.STATIONLAT = the station latitude
%                   REDUCED.STATIONLON = the station longitude
%                   REDUCED.DISTANCE = the distance between source and
%                   station in km
%                   REDUCED.WAVETYPE = the wave type (body or surface)
%                   assumed
%                   REDUCED.F = the frequency used for surface waves
%                   REDUCED.WAVESPEED = the S wave speed
%                   REDUCED.ISREDUCED = True if the data are reduced
%   UNITS:  the units of the data, e.g. nm / sec.
%   USE: use this rsam object in plots?
    p = inputParser;
    p.addParameter('waveSpeed', 2000);
    p.addParamter('f', 2.0);
    p.parse(varargin{:});

    self.reduced.waveSpeed = p.Results.waveSpeed;

    if self.reduced.isReduced
        disp('Data are already reduced');
        return;
    end

    self.reduced.distance = deg2km(distance(sourcelat, sourcelon, stationlat, stationlon)) *1000; % m

    switch self.units
        case 'nm'  % Displacement
            % Do computation in cm
            self.data = self.data / 1e7;
            r = self.reduced.distance * 100; % cm
            ws = waveSpeed * 100; % cm/2
            self.measure = sprintf('%sR%s',self.measure(1),self.measure(2:end));
            switch self.reduced.waveType
                case 'body'
                    self.data = self.data * r; % cm^2
                    self.units = 'cm^2';
                case 'surface'
                    wavelength = ws / p.Results.f; % cm
                    try
                            self.data = self.data .* sqrt(r * wavelength); % cm^2
                    catch
                            debug.print_debug(5, 'mean wavelength instead')
                            self.data = self.data * sqrt(r * mean(wavelength)); % cm^2
                    end
                    self.units = 'cm^2';
                    self.reduced.isReduced = true;
                otherwise
                    error(sprintf('Wave type %s not recognised'), self.reduced.waveType);
            end
        case 'Pa'  % Pressure
            % Do computation in metres
            self.data = self.data * self.reduced.distance; % Pa.m
            self.units = 'Pa m';
            self.reduced.isReduced = true;
            self.measure = sprintf('%sR%s',self.measure(1),self.measure(2:end));
        otherwise
            error(sprintf('Units %s for measure %s not recognised', self.units, self.measure));
    end
    self.reduced.sourcelat = sourcelat;
    self.reduced.sourcelon = sourcelon;
    self.reduced.stationlat = stationlat;
    self.reduced.stationlon = stationlon;

end

author: geoscience-community-codes

repo: GISMO


text:
function tests = test_stackRotate90(varargin)
%TEST_STACKROTATE90  One-line description here, please.
%   output = test_stackRotate90(input)
%
%   Example
%   test_stackRotate90
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-05-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function testRotateX(testCase)

img = createBasicTestImage;

rotX1th = cat(3, [5 6;1 2], [7 8;3 4]);
imgRX1 = stackRotate90(img, 1, 1);
assertEqual(testCase, rotX1th, imgRX1);
imgRX1 = stackRotate90(img, 'x', 1);
assertEqual(testCase, rotX1th, imgRX1);

rotX2th = cat(3, [7 8;5 6], [3 4;1 2]);
imgRX2 = stackRotate90(img, 1, 2);
assertEqual(testCase, rotX2th, imgRX2);
imgRX2 = stackRotate90(img, 'x', 2);
assertEqual(testCase, rotX2th, imgRX2);

rotX3th = cat(3, [3 4;7 8], [1 2;5 6]);
imgRX3 = stackRotate90(img, 1, 3);
assertEqual(testCase, rotX3th, imgRX3);
imgRX3 = stackRotate90(img, 'x', 3);
assertEqual(testCase, rotX3th, imgRX3);

function testRotateY(testCase)

img = createBasicTestImage;

rotY1th = cat(3, [2 6;4 8], [1 5;3 7]);
imgRY1 = stackRotate90(img, 2, 1);
assertEqual(testCase, rotY1th, imgRY1);
imgRY1 = stackRotate90(img, 'y', 1);
assertEqual(testCase, rotY1th, imgRY1);

rotY2th = cat(3, [6 5;8 7], [2 1;4 3]);
imgRY2 = stackRotate90(img, 2, 2);
assertEqual(testCase, rotY2th, imgRY2);
imgRY2 = stackRotate90(img, 'y', 2);
assertEqual(testCase, rotY2th, imgRY2);

rotY3th = cat(3, [5 1;7 3], [6 2;8 4]);
imgRY3 = stackRotate90(img, 2, 3);
assertEqual(testCase, rotY3th, imgRY3);
imgRY3 = stackRotate90(img, 'y', 3);
assertEqual(testCase, rotY3th, imgRY3);

function testRotateZ(testCase)

img = createBasicTestImage;

rotZ1th = cat(3, [3 1;4 2], [7 5;8 6]);
imgRZ1 = stackRotate90(img, 3, 1);
assertEqual(testCase, rotZ1th, imgRZ1);
imgRZ1 = stackRotate90(img, 'z', 1);
assertEqual(testCase, rotZ1th, imgRZ1);

rotZ2th = cat(3, [4 3;2 1], [8 7;6 5]);
imgRZ2 = stackRotate90(img, 3, 2);
assertEqual(testCase, rotZ2th, imgRZ2);
imgRZ2 = stackRotate90(img, 'z', 2);
assertEqual(testCase, rotZ2th, imgRZ2);

rotZ3th = cat(3, [2 4;1 3], [6 8;5 7]);
imgRZ3 = stackRotate90(img, 3, 3);
assertEqual(testCase, rotZ3th, imgRZ3);
imgRZ3 = stackRotate90(img, 'z', 3);
assertEqual(testCase, rotZ3th, imgRZ3);

function testRotateGrayscaleY(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around Y-axis
img2 = stackRotate90(img, 2, 1);

% check dimension
assertEqual(testCase, dim([1 3 2]), size(img2));

function testRotateGrayscaleX(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around X-axis
img2 = stackRotate90(img, 1, 1);

assertEqual(testCase, dim([3 2 1]), size(img2));

function testRotateGrayscaleZ(testCase)

img = createTestImageGray;
dim = size(img);

% rotate around Z-axis
img2 = stackRotate90(img, 3, 1);
assertEqual(testCase, dim([2 1 3]), size(img2));

function testColorImage(testCase)

% Create rgb image with non equal size in each dimension
lx = 1:50;
ly = 1:52;
lz = 1:54;
r = uint8(discreteBall(lx, ly, lz, [20 30 30], 15)*255);
g = uint8(discreteBall(lx, ly, lz, [30 20 30], 15)*255);
b = uint8(discreteBall(lx, ly, lz, [30 30 20], 15)*255);
img = imMergeChannels(r, g, b);

% basic checks
assertEqual(testCase, 'uint8', class(img));
assertEqual(testCase, [52 50 3 54], size(img));

% dimension of base image, without color
dim = size(r);

% rotate around Y-axis
img2 = stackRotate90(img, 2, 1);
assertEqual(testCase, [dim(1) dim(3) 3 dim(2)], size(img2));
assertEqual(testCase, img(1,1,:,1), img2(1, end,:,1));
assertEqual(testCase, img(1,1,:,end), img2(1, 1,:,1));

% rotate around X-axis
img2 = stackRotate90(img, 1, 1);
assertEqual(testCase, [dim(3) dim(2) 3 dim(1)], size(img2));

% rotate around Z-axis
img2 = stackRotate90(img, 3, 1);
assertEqual(testCase, [dim(2) dim(1) 3 dim(3)], size(img2));

function img = createBasicTestImage
img = cat(3, [1 2;3 4], [5 6;7 8]);

function img = createTestImageGray
% create a 3D gray-scale test image
%  ______
% |1    2|\
% | 5    |6|
% |3____4| |
%  \7____\8|
%
lx = 1:10;
ly = 30:5:60;
lz = 50:10:200;
[x, y, z] = meshgrid(lx, ly, lz);
img = uint8(x+y+z);

author: mattools

repo: matImage


text:
function outpoints = icbm_spm2tal(inpoints)
%
% This function converts coordinates from MNI space (normalized
% using the SPM software package) to Talairach space using the
% icbm2tal transform developed and validated by Jack Lancaster
% at the Research Imaging Center in San Antonio, Texas.
%
% http://www3.interscience.wiley.com/cgi-bin/abstract/114104479/ABSTRACT
%
% FORMAT outpoints = icbm_spm2tal(inpoints)
% Where inpoints is N by 3 or 3 by N matrix of coordinates
% (N being the number of points)
%
% ric.uthscsa.edu 3/14/07

% find which dimensions are of size 3
dimdim = find(size(inpoints) == 3);
if isempty(dimdim)
  error('input must be a N by 3 or 3 by N matrix')
end

% 3x3 matrices are ambiguous
% default to coordinates within a row
if dimdim == [1 2]
  disp('input is an ambiguous 3 by 3 matrix')
  disp('assuming coordinates are row vectors')
  dimdim = 2;
end

% transpose if necessary
if dimdim == 2
  inpoints = inpoints';
end

% Transformation matrices, different for each software package
icbm_spm = [0.9254 0.0024 -0.0118 -1.0207
	   	   -0.0048 0.9316 -0.0871 -1.7667
            0.0152 0.0883  0.8924  4.0926
            0.0000 0.0000  0.0000  1.0000];

% apply the transformation matrix
inpoints = [inpoints; ones(1, size(inpoints, 2))];
inpoints = icbm_spm * inpoints;

% format the outpoints, transpose if necessary
outpoints = inpoints(1:3, :);
if dimdim == 2
  outpoints = outpoints';
end

repo: REST


text:
function [ abd, info ] = cpbfa ( abd, lda, n, m )

%*****************************************************************************80
%
%% CPBFA factors a complex hermitian positive definite band matrix.
%
%  Discussion:
%
%    CPBFA is usually called by CPBCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Band storage:
%
%    If A is a hermitian positive definite band matrix,
%    the following program segment will set up the input.
%
%      m = (band width above diagonal)
%      do j = 1, n
%        i1 = max ( 1, j-m )
%        do i = i1, j
%          k = i-j+m+1
%          abd(k,j) = a(i,j)
%        end do
%      end do
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex ABD(LDA,N); the matrix to be factored.
%    The columns of the upper triangle are stored in the columns of ABD
%    and the diagonals of the upper triangle are stored in the rows of ABD.
%
%    Input, integer LDA, the leading dimension of ABD.
%    LDA must be at least M+1.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer M, the number of diagonals above the main diagonal.
%    0 <= M < N.
%
%    Output, integer INFO.
%    0, for normal return.
%    K, if the leading minor of order K is not positive definite.
%
%    Output, complex ABD(LDA,N); an upper triangular matrix R, stored in
%    band form, so that A = hermitian(R)*R.
%
  info = 0;

  for j = 1 : n

    s = 0.0;
    ik = m + 1;
    jk = max ( j - m, 1 );
    mu = max ( m + 2 - j, 1 );

    for k = mu : m
      t = abd(k,j) - abd(ik:ik+k-mu-1,jk)' * abd(mu:mu+k-mu-1,j);
      t = t / abd(m+1,jk);
      abd(k,j) = t;
      s = s + real ( t * conj ( t ) );
      ik = ik - 1;
      jk = jk + 1;
    end

    s = real ( abd(m+1,j) ) - s;

    if ( s <= 0.0 | imag ( abd(m+1,j) ) ~= 0.0 )
      info = j;
      break
    end

    abd(m+1,j) = sqrt ( s );

  end

  return
end