text:
function MeshStruct = read_gmsh(fileName)

T=txtfile2cell(fileName);

start_physicals =  find(contains(T,'$PhysicalNames'));
start_entities = find(contains(T,'$Entities'));
start_nodes = find(contains(T,'$Nodes'));
start_elements = find(contains(T,'$Elements'));

entity_info = str2num(T{start_entities+1});
node_info = str2num(T{start_nodes+1});
element_info = str2num(T{start_elements+1});

num_nodes = node_info(2);
num_node_blocks = node_info(1);
num_elements = element_info(2);
num_element_blocks = element_info(1);
num_entities = sum(entity_info);

%%%read physicals
num_physicals = str2num(T{start_physicals+1});
Physicals = T(start_physicals+2:start_physicals+1+num_physicals);
for i = 1:size(Physicals, 1)
    Physicals{i} = strsplit(Physicals{i},' ');
    Physicals{i}{1} = str2num(Physicals{i}{1});
    Physicals{i}{2} = str2num(Physicals{i}{2});
end

%%%read and sort entities
Entities = cell(4,1);
Points = {};Curves = {};Surfaces = {};Volumes = {};
line = start_entities+2;
for i = 1:entity_info(1)
    Temp = str2num(T{line});
    Points{i,1} = Temp;
    line = line+1;
end
for i = 1:entity_info(2)
    Temp = str2num(T{line});
    Curves{i,1} = Temp;
    line = line+1;
end
for i = 1:entity_info(3)
    Temp =  str2num(T{line});
    Surfaces{i,1} = Temp;
    line = line+1;
end
for i = 1:entity_info(4)
    Temp = str2num(T{line});
    Volumes{i,1} = Temp;
    line = line+1;
end
Entities{1} = Points;
Entities{2} = Curves;
Entities{3} = Surfaces;
Entities{4} = Volumes;

%%%read and sort nodes
node_block_info = zeros(num_node_blocks,4);
node_block_numbers = cell(num_node_blocks,1);
node_block_coords = cell(num_node_blocks,1);
line = start_nodes+2;
for i = 1:num_node_blocks

    node_block_info_i = str2num(T{line});
    node_block_size= node_block_info_i(1,4);

    node_block_info(i,:) = node_block_info_i;
    node_block_numbers_temp = T(line+1:line+node_block_size);
    node_block_coords_temp = T(line+node_block_size +1:line+node_block_size+node_block_size);
    line = line+2*node_block_size+1;

    for j = 1:node_block_size
        node_block_numbers{i}(j,:) = str2num(node_block_numbers_temp{j});
        node_block_coords{i}(j,:) = str2num(node_block_coords_temp{j});
    end

end

%%%Read and sort elements
element_block_info = zeros(num_element_blocks,4);
element_block_numbers = cell(num_element_blocks,1);
line = start_elements+2;
for i = 1:num_element_blocks

    element_block_info_i = str2num(T{line});element_block_size= element_block_info_i(1,4);
    element_block_info(i,:) = element_block_info_i;
    element_block_numbers_temp = T(line+1:line+element_block_size);

    for j = 1:element_block_size
        element_block_numbers{i}(j,:) = str2num(element_block_numbers_temp{j});
    end

    line = line+element_block_size+1;
end

node_ids = [];
nodes = [];
for i = 1:num_node_blocks

    nodes = [nodes;node_block_coords{i}];
    node_ids = [node_ids;node_block_numbers{i}];

end

nodes = double(nodes);
Mesh.Physicals = Physicals;
Mesh.entity_info = entity_info;
Mesh.node_info = node_info;
Mesh.element_info = element_info;
Mesh.nodes = nodes;
Mesh.Node_Blocks = node_block_numbers;
Mesh.Node_ids = node_ids;
Mesh.node_block_info = node_block_info;
Mesh.Entities = Entities;
Mesh.Elements = element_block_numbers;
Mesh.element_block_info = element_block_info;

%%
vol_ct = 0;
surf_ct = 0;
curve_ct= 0;

[~,part_name,~] = fileparts(fileName);

 node_ids = Mesh.Node_ids; %The node id's
 nodes = Mesh.nodes; %The nodel coordinates

 %%

for j = 1:size(Mesh.Physicals,1)

    physical = Mesh.Physicals{j};
    phys_name = physical{3}(2:end-1);
    dim = physical{1};
    tag = physical{2};

    switch dim
        case 3

            vol_ct = vol_ct+1;
            entities = Mesh.Entities{4,1};
            element_ids = [];
            element_mat = [];

            for k = 1:size(entities, 1)
                num_phys = entities{k}(8);
                phys_tags = entities{k}(9:9+num_phys-1);
                if sum(find(phys_tags==tag))
                   entity_tag = entities{k}(1);
                   elements_3 = (Mesh.element_block_info(:,1) == dim).*Mesh.element_block_info;
                   entity_ind = find(elements_3(:,2)==entity_tag);
                   element_ids = [element_ids;Mesh.Elements{entity_ind}(:,1)];
                   element_mat = [element_mat;Mesh.Elements{entity_ind}(:,2:end)];
                end

            end

            %febio_spec.Mesh.Elements{vol_ct}.ATTR.type='tet4'; %Element type
            volumes_names{vol_ct,1} = phys_name; %Name of this part
            elements_IDs_blocks{vol_ct} = element_ids; %Element id's
            elements_blocks{vol_ct} = element_mat; %The element matrix

        case 2

            surf_ct = surf_ct+1;
            entities = Mesh.Entities{3,1};
            element_mat = [];
            element_start = 0;

            for k = 1:size(entities, 1)
                num_phys = entities{k}(8);
                phys_tags = entities{k}(9:9+num_phys-1);
                if sum(find(phys_tags==tag))
                   entity_tag = entities{k}(1);
                   elements_2 = (Mesh.element_block_info(:,1) == dim).*(Mesh.element_block_info(:,2) == entity_tag);
                   entity_ind = find(elements_2);
                   element_mat = [element_mat;Mesh.Elements{entity_ind}(:,2:end)];
                   element_start = element_start+1;
                   element_start = element_start+size(element_mat,1);
                end

            end

            element_ids = [1:size(element_mat,1)]';
            facesBoundary_Names{surf_ct,1} = phys_name; %Name of this surface
            facesBoundary_IDs_Block{surf_ct} = element_ids; %Element id's
            facesBoundary_Block{surf_ct} =element_mat; %The element matrix

            case 1

            curve_ct = curve_ct+1;
            entities = Mesh.Entities{2,1};
            element_mat = [];
            element_start = 0;

            for k = 1:size(entities, 1)
                num_phys = entities{k}(8);
                phys_tags = entities{k}(9:9+num_phys-1);
                if sum(find(phys_tags==tag))
                   entity_tag = entities{k}(1);
                   elements_1 = (Mesh.element_block_info(:,1) == dim).*(Mesh.element_block_info(:,2) == entity_tag);
                   entity_ind = find(elements_1);
                   element_mat = [element_mat;Mesh.Elements{entity_ind}(:,2:end)];
                   element_start = element_start+1;
                   element_start = element_start+size(element_mat,1);
                end

            end

             element_ids_block{curve_ct} = element_mat(:,1);
             curve_names{curve_ct} = phys_name; %Name of this curve
             curve_elements_block{curve_ct} = element_ids; %Element id's

   end
end

%%
facesBoundary = [];
boundaryMarker = [];
elements = [];
elementsMaterialID = [];

%%
for i = 1:vol_ct

    elementMaterialID_block = ones(length(elements_IDs_blocks{i}),1).*i;
    elements = [elements;elements_blocks{i}];
    elementsMaterialID = [elementsMaterialID;elementMaterialID_block];

end

for i = 1:surf_ct

    boundaryMarker_block = ones(length(facesBoundary_IDs_Block{i}),1).*i;
    facesBoundary = [facesBoundary;facesBoundary_Block{i}];
    boundaryMarker = [boundaryMarker;boundaryMarker_block];

end

% for i = 1:curve_ct
%
% end

MeshStruct.node_IDs = node_ids;
MeshStruct.nodes = nodes;
MeshStruct.facesBoundary = facesBoundary;
MeshStruct.boundaryMarker = boundaryMarker;
MeshStruct.elements = elements;
MeshStruct.elementMaterialID = elementsMaterialID;
MeshStruct.loadNameStruct.MeshName = part_name;
MeshStruct.loadNameStruct.VolumeNames = volumes_names;
MeshStruct.loadNameStruct.SurfaceNames = facesBoundary_Names;

%                 nodes: [415×3 double]
%         facesBoundary: [320×3 double]
%        boundaryMarker: [320×1 double]
%                 faces: [8144×3 double]
%              elements: [2036×4 double]
%     elementMaterialID: [2036×1 double]
%        faceMaterialID: [8144×1 double]
%        loadNameStruct: [1×1 struct]

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function lik = lik_multinom(varargin)
%LIK_MULTINOM    Create a multinom likelihood structure
%
%  Description
%    LIK = LIK_MULTINOM creates multinom likelihood for multi-class
%    count data. The observed numbers in each class with C classes is given
%    as 1xC vector.
%
%    The likelihood is defined as follows:
%                              __ n                __ C
%      p(y|f^1, ..., f^C, z) = || i=1 [ gamma(N+1) || c=1 p_i^c^(y_i^c)/gamma(y_i^c+1)]
%
%    where p_i^c = exp(f_i^c)/ (sum_c=1^C exp(f_i^c)) is the succes
%    probability for class c, which is a function of the latent variable
%    f_i^c for the corresponding class and N=sum(y) is the number of trials.
%
%  See also
%    GP_SET, LIK_*

% Copyright (c) 2010 Jaakko Riihim�ki, Pasi Jyl�nki
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2010 Jarno Vanhatalo

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

  ip=inputParser;
  ip.FunctionName = 'LIK_MULTINOM';
  ip.addOptional('lik', [], @isstruct);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Multinom';
    lik.nondiagW = true;
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Multinom')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_multinom_pak;
    lik.fh.unpak = @lik_multinom_unpak;
    lik.fh.ll = @lik_multinom_ll;
    lik.fh.llg = @lik_multinom_llg;
    lik.fh.llg2 = @lik_multinom_llg2;
    lik.fh.llg3 = @lik_multinom_llg3;
    lik.fh.predy = @lik_multinom_predy;
    lik.fh.invlink = @lik_multinom_invlink;
    lik.fh.recappend = @lik_multinom_recappend;
  end

end

function [w,s] = lik_multinom_pak(lik)
%LIK_MULTINOM_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_MULTINOM_PAK(LIK) takes a likelihood structure LIK and
%    returns an empty verctor W. If Multinom likelihood had
%    parameters this would combine them into a single row vector
%    W (see e.g. lik_negbin). This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%
%  See also
%    LIK_MULTINOM_UNPAK, GP_PAK

  w = []; s = {};
end

function [lik, w] = lik_multinom_unpak(lik, w)
%LIK_MULTINOM_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_MULTINOM_UNPAK(W, LIK) Doesn't do anything.
%
%    If Multinom likelihood had parameters this would extracts them
%    parameters from the vector W to the LIK structure. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%
%  See also
%    LIK_MULTINOM_PAK, GP_UNPAK

  lik=lik;
  w=w;
end

function ll = lik_multinom_ll(lik, y, f, z)
%LIK_MULTINOM_LL  Log likelihood
%
%  Description
%    LL = LIK_MULTINOM_LL(LIK, Y, F) takes a likelihood structure
%    LIK, class counts Y (NxC matrix), and latent values F (NxC
%    matrix). Returns the log likelihood, log p(y|f,z). This
%    subfunction is needed when using Laplace approximation or
%    MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria
%    (DIC, WAIC) computations.
%
%  See also
%    LIK_MULTINOM_LLG, LIK_MULTINOM_LLG3, LIK_MULTINOM_LLG2, GPLA_E

  f=reshape(f,size(y));
  expf = exp(f);
  p = expf ./ repmat(sum(expf,2),1,size(expf,2));
  N = sum(y,2);

  ll = sum(gammaln(N+1) - sum(gammaln(y+1),2) + sum(y.*log(p),2) );

end

function llg = lik_multinom_llg(lik, y, f, param, z)
%LIK_MULTINOM_LLG    Gradient of the log likelihood
%
%  Description
%    LLG = LIK_MULTINOM_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F. Returns
%    the gradient of the log likelihood with respect to PARAM. At
%    the moment PARAM can be 'param' or 'latent'. This subfunction
%    is needed when using Laplace approximation or MCMC for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG2, LIK_MULTINOM_LLG3, GPLA_E

  f=reshape(f,size(y));
  C = size(y,2);
  expf2 = exp(f);
  N=sum(y, 2);
  pi2 = (N*ones(1,C)).*expf2./(sum(expf2, 2)*ones(1,C));
  pi_vec=pi2(:);
  llg = y(:)-pi_vec;

end

function [pi_vec, pi_mat] = lik_multinom_llg2(lik, y, f, param, z)
%LIK_MULTINOM_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_MULTINOM_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F. Returns
%    the Hessian of the log likelihood with respect to PARAM. At
%    the moment PARAM can be only 'latent'. LLG2 is a vector with
%    diagonal elements of the Hessian matrix (off diagonals are
%    zero). This subfunction is needed when using Laplace
%    approximation or EP for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG, LIK_MULTINOM_LLG3, GPLA_E

% multinom:
  [n,nout]=size(y);
  N = sum(y,2)*ones(1,nout);
  f=reshape(f,n,nout);

  expf2 = exp(f);
  pi2 = expf2./(sum(expf2, 2)*ones(1,nout));
  pi_vec=pi2(:).*N(:);

  pi_mat=zeros(nout*n, n);
  for i1=1:nout
    pi_mat((1+(i1-1)*n):(nout*n+1):end)=pi2(:,i1).*sqrt(N(:,i1));
  end
  %     D = diag(pi_vec);
  %     llg2 = -D + pi_mat*pi_mat';

end

function [dw_mat] = lik_multinom_llg3(lik, y, f, param, z)
%LIK_MULTINOM_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_MULTINOM_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F and
%    returns the third gradients of the log likelihood with
%    respect to PARAM. At the moment PARAM can be only 'latent'.
%    LLG3 is a vector with third gradients. This subfunction is
%    needed when using Laplace approximation for inference with
%    non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG, LIK_MULTINOM_LLG2, GPLA_E, GPLA_G

  [n,nout] = size(y);
  f2 = reshape(f,n,nout);

  N=sum(y, 2);
  expf2 = exp(f2);
  pi2 = expf2./(sum(expf2, 2)*ones(1,nout));
  pi_vec=pi2(:);

  dw_mat=zeros(nout,nout,nout,n);

  for cc3=1:nout
    for ii1=1:n

      pic=pi_vec(ii1:n:(nout*n));
      for cc1=1:nout
        for cc2=1:nout

          % multinom third derivatives
          cc_sum_tmp=0;
          if cc1==cc2 && cc1==cc3 && cc2==cc3
            cc_sum_tmp=cc_sum_tmp+pic(cc1);
          end
          if cc1==cc2
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc3);
          end
          if cc2==cc3
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc2);
          end
          if cc1==cc3
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc2);
          end
          cc_sum_tmp=cc_sum_tmp+2*pic(cc1)*pic(cc2)*pic(cc3);

          dw_mat(cc1,cc2,cc3,ii1)=cc_sum_tmp.*N(ii1);
        end
      end
    end
  end

end

function [lpy, Ey, Vary] = lik_multinom_predy(lik, Ef, Varf, yt, zt)
%LIK_MULTINOM_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_MULTINOM_PREDY(LIK, EF, VARF YT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | y) = \int p(yt | f) p(f|y) df.
%    This requires also the incedence counts YT. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%    [LPY, EY, VARY] = LIK_MULTINOM_PREDY(LIK, EF, VARF, YT) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  N=sum(yt,2);
  S=10000;
  [ntest, nout]=size(yt);
  pi=zeros(ntest,nout);
  lpy=zeros(ntest,nout);
  Ey=zeros(ntest,nout);
  Vary=zeros(size(Varf));
  Ef=reshape(Ef(:),ntest,nout);
  [notused,notused,c] =size(Varf);
  if c>1
    mcmc=false;
  else
    mcmc=true;
    Varf=reshape(Varf(:), ntest, nout);
  end
  for i1=1:ntest
    if mcmc
      Sigm_tmp = (Varf(i1,:));
      f_star=bsxfun(@plus, Ef(i1,:), bsxfun(@times, sqrt(Sigm_tmp), ...
        randn(S,nout)));
    else
      Sigm_tmp=(Varf(:,:,i1)'+Varf(:,:,i1))./2;
      f_star=mvnrnd(Ef(i1,:), Sigm_tmp, S);
    end

    tmp = exp(f_star);
    tmp = tmp./(sum(tmp, 2)*ones(1,size(tmp,2)));

    if nargout > 1
        Ey(i1,:) = N(i1).*mean(tmp);
        for z1 = 1:nout;
          for z2 = 1:nout
            for z3=1:S
              Var_tmp(:,:,z3) = (diag(tmp(z3,:)) - tmp(z3,:)'*tmp(z3,:));
            end
            if mcmc
              Vary(i1+(0:nout-1)*ntest,:) = diag(N(i1).*mean(Var_tmp,3));
            else
              Vary(:,:,i1) = N(i1).*mean(Var_tmp,3);
            end
          end
        end
    end
    lpy=[];
    if ~isempty(yt)
      ytmp = repmat(yt(i1,:),S,1);
      lpy(i1,:) = log(mean( mnpdf(ytmp,tmp) ));
    end
  end
  lpy=lpy(:);
  Ey=Ey(:);
end

function p = lik_multinom_invlink(lik, f, z)
%LIK_MULTINOM_INVLINK Returns values of inverse link function
%
%  Description
%    P = LIK_MULTINOM_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_MULTINOM_LL, LIK_MULTINOM_PREDY
p = multinominv(f).*z;
end

function reclik = lik_multinom_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_MULTINOM_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    reclik.type = 'Multinom';
    reclik.nondiagW = true;

    % Set the function handles
    reclik.fh.pak = @lik_multinom_pak;
    reclik.fh.unpak = @lik_multinom_unpak;
    reclik.fh.ll = @lik_multinom_ll;
    reclik.fh.llg = @lik_multinom_llg;
    reclik.fh.llg2 = @lik_multinom_llg2;
    reclik.fh.llg3 = @lik_multinom_llg3;
    reclik.fh.predy = @lik_multinom_predy;
    reclik.fh.invlink = @lik_multinom_invlink;
    reclik.fh.recappend = @lik_multinom_recappend;
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = spm_render(dat,brt,rendfile)
% Render blobs on surface of a 'standard' brain
% FORMAT spm_render(dat,brt,rendfile)
%
% dat      - a struct array of length 1 to 3
%            each element is a structure containing:
%            - XYZ - the x, y & z coordinates of the transformed SPM{.}
%                    values in units of voxels.
%            - t   - the SPM{.} values.
%            - mat - affine matrix mapping from XYZ voxels to MNI.
%            - dim - dimensions of volume from which XYZ is drawn.
% brt      - brightness control:
%            If NaN, then displays using the old style with hot
%            metal for the blobs, and grey for the brain.
%            Otherwise, it is used as a ``gamma correction'' to
%            optionally brighten the blobs up a little.
% rendfile - the file containing the images to render on to (see also
%            spm_surf.m) or a surface mesh file.
%
% Without arguments, spm_render acts as its own UI.
%__________________________________________________________________________
%
% spm_render prompts for details of up to three SPM{.}s that are then
% displayed superimposed on the surface of a 'standard' brain.
%
% The first is shown in red, then green then blue.
%
% The blobs which are displayed are the integral of all transformed t
% values, exponentially decayed according to their depth. Voxels that
% are 10mm behind the surface have half the intensity of ones at the
% surface.
%__________________________________________________________________________
% Copyright (C) 1996-2019 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_render.m 7577 2019-04-24 08:59:56Z guillaume $

SVNrev = '$Rev: 7577 $';

global prevrend
if ~isstruct(prevrend)
    prevrend = struct('rendfile','', 'brt',[], 'col',[]);
end

varargout = {};

%-Parse arguments, get data if not passed as parameters
%==========================================================================

%-Get surface
%--------------------------------------------------------------------------
if nargin < 3 || isempty(prevrend.rendfile)
    spm('FnBanner',mfilename,SVNrev);
    [rendfile, sts] = spm_select(1,{'mat','mesh'},'Render file');
    if ~sts, return; end
end
prevrend.rendfile = rendfile;

ext = spm_file(rendfile,'ext');
loadgifti = false;
if strcmpi(ext,'mat')
    load(rendfile);
    if ~exist('rend','var') && ~exist('Matrixes','var')
        loadgifti = true;
    end
end
if ~strcmpi(ext,'mat') || loadgifti
    try
        rend = gifti(rendfile);
    catch
        error('Cannot read  render file "%s".\n', rendfile);
    end
    if ~isfield(rend,'vertices')
        try
            M = rend;
            rend = gifti(rend.private.metadata(1).value);
            try, rend.cdata = full(M.cdata); end
        catch
            error('Cannot find a surface mesh to be displayed.');
        end
    end
    rend = export(rend,'patch');
end

%-Get data
%--------------------------------------------------------------------------
if isfield(rend,'facevertexcdata')
    dat = rend.facevertexcdata;
    num = 1;
elseif nargin < 1
    spm('FigName','Results: render');

    num   = spm_input('Number of sets',1,'1 set|2 sets|3 sets',[1 2 3],1);

    for i = 1:num
        [SPM,xSPM] = spm_getSPM;
        dat(i) = struct('XYZ', xSPM.XYZ,...
                        't',   xSPM.Z',...
                        'mat', xSPM.M,...
                        'dim', xSPM.DIM);
    end
    showbar = 1;
else
    if isstruct(dat)
        num     = length(dat);
        showbar = 0;
    else
        files = cellstr(dat);
        clear dat
        num = numel(files);
        for i=1:num
            V  = spm_vol(files{i});
            d  = spm_read_vols(V);
            id = find(isfinite(d));
            [X,Y,Z]    = ind2sub(V.dim,id);
            dat(i).XYZ = [X Y Z]';
            dat(i).t   = d(id);
            dat(i).mat = V.mat;
            dat(i).dim = V.dim';
            clear X Y Z id
        end
        showbar = 1;
    end
end

%-Get brightness & colours (mesh)
%--------------------------------------------------------------------------
if isfield(rend,'vertices')
    if num == 1
        col = hot(256);
    else
        col = eye(3);
        if spm_input('Which colours?','+1','b',{'RGB','Custom'},[0 1],1)
            for k = 1:num
                col(k,:) = uisetcolor(col(k,:),sprintf('Colour of blob set %d',k));
            end
        end
    end
    surf_rend(dat,rend,col);
    return
end

%-Get brightness & colours (image)
%--------------------------------------------------------------------------
if nargin < 2  || isempty(prevrend.brt)
    brt = 1;
    if num==1
        brt = spm_input('Style',1,'new|old',[1 NaN], 1);
    end

    if isfinite(brt)
        brt = spm_input('Brighten blobs','+1','none|slightly|more|lots',[1 0.75 0.5 0.25], 1);
        col = eye(3);
        % ask for custom colours & get rgb values
        %------------------------------------------------------------------
        if spm_input('Which colours?','+1','b',{'RGB','Custom'},[0 1],1)
            for k = 1:num
                col(k,:) = uisetcolor(col(k,:),sprintf('Colour of blob set %d',k));
            end
        end
    else
        col = [];
    end
elseif isfinite(brt) && isempty(prevrend.col)
    col = eye(3);
elseif isfinite(brt)  % don't need to check prevrend.col again
    col = prevrend.col;
else
    col = [];
end
prevrend.brt = brt;
prevrend.col = col;

%-Perform the rendering
%==========================================================================
spm('Pointer','Watch');

if showbar, spm_progress_bar('Init', size(dat,1)*length(rend),...
            'Formatting Renderings', 'Number completed'); end
for i=1:length(rend)
    rend{i}.max=0;
    rend{i}.data = cell(size(dat,1),1);
    if issparse(rend{i}.ren)
        % Assume that images have been DCT compressed
        % - the SPM99 distribution was originally too big.
        d = size(rend{i}.ren);
        B1 = spm_dctmtx(d(1),d(1));
        B2 = spm_dctmtx(d(2),d(2));
        rend{i}.ren = B1*rend{i}.ren*B2';
        % the depths did not compress so well with
        % a straight DCT - therefore it was modified slightly
        rend{i}.dep = exp(B1*rend{i}.dep*B2')-1;
    end
    rend{i}.ren(rend{i}.ren>=1) = 1;
    rend{i}.ren(rend{i}.ren<=0) = 0;
    if showbar, spm_progress_bar('Set', i); end
end
if showbar, spm_progress_bar('Clear'); end

if showbar, spm_progress_bar('Init', length(dat)*length(rend),...
            'Making pictures', 'Number completed'); end

mx = zeros(length(rend),1)+eps;
mn = zeros(length(rend),1);

for j=1:length(dat)
    XYZ = dat(j).XYZ;
    t   = dat(j).t;
    dim = dat(j).dim;
    mat = dat(j).mat;

    for i=1:length(rend)

        % transform from Talairach space to space of the rendered image
        %------------------------------------------------------------------
        M1  = rend{i}.M*mat;
        zm  = sum(M1(1:2,1:3).^2,2).^(-1/2);
        M2  = diag([zm' 1 1]);
        M  = M2*M1;
        cor = [1 1 1 ; dim(1) 1 1 ; 1 dim(2) 1; dim(1) dim(2) 1 ;
               1 1 dim(3) ; dim(1) 1 dim(3) ; 1 dim(2) dim(3); dim(1) dim(2) dim(3)]';
        tcor= M(1:3,1:3)*cor + M(1:3,4)*ones(1,8);
        off = min(tcor(1:2,:)');
        M2  = spm_matrix(-off+1)*M2;
        M  = M2*M1;
        xyz = (M(1:3,1:3)*XYZ + M(1:3,4)*ones(1,size(XYZ,2)));
        d2  = ceil(max(xyz(1:2,:)'));

        % Calculate 'depth' of values
        %------------------------------------------------------------------
        if ~isempty(d2)
            dep = spm_slice_vol(rend{i}.dep,spm_matrix([0 0 1])*inv(M2),d2,1);
            z1  = dep(round(xyz(1,:))+round(xyz(2,:)-1)*size(dep,1));

            if ~isfinite(brt), msk = find(xyz(3,:) < (z1+20) & xyz(3,:) > (z1-5));
            else msk = find(xyz(3,:) < (z1+60) & xyz(3,:) > (z1-5)); end
        else
            msk = [];
        end

        if ~isempty(msk)

            % Generate an image of the integral of the blob values.
            %--------------------------------------------------------------
            xyz = xyz(:,msk);
            if ~isfinite(brt), t0  = t(msk);
            else
                dst = xyz(3,:) - z1(msk);
                dst = max(dst,0);
                t0  = t(msk).*exp((log(0.5)/10)*dst)';
            end
            X0  = full(sparse(round(xyz(1,:)), round(xyz(2,:)), t0, d2(1), d2(2)));
            hld = 1; if ~isfinite(brt), hld = 0; end
            X   = spm_slice_vol(X0,spm_matrix([0 0 1])*M2,size(rend{i}.dep),hld);
            msk = find(X<0);
            X(msk) = 0;
        else
            X = zeros(size(rend{i}.dep));
        end

        % Brighten the blobs
        %------------------------------------------------------------------
        if isfinite(brt), X = X.^brt; end

        mx(j) = max([mx(j) max(max(X))]);
        mn(j) = min([mn(j) min(min(X))]);

        rend{i}.data{j} = X;

        if showbar, spm_progress_bar('Set', i+(j-1)*length(rend)); end
    end
end

mxmx = max(mx);
mnmn = min(mn);

if showbar, spm_progress_bar('Clear'); end

%-Display
%==========================================================================
if ~spm('CmdLine')
    Fgraph = spm_figure('GetWin','Graphics');
    spm_results_ui('Clear',Fgraph);

    nrow = ceil(length(rend)/2);
    if showbar, hght = 0.95; else hght = 0.5; end
    % subplot('Position',[0, 0, 1, hght]);
    ax = axes('Parent',Fgraph,...
        'units','normalized',...
        'Position',[0, 0, 1, hght],...
        'Visible','off');
    image(0,'Parent',ax);
    set(ax,'YTick',[],'XTick',[]);
end

rgb = cell(1,length(rend));

if ~isfinite(brt)
    % Old style split colourmap display
    %----------------------------------------------------------------------
    split = [gray(64); hot(64)];
    if ~spm('CmdLine'), colormap(split); end
    for i=1:length(rend)
        ren = rend{i}.ren;
        X   = (rend{i}.data{1}-mnmn)/(mxmx-mnmn);
        msk = find(X);
        ren(msk) = X(msk)+(1+1.51/64);
        rgb{i} = ind2rgb(round(ren*64),split);
        if ~spm('CmdLine')
            ax = axes('Parent',Fgraph,...
                'units','normalized',...
                'Position',[rem(i-1,2)*0.5, floor((i-1)/2)*hght/nrow, 0.5, hght/nrow],...
                'Visible','off');
            image(ren*64,'Parent',ax);
            set(ax,'DataAspectRatio',[1 1 1],...
                'PlotBoxAspectRatioMode','auto',...
                'YTick',[],'XTick',[],...
                'XDir','normal','YDir','normal');
        end
    end
else
    % Combine the brain surface renderings with the blobs, and display
    % using 24 bit colour
    %----------------------------------------------------------------------
    for i=1:length(rend)
        ren = rend{i}.ren;
        X = cell(3,1);
        for j=1:length(rend{i}.data)
            X{j} = rend{i}.data{j}/(mxmx-mnmn)-mnmn;
        end
        for j=(length(rend{i}.data)+1):3
            X{j}=zeros(size(X{1}));
        end

        rgb{i} = zeros([size(ren) 3]);
        tmp = ren.*max(1-X{1}-X{2}-X{3},0);
        for k = 1:3
            rgb{i}(:,:,k) = tmp + X{1}*col(1,k) + X{2}*col(2,k) +X{3}*col(3,k);
        end
        rgb{i}(rgb{i}>1) = 1;
        if ~spm('CmdLine')
            ax = axes('Parent',Fgraph,...
                'units','normalized',...
                'Position',[rem(i-1,2)*0.5, floor((i-1)/2)*hght/nrow, 0.5, hght/nrow],...
                'nextplot','add', ...
                'Visible','off');
            image(rgb{i},'Parent',ax);
            set(ax,'DataAspectRatio',[1 1 1],...
                'PlotBoxAspectRatioMode','auto',...
                'YTick',[],'XTick',[],...
                'XDir','normal','YDir','normal');
        end
    end
end

spm('Pointer','Arrow');

if nargout
    for i=1:numel(rgb)
        for ch=1:size(rgb{i},3)
            rgb{i}(:,:,ch) = flipud(rgb{i}(:,:,ch));
        end
    end
    varargout = { rgb };
end

%==========================================================================
% function surf_rend(dat,rend,col)
%==========================================================================
function surf_rend(dat,rend,col)

%-Setup figure and axis
%--------------------------------------------------------------------------
Fgraph = spm_figure('GetWin','Graphics');
spm_results_ui('Clear',Fgraph);

ax0 = axes(...
    'Tag',      'SPMMeshRenderBackground',...
    'Parent',   Fgraph,...
    'Units',    'normalized',...
    'Color',    [1 1 1],...
    'XTick',    [],...
    'YTick',    [],...
    'Position', [-0.05, -0.05, 1.05, 0.555]);

ax = axes(...
    'Parent',   Fgraph,...
    'Units',    'normalized',...
    'Position', [0.05, 0.05, 0.9, 0.4],...
    'Visible',  'off');

H = spm_mesh_render('Disp',rend,struct('parent',ax));
spm_mesh_render('Overlay',H,dat,col);
camlight(H.light);

try
    setAllowAxesRotate(H.rotate3d, setxor(findobj(Fgraph,'Type','axes'),ax), false);
end

%-Register with MIP
%--------------------------------------------------------------------------
try % meaningless when called outside spm_results_ui
    hReg = spm_XYZreg('FindReg',spm_figure('GetWin','Interactive'));
    xyz  = spm_XYZreg('GetCoords',hReg);
    hs   = mydispcursor('Create',ax,dat.mat,xyz);
    spm_XYZreg('Add2Reg',hReg,hs,@mydispcursor);
end

%==========================================================================
function varargout = mydispcursor(varargin)

switch lower(varargin{1})
    %======================================================================
    case 'create'
    %======================================================================
    % hMe = mydispcursor('Create',ax,M,xyz)
    ax  = varargin{2};
    M   = varargin{3};
    xyz = varargin{4};

    [X,Y,Z] = sphere;
    vx = sqrt(sum(M(1:3,1:3).^2));
    X = X*vx(1) + xyz(1);
    Y = Y*vx(2) + xyz(2);
    Z = Z*vx(3) + xyz(3);
    hold(ax,'on');
    hs = surf(X,Y,Z,'parent',ax,...
        'EdgeColor','none','FaceColor',[1 0 0],'FaceLighting', 'gouraud');
    set(hs,'UserData',xyz);

    varargout = {hs};

    %======================================================================
    case 'setcoords'    % Set co-ordinates
    %======================================================================
    % [xyz,d] = mydispcursor('SetCoords',xyz,hMe,hC)
    hMe  = varargin{3};
    pxyz = get(hMe,'UserData');
    xyz  = varargin{2};

    set(hMe,'XData',get(hMe,'XData') - pxyz(1) + xyz(1));
    set(hMe,'YData',get(hMe,'YData') - pxyz(2) + xyz(2));
    set(hMe,'ZData',get(hMe,'ZData') - pxyz(3) + xyz(3));
    set(hMe,'UserData',xyz);

    varargout = {xyz,[]};

    %======================================================================
    otherwise
    %======================================================================
    error('Unknown action string')

end

author: spm


text:
function ret = filterBlock(in)
    in = in + 1e-8;
%     idx = [2 4 5 6 8];
%     in = in(idx, :);
    mid = in(5, :);
%     ret = sum(abs(bsxfun(@minus, log(mid), log(in)))) / 8;
    temp = abs(bsxfun(@minus, log(mid), log(in)));
    ret = sum((temp)) / 8;
end

author: andrewssobral

repo: lrslibrary


text:
function ethogram_export(hFig0)

if ~exist('hFig0','var')
  hFig0 = gcf;
end

v = version('-release');
if ~strcmp(v,'2014b')
  error('ethogram_export:version','Export is intended for use with MATLAB R2014b.');
end

% copy the entire ethogram
hFig = copyobj(hFig0,0);

% delete stuff
TAGS2DELETE = {'udUserSelectionPatch' 'udCustomGroupMarker' 'udGroupStatsAx'};
for i = 1:numel(TAGS2DELETE)
  h = findall(hFig,'userdata',TAGS2DELETE{i});
  delete(h);
end

% massage x-axis, legend
XTICK = 500;
ax = findall(hFig,'type','axes');
assert(isscalar(ax),'Multiple axes found.');
set(ax,'XTick',XTICK);
leg = findall(hFig,'type','legend');
assert(isscalar(leg),'Multiple legends found.');
set(leg,'Location','NorthOutside')

% get a location/file
exportpath = ExpPP.loadConfigVal('exportpath');
if isempty(exportpath)
  exportpath = pwd;
end
filterspec = fullfile(exportpath,'*.fig');
[fname,pname] = uiputfile(filterspec,'Export data');
if isequal(fname,0) || isequal(pname,0)
  return;
end

[~,fbase] = fileparts(fname);
figname = fullfile(pname,fname);
pdfname = fullfile(pname,[fbase '.pdf']);
epsname = fullfile(pname,[fbase '.eps']);

% save fig, print pdf, eps
hgsave(hFig,figname);
hlpExport(hFig,pdfname,'pdf');
hlpExport(hFig,epsname,'epsc2');

ExpPP.saveConfigVal('exportpath',pname);

delete(hFig);

function hlpExport(hFig,fullfname,ext)
if exist(fullfname,'file')
  [fname,pname] = uiputfile(fullfname,['Export ' ext]);
  if isequal(fname,0) || isequal(pname,0)
    fullfname = [];
  else
    fullfname = fullfile(pname,fname);
  end
end
if ~isempty(fullfname)
  print(hFig,['-d' ext],'-noui',fullfname);
end

author: kristinbranson

repo: JAABA


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR:  SSD versus discretization width (midpoint quadrature rule)
%
% - load data (setup2DHNSPData)
% - interpolate on various grids (splineInter)
% - compute SSD and plot it
%==============================================================================

clear; help(mfilename);

setup2DHNSPData; clf; h = []; Q = [];
imgModel('reset','imgModel','splineInter');
[T,R] = imgModel('coefficients',dataT,dataR,omega,'out',0);
for j=1:10,
  m    = 2^j*[1,1];
  h(j) = prod((omega(2:2:end)-omega(1:2:end))./m);
  xc   = getCellCenteredGrid(omega,m);
  res  = imgModel(T,omega,xc) - imgModel(R,omega,xc);
  psi  = 0.5*h(j)*res'*res;
  Q(j) = psi;
end;
figure(1); clf; p1=semilogx(h/h(1),Q+eps,'kx',h/h(1),Q,'k-');
%==============================================================================


text:
classdef prtClassCorrelation < prtClass
    % prtClassCorrelation  Correlation classifier
    %
    %    CLASSIFIER = prtClassCorrelation returns a correlation classifier
    %
    %    CLASSIFIER = prtClassCorrelation(PROPERTY1, VALUE1, ...)
    %    constructs a prtClassCorrelation object CLASSIFIER with properties
    %    as specified by PROPERTY/VALUE pairs.
    %
    %    A prtClassCorrelation object inherits all properties from the
    %    abstract class prtClass.
    %
    %    A prtClassCorrelation object inherits the TRAIN, RUN,
    %    CROSSVALIDATE and KFOLDS methods from prtAction. It also inherits
    %    the PLOT method from prtClass.
    %
    %    Example:
    %
    %   TestDataSet = prtDataGenUnimodal;      % Create some test and
    %   TrainingDataSet = prtDataGenUnimodal;  % training data
    %   classifier = prtClassCorrelation;           % Create a classifier
    %   classifier = classifier.train(TrainingDataSet);    % Train
    %   classified = run(classifier, TestDataSet);         % Test
    %   subplot(2,1,1);
    %   classifier.plot;
    %   subplot(2,1,2);
    %   [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %   h = plot(pf,pd,'linewidth',3);
    %   title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %    See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
    %    prtClassMap, prtClassCap, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
    %    prtClassKnn, prtClassFld, prtClassRvm, prtClassGlrt,  prtClass

    properties (SetAccess=private)
        name = 'Correlation Classifier'
        nameAbbreviation = 'CORR'
        isNativeMary = true;
    end

    properties
        useMeanLibrary = false;
    end
    properties (SetAccess=protected)
        xMeans

    end

    methods
        function Obj = prtClassCorrelation(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)
        function Obj = preTrainProcessing(Obj,DataSet)
            if ~Obj.verboseStorage
                warning('prtClassKnn:verboseStorage:false','prtClassKnn requires verboseStorage to be true; overriding manual settings');
            end
            Obj.verboseStorage = true;
            Obj = preTrainProcessing@prtClass(Obj,DataSet);
        end

        function Obj = trainAction(Obj,DataSet)

            Obj.xMeans = mean(DataSet.getObservations);

        end

        function DataSet = runAction(Obj,DataSet)

            X = DataSet.X;
            X = bsxfun(@minus,X,Obj.xMeans);
            isOneDimension = (size(X,2) == 1);

            Yout = zeros(size(X,1),Obj.dataSetSummary.nClasses);

            for iClass = 1:Obj.dataSetSummary.nClasses
                if Obj.useMeanLibrary
                    cTrainingX = mean(Obj.dataSet.getObservationsByClassInd(iClass),1);
                else
                    cTrainingX = Obj.dataSet.getObservationsByClassInd(iClass);
                end

                cTrainingX = bsxfun(@minus, cTrainingX, Obj.xMeans);

                if isOneDimension
                     cCorrMat = cTrainingX *X';
                else
                     cCorrMat = cTrainingX *X'./((sum(cTrainingX .^2,2).^.5)*(sum(X.^2,2).^.5)');
                end

                Yout(:,iClass) = max(cCorrMat,[],1);
             end

            Yout = (Yout+1)/2; % Make things be between 0 and 1;
            Yout(Yout>1) = 1;
            Yout(Yout<0) = 0;

            DataSet = DataSet.setObservations(Yout);
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% [descriptors, locs] = sift(img)
%
% This function returns IMAGE's SIFT keypoints.
%   Input parameters:
%     img: the image.
%
%   Returned:
%     descriptors: a K-by-128 matrix, where each row gives an invariant
%         descriptor for one of the K keypoints.  The descriptor is a vector
%         of 128 values normalized to unit length.
%     locs: K-by-4 matrix, in which each row has the 4 values for a
%         keypoint location (row, column, scale, orientation).  The
%         orientation is in the range [-PI, PI] radians.
%
% Credits: Thanks for initial version of this program to D. Alvaro and
%          J.J. Guerrero, Universidad de Zaragoza (modified by D. Lowe)

function [descriptors, locs] = sift(img)

% If you have the Image Processing Toolbox, you can uncomment the following
%   lines to allow input of color images, which will be converted to grayscale.
if isrgb(img)
   img = rgb2gray(img);
end

[rows, cols] = size(img);

% Convert into PGM imagefile, readable by "keypoints" executable
f = fopen('tmp.pgm', 'w');
if f == -1
    error('Could not create file tmp.pgm.');
end
fprintf(f, 'P5\n%d\n%d\n255\n', cols, rows);
fwrite(f, img', 'uint8');
fclose(f);

% Call keypoints executable
if isunix
    command = '!./sift ';
else
    command = '!siftWin32 ';
end
command = [command ' <tmp.pgm >tmp.key'];
eval(command);

% Open tmp.key and check its header
g = fopen('tmp.key', 'r');
if g == -1
    error('Could not open file tmp.key.');
end
[header, count] = fscanf(g, '%d %d', [1 2]);
if count ~= 2
    error('Invalid keypoint file beginning.');
end
num = header(1);
len = header(2);
if len ~= 128
    error('Keypoint descriptor length invalid (should be 128).');
end

% Creates the two output matrices (use known size for efficiency)
locs = double(zeros(num, 4));
descriptors = double(zeros(num, 128));

% Parse tmp.key
for i = 1:num
    [vector, count] = fscanf(g, '%f %f %f %f', [1 4]); %row col scale ori
    if count ~= 4
        error('Invalid keypoint file format');
    end
    locs(i, :) = vector(1, :);

    [descrip, count] = fscanf(g, '%d', [1 len]);
    if (count ~= 128)
        error('Invalid keypoint file value.');
    end
    % Normalize each input vector to unit length
    descrip = descrip / sqrt(sum(descrip.^2));
    descriptors(i, :) = descrip(1, :);
end
fclose(g);

delete('tmp.pgm');

author: Sable

repo: mcbench-benchmarks


text:
%%
example = qMRusage(qMT_bSSFP,'fit');

% Duplicate %
prctpos=strfind(example,'%');
for ia = 1:length(prctpos)
    example=[example(1:((prctpos(ia)+ia-1)-1)) '%' example((prctpos(ia)+ia-1):end)];
end

fid = fopen('bSSFP_batch.m','w');
fprintf(fid,example);
fclose(fid);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
addpath(genpath('../mutils/My/'));
addpath(genpath('../ptv'));

%%
basepth = '../../data_prj/dir_dataset/DIR_files/';

TRE = zeros(10,1);
TIME = zeros(10,1);
TREstd = zeros(10,1);
TREc = zeros(10, 3);
TREcstd = zeros(10, 3);
opts = [];
opts.grid_spacing = [4, 4, 3];  % grid spacing in pixels
opts.display = 'off';
opts.k_down = 0.7;
opts.interp_type = 0;
opts.metric = 'nuclear';
opts.max_iters =  10;
opts.max_iters = 110;
opts.check_gradients = 100*0;
opts.interp_type = 0;

opts.spline_order = 1;

% opts.singular_coefs = sqrt(1:6);
% opts.singular_coefs(1) = 0;
opts.spat_reg_p_val = 0.75;
opts.jac_reg = 10;

opts.D1Lp = 1e-5;
opts.D1Lp = 8e-6;
% opts.D1Lp = 3e-6*0;

% opts.D1Lp = 1e-4; opts.opt_method = 'adam'; % not distributed
opts.opt_method = 'lbfgs';
% opts.opt_method = 'adam';

opts.nlvl = [];

for idx = 8 %1:10
    pts_struct = DIR_get_all_points_for_the_case(idx, basepth);
    % read images
    [volmov, spc] = read_DIR_volume_4dCT(idx, 5, basepth);
    pts_mov = pts_struct.extreme.e;

    init_size = size(volmov);
    min_max1 = [ min(pts_struct.extreme.b, [], 1)', max(pts_struct.extreme.b, [], 1)'];
    min_max2 = [ min(pts_struct.extreme.e, [], 1)', max(pts_struct.extreme.e, [], 1)'];
    min_max = [ min(min_max1(:, 1), min_max2(:, 1)), max(min_max1(:, 2), min_max2(:, 2))];
    d = [10, 10, 5];
    crop_v = [ max(1, min_max(1,1) - d(1)), min(size(volmov, 1), min_max(1,2) + d(1)); ...
               max(1, min_max(2,1) - d(2)), min(size(volmov, 2), min_max(2,2) + d(2)); ...
               max(1, min_max(3,1) - d(3)), min(size(volmov, 3), min_max(3,2) + d(3));];

    spc_orig = spc;
    for i = 1 : 6
        [tmp, spc] = read_DIR_volume_4dCT(idx, i-1, basepth);
        tmp = crop_data(tmp, crop_v);
        bszv = size(tmp);
        spc_tmp = [1, 1, 1];
        tmpu = volresize(tmp, round(bszv .* spc .* spc_tmp), 1);
        spc = [1,1,1] ./ spc_tmp;

%         tmpu = tmp;
        if i == 1
            vols = zeros([size(tmpu), 6]);
        end
        vols(:,:,:, i) = tmpu;
    end
    pts_fix = pts_struct.extreme.b;
    vols = img_thr(vols, 50, 1100, 1);
    opts.pix_resolution = spc;
    opts.border_mask = 5;
    display('run')

    vols = reshape(vols, [size(vols,1),size(vols,2),size(vols,3), 1, size(vols, 4)]);
    tic
    [voldef, Tmin3_out, Kmin, itinfo] = ptv_register(vols, [], opts);
    toc
end

%
id1 = 1;
id2 = 6;
[Tnew, imgs_1] = remap_displacements( Tmin3_out, id1, vols, opts.pix_resolution);

Tptv_rsz = cat(4, volresize(Tnew(:,:,:,1,end), bszv), volresize(Tnew(:,:,:,2,end), bszv), volresize(Tnew(:,:,:,3,end), bszv));
voldef_rsz = volresize(voldef(:,:,:, 1), bszv);

[~, Tr] = uncrop_data(voldef_rsz(:,:,:,1), Tptv_rsz(:,:,:,:), crop_v, init_size);
[pt_errs_phys, pts_moved_pix, TRE_phys, TREstd_phys] = DIR_movepoints(pts_mov, pts_fix, Tr, spc_orig, [0,0,0]);
TRE(idx) = TRE_phys;
TREstd(idx) = TREstd_phys;
fprintf('$$$$$$$$ For idx=%d error = %.3f  std = %.3f,  time=%.1f\n', idx, TRE(idx), TREstd(idx), TIME(idx));

%%
imgs_exp = [];
Tpix = conv_3d_T_from_phys_to_pix(Tmin3_out(:,:,:,:, :), opts.pix_resolution);
for i = 1 : size(Tpix, 5)
    hold off;
    imagesc(squeeze(vols(:, round(end/2), :, :, i))');
    colormap gray;
    daspect([1,1,1]);
    hold on;
    plot_wiremesh(permute(squeeze(Tpix(:, round(end/2), :, [3,2], i)), [2,1,3]), 6, 'r-')

    F = getframe;
    imgs_exp = cat(4, imgs_exp, double(F.cdata));

    pause(0.1);
end
%%
implay_2d_deform( permute(squeeze(vols(:, round(end/2), :, :, :)), [2, 1, 3]),  permute(squeeze(imgs_1(:, round(end/2), :, :, :)), [2, 1, 3]), ...
    permute(squeeze(Tpix(:, round(end/2), :, :, :)), [2, 1, 3, 4]), 5, false, [0,1]);
%%
ima = permute(squeeze(voldef(:, round(end/2), :, :, :)), [2,1,3]);
savegif('results/DIR_group_track.gif', squeeze([imgs_exp]/250), 1/5);
savegif('results/DIR_group_stab.gif', ima, 1/5);

repo: pTVreg


text:
function [dim] = pos2dim(pos)

% POS2DIM reconstructs the volumetric dimensions from an ordered list of
% positions.
%
% Use as
%   [dim] = pos2dim(pos)
% where pos is an ordered list of positions.
%
% The output dim is a 3-element vector which correspond to the 3D
% volumetric dimensions
%
% See also POS2TRANSFORM

% Copyright (C) 2009, Jan-Mathijs Schoffelen

if isstruct(pos)
  % the input is a FieldTrip data structure
  pos = pos.pos;
end

% this part depends on the assumption that the list of positions is describing a full 3D volume in
% an ordered way which allows for the extraction of a transformation matrix, i.e. slice by slice
npos = size(pos, 1);
dpos = zscore(abs(diff(pos, [], 1)));

dim        = nan(1, 3);
[tmp, ind] = max(dpos, [], 2);
dim(1)     = find(tmp>1.5, 1, 'first');
dpos       = dpos(dim:dim:npos-1, :);
[tmp, ind] = max(dpos(:, setdiff(1:3, ind(dim(1)))), [], 2);
dim(2)     = find(tmp>1.1*min(tmp), 1, 'first'); % this threshold seems to work on what I tried out
dim(3)     = npos./prod(dim(1:2));


text:
function load_existing_bgrid_version_A()
    % load existing b-value grid
    % VERSION A
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();
    % extracted from the various b cross routines

    [file1,path1] = uigetfile('*.mat','b-value gridfile');
    if length(path1) > 1

        load([path1 file1])
        xsecx = newa(:,end)';
        xsecy = newa(:,7);
        xvect = gx; yvect = gy;
        tmpgri=zeros((length(xvect)*length(yvect)),2);

        normlap2=nan(length(tmpgri(:,1)),1)
        normlap2(ll)= bvg(:,1);
        valueMap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,5);
        r=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,6);
        meg=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,2);
        old1 =reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,7);
        pro=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,8);
        avm=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,9);
        stanm=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,10);
        maxm=reshape(normlap2,length(yvect),length(xvect));

        old = valueMap;

        view_bv2
    else
        return
    end
end

author: CelsoReyes


text:
function [data] = ft_denoise_ssp(cfg, data)

% FT_DENOISE_SSP projects out topographies based on ambient noise on
% Neuromag/Elekta/MEGIN systems. These topographies are estimated during maintenance
% visits from the engineers of MEGIN
%
% Use as
%   [data] = ft_denoise_ssp(cfg, data)
% where data should come from FT_PREPROCESSING and the configuration
% should contain
%   cfg.ssp        = 'all' or a cell array of SSP names to apply (default = 'all')
%   cfg.trials     = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.updatesens = 'no' or 'yes' (default = 'yes')
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_PREPROCESSING, FT_DENOISE_SYNTHETIC, FT_DENOISE_PCA

% Copyright (C) 2004-2022, Gianpaolo Demarchi, Lau
% Møller Andersen, Robert Oostenveld, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'required',   {'ssp'});

% set the defaults
cfg.ssp        = ft_getopt(cfg, 'ssp', 'all');
cfg.trials     = ft_getopt(cfg, 'trials', 'all', 1);
cfg.updatesens = ft_getopt(cfg, 'updatesens', 'yes');

% store the original type of the input data
dtype = ft_datatype(data);

% check if the input data is valid for this function
% this will convert timelocked input data to a raw data representation if needed
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes', 'hassampleinfo', 'yes');

% check whether it is neuromag data
if ~ft_senstype(data, 'neuromag')
  ft_error('SSP vectors can only be applied to neuromag data');
end

% select trials of interest
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
data   = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% remember the original channel ordering
labelold = data.label;

% apply the balancing to the MEG data and to the gradiometer definition
current = data.grad.balance.current;
desired = cfg.ssp;

if ~strcmp(current, 'none')
  % first undo/invert the previously applied balancing
  try
    current_montage = data.grad.balance.(current);
  catch
    ft_error('unknown balancing for input data');
  end
  fprintf('converting the data from "%s" to "none"\n', current);
  data = ft_apply_montage(data, current_montage, 'keepunused', 'yes', 'inverse', 'yes');
  if istrue(cfg.updatesens)
    fprintf('converting the sensor description from "%s" to "none"\n', current);
    data.grad = ft_apply_montage(data.grad, current_montage, 'keepunused', 'yes', 'inverse', 'yes');
    data.grad.balance.current = 'none';
  end
end % if current

if ~strcmp(desired, 'none')
  % then apply the desired balancing
  if strcmp(desired, 'all')
    desireds = fieldnames(data.grad.balance);
  else
      desireds = cfg.ssp;
    if ~iscell(desireds)
        ft_error('cfg.ssp must be a cell array of projector names')
    end

  end
  for desired_index = 1:length(desireds)
    desired = desireds{desired_index};
    if ~strcmp(desired, 'current')
      try
        desired_montage = data.grad.balance.(desired);
      catch
        ft_error('unknown balancing for input data');
      end
      fprintf('converting the data from "none" to "%s"\n', desired);
      data = ft_apply_montage(data, desired_montage, 'keepunused', 'yes', 'balancename', desired);
      if istrue(cfg.updatesens)
        fprintf('converting the sensor description from "none" to "%s"\n', desired);
        data.grad = ft_apply_montage(data.grad, desired_montage, 'keepunused', 'yes', 'balancename', desired);
      end
    end % if desired
  end
end

% reorder the channels to stay close to the original ordering
[selold, selnew] = match_str(labelold, data.label);
if numel(selnew)==numel(labelold)
  for i=1:numel(data.trial)
    data.trial{i} = data.trial{i}(selnew,:);
  end
  data.label = data.label(selnew);
else
  ft_warning('channel ordering might have changed');
end

% convert back to input type if necessary
switch dtype
  case 'timelock'
    data = ft_checkdata(data, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug

ft_postamble previous   data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data


text:
%% Sharp Color Keys
%
%%
% In this section we discuss color keys that are particular useful when
% analyzing data with very small deviation in orientation. Let us consider
% the following calcite data set

% plotting conventions
plotx2east, plotb2east
mtexdata sharp

ebsd = ebsd('calcite');

ipfKey = ipfColorKey(ebsd);

close all;
plot(ebsd,ipfKey.orientation2color(ebsd.orientations))

%%
% and have a look into the 100 inverse pole figure.

% compute the positions in the inverse pole figure
h = ebsd.orientations .\ vector3d.X;
h = project2FundamentalRegion(h);

% compute the azimuth angle in degree
color = h.rho ./ degree;

plotIPDF(ebsd.orientations,vector3d.X,'property',color,'MarkerSize',3,'grid')
mtexColorbar

%%
% We see that all individual orientations are clustered around azimuth
% angle 115 degrees with some outliers at 65 degree. In order to
% increase the contrast for the main group, we restrict the color range from
% 110 degree to 120 degree.

caxis([110 120]);

% by the following lines we colorcode the outliers in purple.
cmap = colormap;
cmap(end,:) = [1 0 1]; % make last color purple
cmap(1,:) = [1 0 1];   % make first color purple
colormap(cmap)

%%
% The same colorcoding we can now apply to the EBSD map.

% plot the data with the customized color
plot(ebsd,color)

% set scaling of the angles to 110 - 120 degree
caxis([110 120]);

% colorize outliers in purple.
cmap = colormap;
cmap(end,:) = [1 0 1];
cmap(1,:) = [1 0 1];
colormap(cmap)

%% Sharpening the default colorcoding
% Next, we want to apply the same ideas as above to the default MTEX color
% key, i.e. we want to stretch the colors such that they cover just the
% orientations of interest.

ipfKey = ipfHSVKey(ebsd.CS.properGroup);

% To this end, we first compute the inverse pole figure direction such that
% the mean orientation is just at the gray spot of the inverse pole figure
ipfKey.inversePoleFigureDirection = mean(ebsd.orientations,'robust') * ipfKey.whiteCenter;

close all;
plot(ebsd,ipfKey.orientation2color(ebsd.orientations))

%%
% We observe that the orientation map is almost completely gray, except for
% the  outliers which appears black. Next, we use the option |maxAngle| to
% increase contrast in the grayish part

ipfKey.maxAngle = 7.5*degree;
plot(ebsd,ipfKey.orientation2color(ebsd.orientations))

%%
% You may play around with the option |maxAngle| to obtain better
% results. As for interpretation keep in mind that white color represents
% the mean orientation and the color becomes more saturated and later dark
% as the orientation to color diverges from the mean orientation.
%
% Let's have a look at the corresponding color map.

plot(ipfKey,'resolution',0.25*degree)

% plot orientations into the color key
hold on
plotIPDF(ebsd.orientations,'points',10,'MarkerSize',1,'MarkerFaceColor','w','MarkerEdgeColor','w')
hold off
%%
% observe how in the inverse pole figure the orientations are scattered
% closely around the white center. Together with the fact that the
% transition from white to color is quite rapidly, this gives a high
% contrast.

%% The axis angle color key
%

[grains,ebsd.grainId] = calcGrains(ebsd,'angle',1.5*degree);
ebsd(grains(grains.grainSize<=3)) = [];
[grains,ebsd.grainId] = calcGrains(ebsd,'angle',1.5*degree);

grains = smooth(grains,5);

%%

ipfKey = axisAngleColorKey(ebsd('indexed'));

% use for the reference orientation the grain mean orientation
ipfKey.oriRef = grains.meanOrientation(ebsd.grainId);

plot(ebsd,ipfKey.orientation2color(ebsd.orientations))

hold on
plot(grains.boundary,'lineWidth',4)
hold off

%%

F = halfQuadraticFilter;

ebsdS = smooth(ebsd,F,'fill',grains);

% use for the reference orientation the grain mean orientation
ipfKey.oriRef = grains.meanOrientation(ebsdS('indexed').grainId);

plot(ebsdS('indexed'),ipfKey.orientation2color(ebsdS('indexed').orientations))

hold on
plot(grains.boundary,'lineWidth',4)
hold off

%%

F = infimalConvolutionFilter;
F.lambda = 0.01;
F.mu = 0.02;

ebsdS = smooth(ebsd,F);

[grains,ebsdS.grainId] = calcGrains(ebsdS,'angle',1*degree);
%ebsdS(grains(grains.grainSize<=3)) = [];
%[grains,ebsdS.grainId] = calcGrains(ebsdS,'angle',1.5*degree);

grains = smooth(grains,5);

% use for the reference orientation the grain mean orientation
ipfKey.oriRef = grains(ebsdS('indexed').grainId).meanOrientation;
%ipfKey.oriRef = mean(ebsdS('indexed').orientations);

plot(ebsdS('indexed'),ipfKey.orientation2color(ebsdS('indexed').orientations),'micronbar','off')

hold on
plot(grains.boundary,'lineWidth',4)
hold off

%%
% Another example is when analyzing the orientation distribution within
% grains

mtexdata forsterite
ebsd = ebsd('indexed');

% segment grains
[grains,ebsd.grainId] = calcGrains(ebsd);

% find largest grains
largeGrains = grains(grains.grainSize > 800)

ebsd = ebsd(largeGrains(1))

%%
% When plotting one specific grain with its orientations we see that they
% all are very similar and, hence, get the same color

% plot a grain
close all
plot(largeGrains(1).boundary,'linewidth',2)
hold on
plot(ebsd,ebsd.orientations)
hold off

%%
% when applying the option sharp MTEX colors the mean orientation as white
% and scales the maximum saturation to fit the maximum misorientation
% angle. This way deviations of the orientation within one grain can be
% visualized.

% plot a grain
plot(largeGrains(1).boundary,'linewidth',2)
hold on
ipfKey = ipfHSVKey(ebsd);
ipfKey.inversePoleFigureDirection = mean(ebsd.orientations) * ipfKey.whiteCenter;
ipfKey.maxAngle = 2*degree;
plot(ebsd,ipfKey.orientation2color(ebsd.orientations))
hold off

author: mtex-toolbox

repo: mtex

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


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%

for i=1:p

   X_poly(:,i) = X.^i;
end

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

end

author: zzlyw

repo: machine-learning-exercises


text:
%% subfunction: set slice timing related info
function [h, hdr] = sliceTiming(h, hdr)
s = h{1};
TR = tryGetField(s, 'RepetitionTime'); % in ms
if isempty(TR), TR = tryGetField(s, 'TemporalResolution'); end
if isempty(TR), return; end
hdr.pixdim(5) = TR / 1000;
hdr.xyzt_units = 8; % seconds
if hdr.dim(5)<3 || tryGetField(s, 'isDTI', 0) || ...
        strncmp(tryGetField(s, 'MRAcquisitionType'), '3D', 2)
    return; % skip 3D, DTI, fieldmap, short EPI etc
end

nSL = hdr.dim(4);
delay = asc_header(s, 'lDelayTimeInTR', 0)/1000; % in ms now
if delay ~= 0, h{1}.DelayTimeInTR = delay; end
TA = TR - delay;

% Siemens mosaic
t = csa_header(s, 'MosaicRefAcqTimes'); % in ms
if ~isempty(t) && isfield(s, 'LastFile') && max(t)-min(t)>TA % MB wrong vol 1
    try t = mb_slicetiming(s, TA); end %#ok<*TRYNC>
end

if isempty(t) && strncmpi(s.Manufacturer, 'UIH', 3)
    t = zeros(nSL, 1);
    if isfield(s, 'MRVFrameSequence') % mosaic
        for j = 1:nSL
            item = sprintf('Item_%g', j);
            str = s.MRVFrameSequence.(item).AcquisitionDateTime;
            t(j) = datenum(str, 'yyyymmddHHMMSS.fff');
        end
    else
        dict = dicm_dict('', 'AcquisitionDateTime');
        for j = 1:nSL
            s1 = dicm_hdr(h{j}.Filename, dict);
            t(j) = datenum(s1.AcquisitionDateTime, 'yyyymmddHHMMSS.fff');
        end
    end
    t = (t - min(t)) * 24 * 3600 * 1000; % day to ms
end

if isempty(t) && any(isfield(s, {'TriggerTime' 'RTIA_timer'})) % GE
    ind = numel(h) + (1-nSL:0); % seen problem for 1st vol, so use last vol
    t = cellfun(@(c)tryGetField(c, 'TriggerTime', 0), h(ind));
    if all(diff(t)==0), t = cellfun(@(c)tryGetField(c, 'RTIA_timer', 0), h(ind)); end
    if all(diff(t)==0), t = [];
    else
        t = t - min(t);
        ma = max(t) / TA;
        if ma>1, t = t / 10; % was ms*10, old dicom
        elseif ma<1e-3, t = t * 1000; % was sec, new dicom?
        end
    end
end

if isempty(t) && isfield(s, 'ProtocolDataBlock') && ...
        isfield(s.ProtocolDataBlock, 'SLICEORDER') % GE with invalid RTIA_timer
    SliceOrder = s.ProtocolDataBlock.SLICEORDER;
    t = (0:nSL-1)' * TA/nSL;
    if strcmp(SliceOrder, '1') % 0/1: sequential/interleaved based on limited data
        t([1:2:nSL 2:2:nSL]) = t;
    elseif ~strcmp(SliceOrder, '0')
        errorLog(['Unknown SLICEORDER (' SliceOrder ') for ' s.Filename]);
        return;
    end
end

% Siemens multiframe: read TimeAfterStart from last file
if isempty(t) && tryGetField(s, 'NumberOfFrames', 1)>1 &&  ...
        ~isempty(csa_header(s, 'TimeAfterStart'))
    % Use TimeAfterStart, not FrameAcquisitionDatetime. See
    % https://github.com/rordenlab/dcm2niix/issues/240#issuecomment-433036901
    % s2 = struct('FrameAcquisitionDatetime', {cell(nSL,1)});
    % s2 = dicm_hdr(h{end}, s2, 1:nSL); % avoid 1st volume
    % t = datenum(s2.FrameAcquisitionDatetime, 'yyyymmddHHMMSS.fff');
    % t = (t - min(t)) * 24 * 3600 * 1000; % day to ms
    s2 = struct('TimeAfterStart', nan(1, nSL));
    s2 = dicm_hdr(h{end}, s2, 1:nSL); % avoid 1st volume
    t = s2.TimeAfterStart; % in secs
    t = (t - min(t)) * 1000;
end

% Get slice timing for non-mosaic Siemens file. Could remove Manufacturer
% check, but GE/Philips AcquisitionTime seems useless
if isempty(t) && ~tryGetField(s, 'isMos', 0) && strncmpi(s.Manufacturer, 'SIEMENS', 7)
    dict = dicm_dict('', {'AcquisitionDateTime' 'AcquisitionDate' 'AcquisitionTime'});
    t = zeros(nSL, 1);
    for j = 1:nSL
        s1 = dicm_hdr(h{j}.Filename, dict);
        try str = s1.AcquisitionDateTime;
        catch
            try str = [s1.AcquisitionDate s1.AcquisitionTime];
            catch, t = []; break;
            end
        end
        t(j) = datenum(str, 'yyyymmddHHMMSS.fff');
    end
    t = (t - min(t)) * 24 * 3600 * 1000; % day to ms
end

if isempty(t) % non-mosaic Siemens: create 't' based on ucMode
    ucMode = asc_header(s, 'sSliceArray.ucMode'); % 1/2/4: Asc/Desc/Inter
    if isempty(ucMode), return; end
    t = (0:nSL-1)' * TA/nSL;
    if ucMode==2
        t = t(nSL:-1:1);
    elseif ucMode==4
        if mod(nSL,2), t([1:2:nSL 2:2:nSL]) = t;
        else, t([2:2:nSL 1:2:nSL]) = t;
        end
    end
    if asc_header(s, 'sSliceArray.ucImageNumb'), t = t(nSL:-1:1); end % rev-num
end

if numel(t)<2, return; end
t = t - min(t); % it may be relative to 1st slice

t1 = sort(t);
if t1(1)==t1(2) || (t1(end)>TA), sc = 0; % no useful info, or bad timing MB
elseif t1(1) == t1(2), sc = 0; t1 = unique(t1); % was 7 for MB but error in FS
elseif isequal(t, t1), sc = 1; % ascending
elseif isequal(t, flip(t1)), sc = 2; % descending
elseif t(1)<t(3) % ascending interleaved
    if t(1)<t(2), sc = 3; % odd slices first
    else, sc = 5; % Siemens even number of slices
    end
elseif t(1)>t(3) % descending interleaved
    if t(1)>t(2), sc = 4;
    else, sc = 6; % Siemens even number of slices
    end
else, sc = 0; % unlikely to reach
end

h{1}.SliceTiming = 0.5 - t/TR; % as for FSL custom timing
hdr.slice_code = sc;
hdr.slice_end = nSL-1; % 0-based, slice_start default to 0
hdr.slice_duration = min(diff(t1))/1000;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function L=lebesque(x)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Computes the Lebesgue coefficient for a set of nodes on an interval
%
% References:
%
% (1) Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange
%     Interpolation"
%     http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/berrut.ps.gz
% (2) Walter Gaustschi, "Numerical Analysis, An Introduction" (1997) p. 76+
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get number of nodes
n=length(x);

% use this finer mesh for interpolating between nodes
N=8*n+1;

x=x(:);

X=repmat(x,1,n);

% Compute the weights
w=1./prod(X-X'+eye(n),2);

% Fine mesh for interpolating
xp=linspace(min(x),max(x),N)';

xdiff=repmat(xp.',n,1)-repmat(x,1,N);

% find all the points where the difference is zero
zerodex=(xdiff==0);

% See eq. 3.1 in Ref (1)
lfun=prod(xdiff,1);

% kill zeros
xdiff(zerodex)=eps;

% Compute lebesgue function
lebfun=sum(abs(diag(w)*repmat(lfun,n,1)./xdiff));
plot(xp,lebfun);
L=max(lebfun);


text:
function [en,S] = DSAFadapt_Morgan_cllp(un,dn,S)

% DSAFadapt_Morgan_cllp     Closed-Loop Delayless Subband Adaptive Filter Using the
%                           Weight Transformation Proposed by Morgan and Thi (Section 4.4.3.1)
% Arguments:
% un                        Input signal
% dn                        Desired signal
% S                         Adptive filter parameters as defined in MSAFinit.m
% en                        History of error signal
%
% Note:     Only N is even is implemented, where the conjugate symmetric of H_i(z)
%           and H_{N-1}(z) is exploited!!!
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

M = length(S.FULLcoeffs);
Ms = size(S.SUBcoeffs,1);
[L,N] = size(S.analysis);
D = S.decfac;
mu = S.step;
alpha = S.alpha;
AdaptStart = S.AdaptStart;
UpdateRate = S.UpdateRate;
H = S.analysis;
F = S.synthesis;
W = S.SUBcoeffs;                   % Adaptive subfilters
U = zeros(size(W));                % Tapped-delay lines of adaptive subfilters
w = S.FULLcoeffs;                  % Fullband adaptive filter
u = zeros(M,1);                    % Fullband tap-input vector
H = H(:,1:N/2+1);                  % Analysis filter (only the filrst N/2+1 is taken)
x = zeros(L,1);                    % Tapped-delay line of input signal (Analysis FB)
e = zeros(L,1);                    % Tapped-delay line of error signal (Analysis FB)

ITER = length(un);
en = zeros(1,ITER);

if isfield(S,'unknownsys')
    b = S.unknownsys;
    norm_b = norm(b);
    eml = zeros(1,ITER);
    ComputeEML = 1;
else
    ComputeEML = 0;
end

for n = 1:ITER
    x = [un(n); x(1:end-1)];       % Fullband input vector for subband filtering
    u = [un(n); u(1:end-1)];

    if ComputeEML == 1;
        eml(n) = norm(b-w)/norm_b; % System error norm (normalized)
    end

    en(n) = dn(n) - w.'*u;         % Error signal (open loop structure)
    e = [en(n); e(1:end-1)];

    if (mod(n-1,D)==0)             % Tap-weight adaptation at lowest sampling rate
        U = [x.'*H; U(1:end-1,:)]; % Each colums hold a subband regression vector
        eD = e.'*H;                % Row vector, each element is subband estimtion error
        if n >= AdaptStart
            W = W + conj(U)*diag(eD./(sum(U.*conj(U))+alpha))*mu;
            S.iter = S.iter + 1;
        end
    end
    if (n >= AdaptStart)& (mod(n-1,UpdateRate)==0)
        w = WeightTransform_Morgan(W,N,D,Ms);
                                   % Weight transformation (modified Morgan's method)
    end
end

S.FULLcoeffs = w;
S.SUBcoeffs = W;
if ComputeEML == 1;
    S.eml = eml;
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function xfat = r8vec_expand_linear ( n, x, fat )

%*****************************************************************************80
%
%% R8VEC_EXPAND_LINEAR linearly interpolates new data into a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of input data values.
%
%    Input, real X(N), the original data.
%
%    Input, integer FAT, the number of data values to interpolate
%    between each pair of original data values.
%
%    Output, real XFAT((N-1)*(FAT+1)+1), the "fattened" data.
%
  k = 0;

  for i = 1 : n-1

    k = k + 1;
    xfat(k) = x(i);

    for j = 1 : fat
      k = k + 1;
      xfat(k) = ( ( fat - j + 1 ) * x(i)     ...
                + (       j     ) * x(i+1) ) ...
                / ( fat     + 1 );
    end

  end

  k = k + 1;
  xfat(k) = x(n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% run_simExampleBasisSet.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% This script is run simply by editing the input parameters and then
% clicking "Run".
%
% DESCRIPTION:
% Script to generate simulated basis spectra for all metabolites of interest
% in the human brain.  The script will generate an LCModel format .RAW file
% for each metabolite basis spectrum, which can then be passed into
% LCModel's "makebasis" function to generate a complete LCModel basis set.
%
% INPUTS:
% To run this script, edit the following parameters as desired and then click
% run:
%  lb       = linewidth (Hz)
%  np       = Spectral points
%  sw       = Spectral width (Hz)
%  Bo       = Magnetic Field Strength (Tesla)
%  te1      = First PRESS echo time, or SPECIAL echo time (ms)
%  te2      = Second PRESS echo time (if applicable) (ms).
%  seq      = Pulse sequence ('se'= SPECIAL, 'p'=press, 'st'=steam);
%  ref      = Add reference peak at 0ppm (used in LCModel, y or n);
%
% OUTPUTS:
% H2O       = Simulated water spectrum
% Ala       = Simulated alanine spectrum
% Asp       = Simulated aspartate spectrum
% PCh       = Simulated phosphocholine spectrum
% Cr        = Simulated creatine spectrum
% PCr       = Simulated phosphochreatine spectrum
% GABA      = Simulated GABA spectrum
% Gln       = Simulated glutamine spectrum
% Glu       = Simulated glutamate spectrum
% GSH       = Simulated glutathione spectrum
% Gly       = Simulated glycine spectrum
% Ins       = Simulated myo-inositol spectrum
% Lac       = Simulated Lactate spectrum
% NAA       = Simulated NAA spectrum
% Scyllo    = Simulated scyllo-inositol spectrum
% Tau       = Simulated taurine spectrum
% Asc       = Simulated ascorbate spectrum
% bHB       = Simulated beta-hydroxybutyrate spectrum
% bHG       = Simulated 2-hydroxyglutyrate spectrum
% Glc       = Simulated glucose spectrum
% NAAG      = Simulated N-acetylaspartylglutamate spectrum
% GPC       = Simulated glycerophosphocholine spectrum
% PE        = Simulated phosphoethanolamine spectrum
% Ser       = Simulated serine spectrum
% EtOH      = Simulated ethanol spectrum
%
%
% ************INPUT PARAMETERS**********************************
 lb=2;         %linewidth (Hz)
 np=8192;      %Spectral points
 sw=4000;      %Spectral width (Hz)
 Bo=7;         %Magnetic Field Strength (Tesla)
 te1=10;       %First PRESS echo time, or SPECIAL echo time (ms)
 te2=125;      %Second PRESS echo time (if applicable) (ms).
 seq='p'       %Pulse sequence ('se'= SPECIAL, 'p'=press, 'st'=steam, 'l'=laser);
 ref='n'       %Add reference peak at 0ppm (used in LCModel, y or n);
% *************END OF INPUT PARAMETERS**************************

    [RF,H2O]=sim_lcmrawbasis(np,sw,Bo,lb,'H2O',te1,te2,ref,'y',seq);
    [RF,Ala]=sim_lcmrawbasis(np,sw,Bo,lb,'Ala',te1,te2,ref,'y',seq);
    [RF,Asp]=sim_lcmrawbasis(np,sw,Bo,lb,'Asp',te1,te2,ref,'y',seq);
    [RF,PCh]=sim_lcmrawbasis(np,sw,Bo,lb,'PCh',te1,te2,ref,'y',seq);
    [RF,Cr]=sim_lcmrawbasis(np,sw,Bo,lb,'Cr',te1,te2,ref,'y',seq);
    [RF,PCr]=sim_lcmrawbasis(np,sw,Bo,lb,'PCr',te1,te2,ref,'y',seq);
    [RF,GABA]=sim_lcmrawbasis(np,sw,Bo,lb,'GABA',te1,te2,ref,'y',seq);
    [RF,Gln]=sim_lcmrawbasis(np,sw,Bo,lb,'Gln',te1,te2,ref,'y',seq);
    [RF,Glu]=sim_lcmrawbasis(np,sw,Bo,lb,'Glu',te1,te2,ref,'y',seq);
    [RF,GSH]=sim_lcmrawbasis(np,sw,Bo,lb,'GSH',te1,te2,ref,'y',seq);
    [RF,Gly]=sim_lcmrawbasis(np,sw,Bo,lb,'Gly',te1,te2,ref,'y',seq);
    [RF,Ins]=sim_lcmrawbasis(np,sw,Bo,lb,'Ins',te1,te2,ref,'y',seq);
    [RF,Lac]=sim_lcmrawbasis(np,sw,Bo,lb,'Lac',te1,te2,ref,'y',seq);
    [RF,NAA]=sim_lcmrawbasis(np,sw,Bo,lb,'NAA',te1,te2,ref,'y',seq);
    [RF,Scyllo]=sim_lcmrawbasis(np,sw,Bo,lb,'Scyllo',te1,te2,ref,'y',seq);
    [RF,Tau]=sim_lcmrawbasis(np,sw,Bo,lb,'Tau',te1,te2,ref,'y',seq);
    [RF,Asc]=sim_lcmrawbasis(np,sw,Bo,lb,'Asc',te1,te2,ref,'y',seq);
    [RF,bHB]=sim_lcmrawbasis(np,sw,Bo,lb,'bHB',te1,te2,ref,'y',seq);
    [RF,bHG]=sim_lcmrawbasis(np,sw,Bo,lb,'bHG',te1,te2,ref,'y',seq);
    [RF,Glc]=sim_lcmrawbasis(np,sw,Bo,lb,'Glc',te1,te2,ref,'y',seq);
    [RF,NAAG]=sim_lcmrawbasis(np,sw,Bo,lb,'NAAG',te1,te2,ref,'y',seq);
    [RF,GPC]=sim_lcmrawbasis(np,sw,Bo,lb,'GPC',te1,te2,ref,'y',seq);
    [RF,PE]=sim_lcmrawbasis(np,sw,Bo,lb,'PE',te1,te2,ref,'y',seq);
    [RF,Ser]=sim_lcmrawbasis(np,sw,Bo,lb,'Ser',te1,te2,ref,'y',seq);
    [RF,EtOH]=sim_lcmrawbasis(np,sw,Bo,lb,'EtOH',te1,te2,ref,'y',seq);

%LEGEND:
%   'Ala'    = Alanine
%   'Asp'    = Aspartate
%   'PCh'    = PhosphoCholine
%   'Cr'     = Creatine
%   'PCr'    = PhosphoCreatine
%   'GABA'   = Gamma-aminobutyric acid
%   'Gln'    = Glutamine
%   'Glu'    = Glutamate
%   'GSH'    = Glutathione
%   'Gly'    = Glycine
%   'Ins'    = Myo-inositol
%   'Lac'    = Lactate
%   'NAA'    = N-acetyl aspartate
%   'Scyllo' = Scyllo-inositol
%   'Tau'    = Taurine
%   'Asc'    = Ascorbate (Vitamin C)
%   'bHB'    = beta-Hydroxybutyrate
%   'bHG'    = beta-Hydroxyglutarate
%   'Glc'    = Glucose
%   'NAAG'   = N-acetyl aspartyl glutamate
%   'GPC'    = Glycero-phosphocholine
%   'PE'     = Phosphoryl ethanolamine
%   'Ser'    = Serine
%   'EtOH'   = Ethanol

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%%%  e8.2_Sim
clear
path('./icon/',path);
%Run the model initialization file icon/init.m
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV hover
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE =6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE =6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_ANGLE = 3;
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;

%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 3.0; Kvxi = 0.1; Kvxd = 0.01;
Kvyp = 3.0; Kvyi = 0.1; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;

Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%Maximum navigation angle, used in position control to avoid large attitude angles
MAX_CONTROL_ANGLE_NAV_ROLL = 15;
MAX_CONTROL_ANGLE_NAV_PITCH  = 15;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
e8_2_sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [c, r, t, cpulse, acq_codes, verbose, gsr] = tapas_physio_read_physlogfiles_biopac_txt(...
    log_files, cardiac_modality, verbose, varargin)
% Reads in 4-column txt-export from BioPac Data (resp, ppu, gsr, trigger)
%
% [c, r, t, cpulse, acq_codes, verbose] = tapas_physio_read_physlogfiles_biopac_txt(...
%    log_files, cardiac_modality, verbose, varargin)
%
% IN    log_files
%       .log_cardiac        *.txt file, contains 4 columns of the form
%                           RESP - RSP100C	GSR - EDA100C-MRI	PPG - PPG100C	Marker
%                             -0.949402	-0.00610382	0.0134277	0
%                             -0.949402	-0.00610382	0.0134277	0
%                             -0.951233	-0.00915558	0.0204468	0
%                             -0.951233	-0.00915558	0.0204468	0
%                             -0.953064	-0.0122073	0.0259399	0
%                             -0.953064	-0.0122073	0.0259399	0
%                             -0.95459	-0.0076297	0.0296021	0
%                             -0.95459	-0.0076297	0.0296021	0
%       .log_respiration    same as .log_cardiac
%       .sampling_interval  sampling interval (in seconds)
%                           default: 1 ms (1000 Hz)
%       cardiac_modality    'ECG' or 'PULS'/'PPU'/'OXY' to determine
%                           which channel data to be returned
%                           UNUSED, is always pulse plethysmographic unit
%                           for BioPac
%       verbose
%       .level              debugging plots are created if level >=3
%       .fig_handles        appended by handle to output figure
%
% OUT
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%                       <UNUSED>, since not written to logfile
%   r                   respiratory time series
%   t                   vector of time points (in seconds)
%   c                   cardiac time series (PPU)
%   acq_codes           slice/volume start events marked by number <> 0
%                       for time points in t
%                       10/20 = scan start/end;
%                       1 = ECG pulse; 2 = OXY max; 4 = Resp trigger;
%                       8 = scan volume trigger (on)
%                       16 = scan volume trigger (off)
%   gsr                 galvanic skin response (not used)
%
% EXAMPLE
%   tapas_physio_read_physlogfiles_biopac_txt
%
%   See also tapas_physio_read_physlogfiles_siemens tapas_physio_plot_raw_physdata_siemens_hcp

% Author: Lars Kasper
% Created: 2018-09-27
% Copyright (C) 2018 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.

% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% read out values
DEBUG = verbose.level >= 2;

hasRespirationFile = ~isempty(log_files.respiration);
hasCardiacFile = ~isempty(log_files.cardiac);

hasRespirationFile = ~isempty(log_files.respiration);
hasCardiacFile = ~isempty(log_files.cardiac);

if hasCardiacFile
    fileName = log_files.cardiac;
elseif hasRespirationFile
    fileName = log_files.respiration;
end

[C, columnNames] = tapas_physio_read_columnar_textfiles(fileName, 'BIOPAC_TXT');
c = double(C{3});
r = double(C{1});
gsr = double(C{2});
iAcqOn = (double(C{4})~=0); % trigger has 11, rest is 0;

%% Create timing vector from samples

dt = log_files.sampling_interval;

if isempty(dt)
    dt = 1/1000; % 1000 Hz sampling interval
end

nSamples = max(numel(c), numel(r));
t = -log_files.relative_start_acquisition + ((0:(nSamples-1))*dt)';

%% Recompute acq_codes as for Siemens (volume on/volume off)
acq_codes = [];

if ~isempty(iAcqOn) % otherwise, nothing to read ...
    % iAcqOn is a column of 1s and 0s, 1 whenever scan acquisition is on
    % Determine 1st start and last stop directly via first/last 1
    % Determine everything else in between via difference (go 1->0 or 0->1)
    iAcqStart   = find(iAcqOn, 1, 'first');
    iAcqEnd     = find(iAcqOn, 1, 'last');
    d_iAcqOn    = diff(iAcqOn);

    % index shift + 1, since diff vector has index of differences i_(n+1) - i_n,
    % and the latter of the two operands (i_(n+1)) has sought value +1
    iAcqStart   = [iAcqStart; find(d_iAcqOn == 1) + 1];
    % no index shift, for the same reason
    iAcqEnd     = [find(d_iAcqOn == -1); iAcqEnd];

    acq_codes = zeros(nSamples,1);
    acq_codes(iAcqStart) = 8; % to match Philips etc. format
    acq_codes(iAcqEnd) = 16; % don't know...

    % report estimated onset gap between last slice of volume_n and 1st slice of
    % volume_(n+1)
    nAcqStarts = numel(iAcqStart);
    nAcqEnds = numel(iAcqEnd);
    nAcqs = min(nAcqStarts, nAcqEnds);

    if nAcqs >= 1
        % report time of acquisition, as defined in SPM
        TA = mean(t(iAcqEnd(1:nAcqs)) - t(iAcqStart(1:nAcqs)));
        verbose = tapas_physio_log(...
            sprintf('TA = %.4f s (Estimated time of acquisition during one volume TR)', ...
            TA), verbose, 0);
    end
end

%% Plot, if wanted

if DEBUG
    stringTitle = 'Read-In: Raw BioPac physlog data (TXT Export)';
    verbose.fig_handles(end+1) = ...
        tapas_physio_plot_raw_physdata_siemens_hcp(t, c, r, acq_codes, ...
        stringTitle);
end

%% Undefined output parameters

cpulse = [];

author: translationalneuromodeling

repo: tapas


text:
function visualizeBoundary(X, y, model, varargin)
%VISUALIZEBOUNDARY plots a non-linear decision boundary learned by the SVM
%   VISUALIZEBOUNDARYLINEAR(X, y, model) plots a non-linear decision
%   boundary learned by the SVM and overlays the data on it

% Plot the training data on top of the boundary
plotData(X, y)

% Make classification predictions over a grid of values
x1plot = linspace(min(X(:,1)), max(X(:,1)), 100)';
x2plot = linspace(min(X(:,2)), max(X(:,2)), 100)';
[X1, X2] = meshgrid(x1plot, x2plot);
vals = zeros(size(X1));
for i = 1:size(X1, 2)
   this_X = [X1(:, i), X2(:, i)];
   vals(:, i) = svmPredict(model, this_X);
end

% Plot the SVM boundary
hold on
contour(X1, X2, vals, [0 0], {'LineWidth', 1, 'LineColor', 'b'});
hold off;

end

author: vugsus

repo: coursera-machine-learning


text:
function handles = tight_subplots(rows, columns, gap, row_margin, column_margin)
% tight_subplots Initializes a grid of axes
%
% Initializes a grid of axes with adjustable gaps and margins
% and returns their handles.
%
% Credit: Pekka Kumpulainen (2010) @tut.fi Tampere University of Technology / Automation Science and Engineering
%
% Input:
% - rows (integer): Number of axes in hight (vertical direction).
% - columns (integer): Number of axes in width (horizontal direction).
% - gap (double): Gaps between the axes in normalized units (0...1)
%   or [gap_h gap_w] for different gaps in height and width.
% - row_margin (<type>): Margins in height in normalized units (0...1)
%   or [lower upper] for different lower and upper margins.
% - column_margin (<type>): Margins in width in normalized units (0...1)
%   or [left right] for different left and right margins
%
% Output:
% - handles (matrix): Matrix of handles for the axes objects
%   starting from upper left corner, going row-wise as in
%   going row-wise as in.
%

if nargin<3; gap = .02; end
if nargin<4 || isempty(row_margin); row_margin = .05; end
if nargin<5; column_margin = .05; end

if numel(gap)==1;
    gap = [gap gap];
end
if numel(column_margin)==1;
    column_margin = [column_margin column_margin];
end
if numel(row_margin)==1;
    row_margin = [row_margin row_margin];
end

axh = (1-sum(row_margin)-(rows-1)*gap(1))/rows;
axw = (1-sum(column_margin)-(columns-1)*gap(2))/columns;

py = 1-row_margin(2)-axh;

handles = zeros(rows*columns,1);
ii = 0;
for ih = 1:rows
    px = column_margin(1);

    for ix = 1:columns
        ii = ii+1;
        handles(ii) = axes('Units','normalized', ...
            'Position',[px py axw axh], ...
            'XTickLabel','', ...
            'YTickLabel',''); %#ok<LAXES>
        px = px+axw+gap(2);
    end
    py = py-axh-gap(1);
end

author: votchallenge

repo: toolkit-legacy


text:
function a = r8vec_even ( n, alo, ahi )

%*****************************************************************************80
%
%% R8VEC_EVEN returns N real values, evenly spaced between ALO and AHI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values.
%
%    Input, real ALO, AHI, the low and high values.
%
%    Output, real A(N), N evenly spaced values.
%    Normally, A(1) = ALO and A(N) = AHI.
%    However, if N = 1, then A(1) = 0.5*(ALO+AHI).
%
  if ( n == 1 )

    a(1) = 0.5 * ( alo + ahi );

  else

    a(1:n) = ( (n-1:-1:0) * alo + (0:n-1) * ahi ) / ( n - 1 );

  end

  return
end


text:
%% Defining Tensorial Properties
%
% Physical laws describe the relationship between physical properties. The
% most simplest laws are linear ones and are of the form
%
% $$ y = \mathbf A x $$
%
% where $x$ and $y$ are the physical properties and $\mathbf A$ is a
% material constant. In a typical example $y$ could be the force applied to
% a spring, $x$ the displacement and $A$ describes the stiffnes of the
% spring which is essentially Hooks law.
%
% As soon as we consider more general forces and displacements they can not
% be described anymore by scalar numbers $x$ and $y$ but vectors or
% matrices are required. In its most general form the displacment is
% describes by a <strainTensor.html strain matrix> $\sigma_{ij}$ and the
% force is described by a stiffness matrix $\varepsilon_{kl}$. In this setting
% the linear relationship between the two matrices is described by the
% <complianceTensor.complianceTensor.html compliance tensor> $\mathbf
% C_{ijkl}$ which can be seen as a 4 dimensional generalization of a matrix.
%
% More, general a tensor of rank $r$ is a "matrix" of dimension $r$. If
% $r=0$ we speek of scalars, if $r=1$ these are vectors and for $r=2$ they
% are classical $3 \times 3$ matrices.
%
% In the following we explain how tensors of arbitrary rank can be defined
% in MTEX. Independent of the rank any tensor is represented in MTEX by a
% variable of type @tensor.
%
%
%% Scalars (tensors of zero rank)
%
% In physics, properties like temperature or density are not connected
% to any specific direction of the body they are measured. These
% non-directional physical quantities are called scalars, and they are
% defined by a single number. In MTEX, a scalar is defined by:

M = 5;
t = tensor(M,'rank',0)

%% Vectors (tensors of first rank)
%
% In contrast to scalars, other physical quantities can only be defined in
% reference to a specific direction. If we need to specify completely the
% mechanical force acting into a point for example, we need to specify
% the magnitude and its direction. As an alternative, we can choose three
% mutually perpendicular axes (A1,A2 and A3) and give the vector components
% along them. In MTEX, this is done by:

t = tensor([1;2;3],'rank',1)

%%
% where 1, 2 and 3 are the components related to the axes A1, A2 and A3.
% As rank 1 tensors are essentialy vectors we can freely convert tensors to
% @vector3d and vice verca.

% define a tensor from a vector
t = tensor(vector3d.X)

% convert a tensor into a vector
vector3d(t)

%% Matrices (tensors of second rank)
%
% We have now to expand the idea of a vector to three-dimensional space.
% Let's take the example of stress (force per unit of area). Imagine a cube
% of material subjected to load as shown below. As can be seen, one can
% measure ther stresses in this cube in various directions, and in various
% planes. These measurements will for a second rank sensor, where each
% component is associated with a pair of axes, taken in an specific order.
% The generalized second rank stress tensor can be written as
%
% $$
% \sigma_{ij} =
% \left[\begin{array}{ccc}
% \sigma_{11} & \sigma_{12} & \sigma_{13}  \\
% \sigma_{21} & \sigma_{22} & \sigma_{23}  \\
% \sigma_{31} & \sigma_{32} & \sigma_{33}  \\
% \end{array}\right]
% $$
%
% In MTEX, a second-rank tensor where only the main diagonal components are
% of interest is defined as

t = tensor(diag([1,2,3]), 'rank',2)

%%
% If all the components are of interest, the definition is as follow

M = [1     0.75  0.5;...
     0.75  1     0.25;...
     0.5   0.25  1];

t = tensor(M,'rank',2)

%% Tensors (tensors of third rank)
%
% Smart materials are materials that have one or more properties that
% change significantly under external stimuli. A typical example is the
% voltage resulting to applied stress that certain materials have, named
% piezoeletric effect. This property is described as a third rank tensor
% that relates induced electric displacement vector to the second-order
% stress tensor. This is expressed in the form $P_i=d_{ijk} \sigma_{jk}$.
% In MTEX, a third rank tensor can be described as

M =[[-1.9222  1.9222    0   -0.1423     0         0    ];...
    [   0        0      0       0     0.1423    3.8444];...
    [   0        0      0       0       0         0    ]];

t = tensor(M,'rank',3)

%% Tensors (tensors of fourth rank)
%
% Fourth rank tensors are tensors that describe the relation between 2
% second rank tensors. A typical example is the tensor describing the
% elastic properties of materials, which translate the linear relationship
% between the second rank stress and infinitesimal strain tensors. The
% Hooke's Law describing the shape changes in a material subject to stress
% can be written as $\sigma_{ij}=c_{ijkl} \epsilon_{kl}$, where $c_{ijkl}$
% is a fourth rank tensor.
%
% The four indices (ijkl) of the elastic tensor have values between 1 and
% 3, so that there are $3^4=81$ coefficients. As the stress and strain
% tensors are symmetric, both stress and strain second rank tensors only
% have 6 independent values rather than 9. In addition, crystal symmetry
% reduces even more the number of independent components on the elastic
% tensor, from 21 in the case of triclinic phases, to 3 in the case of
% cubic materials. In MTEX, a fourth rank tensor can be defined as:

M = [[320   50  50   0     0     0];...
    [  50  320  50   0     0     0];...
    [  50   50 320   0     0     0];...
    [   0    0   0  64     0     0];...
    [   0    0   0   0    64     0];...
    [   0    0   0   0     0    64]];

C = tensor(M,'rank',4)

%%
% Note the repetition in values in this matrix is related to crystal
% symmetry, in this case, a cubic example, where only $C_{11}$, $C_{12}$
% and $C_{44}$ are independent components.
%
%% Specific tensors
%
% MTEX includes specific classes for the following tensors.
%
% || *name* || *rank* || *symbol* || *name* || *rank* || *symbol* ||
% || @complianceTensor || 4 || $S_{ijkl}$ || @stiffnessTensor || 4 || $C_{ijkl}$ ||
% || @strainTensor || 2 || $\sigma_{ij}$  || @stressTensor || 2 || $\varepsilon_{ij}$ ||
% || @strainRateTensor || 2 || $E$ || @velocityGradientTensor || 2 || $L$ ||
% || @curvatureTensor || 2 || $\kappa_{ij}$ || @deformationGradientTensor || 2 || $F$ ||
% || @refractiveIndexTensor || 2 || $\chi$ || @ChristoffelTensor || 2 || $M_{ij}$ ||
% || @dislocationDensityTensor || 2 || $\alpha$  || <SchmidTensor.html |SchmidTensor|> || 2 || $M_{ij}$ ||
% || <tensor.leviCivita.html |leviCivita|> || 3 || $\varepsilon_{ijk}$ || @spinTensor || 2 || $\Omega$ ||
%
% Those specific tensors are defined by the syntax

M = [0 0 0;...
  0 0 0; ...
  0 0 1];

e = strainTensor(M)

%%
% In many cases shortcuts exist like

e = stressTensor.uniaxial(vector3d.Z)

%%
% The advantage of using these specific tensor classes is that some tensor
% operations like <stressTensor/calcShearStress.html |calcShearStress(e)|>
% are defined only for specific tensor classes.
%
%% Predefined tensors
%
% For certain applications, one may want to have a tensor where all the
% components are 1. In MTEX this is computed as

t = tensor.ones('rank',2)

%%
% *Identity tensor*
%
% The Identity tensor is a second order tensor that has ones n the main
% diagonal and zeros otherwise. The identity matrix has some special
% properties, including (i) When multiplied by itself, the result is itself
% and (ii) rows and columns are linearly independent. In MTEX, this matrix
% can be computed as

t = tensor.eye('rank',2)

%%
% *Random tensors*
%
% One can also define a tensor in which the components are pseudorandom, by
% using the function |<tensor.rand.html tensor.rand>|

t = tensor.rand('rank',2)

%%
% *The Levi Civita tensor*
%
% The Levi-Civita symbol $\epsilon_{ijk}$ is a third rank tensor and is
% defined by 0, if $i=j$, $j=k$ or $k=1$, by 1, if $(i,j,k)=(1,2,3)$,
% $(2,3,1)$ or $(3,1,2)$ and by $-1$, if $(i,j,k)=(3,2,1)$, $(1,3,2)$ or
% $(2,1,3)$. The Levi-Civita symbol allows the cross product of two vectors
% in 3D Euclidean space and the determinant of a square matrix to be
% expressed in Einstein's index notation. With MTEX the Levi Civita tensor
% is expressed as

t = tensor.leviCivita


text:
%CIRCLE Compute points on a circle
%
% CIRCLE(C, R, OPTIONS) plots a circle centred at C (1x2) with radius R on the current
% axes.
%
% X = CIRCLE(C, R, OPTIONS) is a matrix (2xN) whose columns define the
% coordinates [x,y] of points around the circumference of a circle
% centred at C (1x2) and of radius R.
%
% C is normally 2x1 but if 3x1 then the circle is embedded in 3D, and X is Nx3.
% The circle is always in the xy-plane with a z-coordinate of C(3).
%
% Options::
%  'n',N   Specify the number of points (default 50)

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math
function out = circle(centre, rad, varargin)

	opt.n = 50;

    [opt,arglist] = tb_optparse(opt, varargin);

    % compute points on circumference
	th = [0:opt.n-1]'/ opt.n*2*pi;
    x = rad*cos(th) + centre(1);
    y = rad*sin(th) + centre(2);

    % add extra row if z-coordinate is specified, but circle is always in xy plane
    if length(centre) > 2
        z = ones(size(x))*centre(3);
        p = [x y z]';
    else
        p = [x y]';
    end

    if nargout > 0
        % return now
        out = p;
        return;
    else
        % else plot the circle
        p = [p p(:,1)]; % make it close
        if numrows(p) > 2
            plot3(p(1,:), p(2,:), p(3,:), arglist{:});
        else
            plot(p(1,:), p(2,:), arglist{:});
        end
    end

author: petercorke

repo: spatialmath-matlab