%==============================================================================
% 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/
%==============================================================================
% This function manages the solution of H * dy = - dJ' for various settings
%
%   'backslash',         simply solve dy = H\rhs;
%   'mbCG',               matrixBased, plain CG
%   'mbPCG-SGS',         matrixBased, PCG with Symmetric Gauss Seidel preconditioner
%   'mbPCG-ICHOL',     matrixBased, PCG with incomplete Cholesky preconditioner
%   'mbPCG-Jacobi', matrixBased, PCG with Jacobi preconditioner
%
%   'MG-elastic',           matrixFree, PCG Multigrid for elastic regularizer
%   'CG-elastic',       plain CG with H = A as in (1) below
%   'CG-curvature',          plain CG with H = A as in (1) below
%   'PCG-elastic',      Jacobi preconditioned with H = A as in (1,) below
%   'PCG-curvature',      Jacobi preconditioned with H = A as in (1,) below
%   'PCG-hyperElastic', ???
%   'CG'                      plain CG with H full or as operator
%
% the operator H can be a
%    - matrix (matrixBased)
%    - function (coding the action ofH)
%    - struct where H describes the pieces of a complex H
%
%  (1) action of H: distance + regularizer,
%       H = P'*dr'*d2psi*dr*P + d2S
%
%      Hoperator = @(x) ...
%       H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)) ...
%       + H.d2S.d2S(x,H.omega,H.m);
%
%  (2) Jacobi (diagonal) preconditing
%      Ddiag = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr);
%      D = H.d2D.P(full(Ddiag))  +  H.d2S.diag(H.omega,H.m);
%      PC = @(x) D.\x; % Jacobi preconditioner
%==============================================================================

function [dy,solver] = solveLinearSystem(rhs,H,solver,varargin)

if nargin == 0,
  testSolver
  help(mfilename);
  runMinimalExample;
  dy = 'endOfMinimalExample';
  return;
end;

maxIterCG = 500;
tolCG     = 1e-1;

for k=1:2:length(varargin), % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

if isempty(solver) && isnumeric(H),
  warning('no solver specified, try backslash');
  solver = 'backslash';
end;

if isa(solver,'function_handle')
  dy     = solver(rhs,H,maxIterCG,tolCG);
  return;
end

if isstruct(H), % matrixFree mode, configure operator
  Hoperator = @(x) ...
    H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)) ...
    + H.d2S.d2S(x,H.omega,H.m);
end;

%------------------------------------------------------------------------------
switch solver
  %------------------------------------------------------------------------------

  % ---------------------------------------------------------------------------
  % matrix based
  % ---------------------------------------------------------------------------

  case 'backslash',
    dy = H\rhs;

  case 'mbCG',
    [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG);

  case 'mbPCG-SGS',
    L   = tril(H); % Symmetric Gauss Seidel Preconditioning,
    D   = diag(H); % L is lower, D is diagonal, U = L'
    SGS = @(x) L\(D.*(L'\x));
    [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG,SGS);

  case 'mbPCG-ICHOL',
    %L1   = cholinc(sparse(H),'0');
    L = ichol(H);
    [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG,L,L');

  case 'mbPCG-Jacobi',
    D   = diag(H); % D is diagonal
    PC = @(x) D.\x;
    [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG,PC);

    % ---------------------------------------------------------------------------
    % matrix free
    % ---------------------------------------------------------------------------

    % In all matrix free solvers the approximates Hessian is supplied by function handle.
    % Note, that the functionals have the common form
    %
    %   Jc = D + S (+ P) ==> d_2 Jc = d2D + d2S (+d2P)
    %
    % regularization and penalization (S,P) supply a struct stored in H.d2S or H.d2P
    % which already has a function handle describing the operator.
    % The  distance term (d2D) needs more care and the approximation to the Hessian
    % depends on the solver as well. Therefore the objective function does NOT supply
    % the function handle, but all variables needed to built such.
    % For multi-grid (available for elastic and curvature), only the diagonal of the distance
    % term is used
    % whereas in conjugate gradient (CG, available for elastic, curvature and hyperelastic)
    % methods we can model also the off-diagonals d2D(x) = (P'*dr'*d2psi*dr*P) * x
    %
    % The default solver for a chosen regularizer is parameterized by d2S.solver.
    % However, one can overload this setting using
    %    >> MLIR(..., 'solverNPIR','myFavoriteSolver', ... );

  case {'MG-elastic','MG'}, % multigrid for elastic regularizer
    dy = MGsolver(rhs,H);

  case {'CG-elastic','CG-curvature'},
    if isstruct(H)
      %       A         = @(x) ...
      %         H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)) ...
      %         + H.d2S.d2S(x,H.omega,H.m);
      [dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG);
    else
      [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG);
    end

  case {'PCG-elastic','PCG-curvature'}

    % operator
    %     A         = @(x) ...
    %       H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)) ...
    %       + H.d2S.d2S(x,H.omega,H.m);
    % preconditioner
    Ddiag = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr);
    D = H.d2D.P(full(Ddiag))  +  H.d2S.diag(H.omega,H.m);
    PC = @(x) D.\x; % Jacobi preconditioner
    [dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG,PC);

  case 'CG'
    if isstruct(H)
      [dy,flag,relres,iter] = pcg(H.operator,rhs,tolCG,maxIterCG);
    else
      [dy,flag,relres,iter] = pcg(H,rhs,tolCG,maxIterCG);
    end

  case {'PCG-hyperElastic'}
    M         = @(x) H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x));
    Hoperator = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
    Ddiag     = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr);
    D         = H.d2D.P(full(Ddiag))  +  H.d2S.diag(H.d2S.yc);
    Preconditioner = @(x) D.\x; % Jacobi preconditioner
    [dy,flag,relres,iter] = pcg(Hoperator,rhs,tolCG,maxIterCG,Preconditioner);

  otherwise,
    keyboard
    error(1)

    %------------------------------------------------------------------------------
end
%------------------------------------------------------------------------------

if exist('flag','var')
  flags = {
    sprintf('iter=%d>maxIter=%d but relres=%e>relResTol=%s',...
    iter,maxIterCG,relres,tolCG)
    'preconditioner is ill-conditioned'
    'stagnation: x_k=x_{k+1}'
    'scalars too small/large to continue computation'
    };
  switch flag
    case {1,2,3,4}
      fprintf('%s // PCG: %s\n',flags{flag});
    otherwise
      fprintf('%s // PCG :iter = %d of %d, relRes= %1.2e, tolCG= %1.2e\n',...
        mfilename,iter,maxIterCG,relres,tolCG);
  end
end

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

function runMinimalExample

%% prepare 3D data
setup3DbrainData
level = 4;
omega = ML{level}.omega;
m     = ML{level}.m;

viewImage('reset','viewImage','viewImage2D','colormap',bone(256),'axis','off');
imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',1e-2);
distance('reset','distance','SSD');
trafo('reset','trafo','affine3Dsparse');

[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0);
xc = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);
w0   = trafo('w0');
beta = 0; M = []; wRef = [];
[Jc,para,dJPIRsparse,HPIRsparse] = PIRobjFctn(T,Rc,omega,m,beta,M,wRef,xc,w0);

%% prepare 2D data and approximations to Hessian
setup2DhandData

viewImage('reset','viewImage','viewImage2D','colormap',bone(256),'axis','off');
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-2);
distance('reset','distance','SSD');

level = 4;
omega = ML{level}.omega;
m     = ML{level}.m;
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0);

% initialize the image model, distance measure and regularizer
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

% setup PIR objective
w0   = trafo('w0');
beta = 1e3;
M    = diag([4,4,4,4,1,1]);

trafo('reset','trafo','affine2D');
w0   = trafo('w0');
wRef = w0;
[Jc,para,dJPIR,HPIR] = PIRobjFctn(T,Rc,omega,m,beta,M,wRef,xc,w0);

% setup NPIR objective
yc   = getStaggeredGrid(omega,m);
yRef = 0*yc;

regularizer('reset','regularizer','mbElastic','alpha',1e4,'mu',1,'lambda',0);
[Jc,para,dJmbElas,HmbElas] = NPIRobjFctn(T,Rc,omega,m,yRef,yc);

regularizer('reset','regularizer','mfElastic','alpha',1e4,'mu',1,'lambda',0);
[Jc,para,dJmfElas,HmfElas] = NPIRobjFctn(T,Rc,omega,m,yRef,yc);

HxElas = @(x) HmfElas.d2D.P((HmfElas.d2D.dr'*HmfElas.d2D.d2psi*HmfElas.d2D.dr)*HmfElas.d2D.P(x)) ...
  + HmfElas.d2S.d2S(x,HmfElas.omega,HmfElas.m);

yc   = getCellCenteredGrid(omega,m);
yRef = 0*yc;

regularizer('reset','regularizer','mbCurvature','alpha',1e4);
[Jc,para,dJmbCurv,HmbCurv] = NPIRobjFctn(T,Rc,omega,m,yRef,yc);

regularizer('reset','regularizer','mfCurvature','alpha',1e4);
[Jc,para,dJmfCurv,HmfCurv] = NPIRobjFctn(T,Rc,omega,m,yRef,yc);
HxCurv = @(x) HmfCurv.d2D.P((HmfCurv.d2D.dr'*HmfCurv.d2D.d2psi*HmfCurv.d2D.dr)*HmfCurv.d2D.P(x)) ...
  + HmfCurv.d2S.d2S(x,HmfCurv.omega,HmfCurv.m);

problems = {
  'PIR-backslash'
  'PIR-sparse'
  'mbElastic-backslash'
  'mbElastic-PCG-SGS'
  'mbElastic-PCG-ICHOL'
  'mbElastic-PCG-Jacobi'
  'mbElastic-CG'
  'mfElastic-MG'
  'mfElastic-CG'
  'mfElastic-PCG'
  'mbCurvature-backslash'
  'mbCurvature-PCG-SGS'
  'mbCurvature-PCG-ICHOL'
  'mbCurvature-PCG-Jacobi'
  'mbCurvature-CG'
  'mfCurvature-PCG'
  'mfCurvature-CG'
  }

for j = 1:length(problems)

  para = {};

  switch problems{j},

    case 'PIR-backslash',
      solver = 'backslash';
      H = HPIR; b = -dJPIR'; Hx = @(y) H*y;

    case 'PIR-sparse',
      solver = 'backslash';
      H = HPIRsparse; b = -dJPIRsparse'; Hx = @(y) H*y;

      %     case 2, % PRIR with conensed matrix
      % % fprintf('test PIR with condensed matrix\n');
      % % trafo('reset','trafo','affine2Dsparse');
      % % [Jc,para,dJ,H] = PIRobjFctn(T,Rc,omega,m,beta,M,wRef,xc,w0)
      % % dy = solveLinearSystem(-dJ',H,'backslash');
      % % test2 = norm(H*dy + dJ')

    case 'mbElastic-backslash'
      solver = 'backslash';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;

    case 'mbElastic-PCG-SGS',
      solver = 'mbPCG-SGS';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbElastic-PCG-ICHOL',
      solver = 'mbPCG-ICHOL';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbElastic-PCG-Jacobi'
      solver = 'mbPCG-Jacobi';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbElastic-CG'
      solver = 'mbCG';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mfElastic-MG'
      solver = 'MG-elastic';
      H = HmfElas; b = -dJmfElas'; Hx = HxElas;

    case 'mfElastic-CG'
      solver = 'CG-elastic';
      H = HmfElas; b = -dJmfElas'; Hx = HxElas;

    case 'mfElastic-PCG'
      solver = 'PCG-elastic';
      H = HmfElas; b = -dJmfElas'; Hx = HxElas;

    case 'mbCurvature-backslash'
      solver = 'backslash';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;

    case 'mbCurvature-PCG-SGS',
      solver = 'mbPCG-SGS';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbCurvature-PCG-ICHOL',
      solver = 'mbPCG-ICHOL';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbCurvature-PCG-Jacobi'
      solver = 'mbPCG-Jacobi';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mbCurvature-CG'
      solver = 'mbCG';
      H = HmbElas; b = -dJmbElas'; Hx = @(y) H*y;
      para = {'tolCG',1e-5};

    case 'mfCurvature-PCG'
      solver = 'PCG-curvature';
      H = HmfElas; b = -dJmfElas'; Hx = HxElas;

    case 'mfCurvature-CG'
      solver = 'CG-curvature';
      H = HmfElas; b = -dJmfElas'; Hx = HxElas;

    case 'MG'

    otherwise,
      error('1');
  end;

  dy = solveLinearSystem(b,H,solver,para{:});
  test = norm(Hx(dy)-b)/norm(b);
  fprintf('%-2d of %2d: %-40s test = %s\n',j,length(problems),...
    sprintf('test <%s>:',problems{j}),num2str(test));
end;

return;
%==============================================================================

%   case 'mfCurvature',
%     Afun     = @(dy) mfAy(dy,H);
%     [dy,FLAG] = pcg(Afun,rhs,tolCG,maxIterCG);
%     %

%
%   if isnumeric(H),
%
%   elseif isstruct(H),
%     if isfield(H,'solver'),
%       solver = H.solver;
%     elseif isfield(H,'d2S') && isfield(H.d2S,'solver'),
%       solver = H.d2S.solver;
%     else
%       error('solver has not been defined')
%     end;
%   else
%     keyboard
%   end;
%   fprintf('[set solver to <%s> in %s]\n',solver,mfilename);
% end;
% if isempty(solver)
%   if isstruct(H),
%   else % no regularizer initialized, assuming PIR
%     dy = H\rhs;
%     return;
%   end;
% end;
%
% if isa(solver, 'function_handle')
%     dy = feval(solver, rhs, H, maxIterCG, tolCG);
%     return
% end

%
%     %     case {'Joint-CG-hyperElastic'}
%     %         Afctn = @(x) H.M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
%     %         [dy,flag,relres,iter] = pcg(Afctn,rhs,tolCG,maxIterCG);
%     %     case {'CG-hyperElastic'}
%     %         M =@(x) H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x));
%     %         Afctn = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
%     %         [dy,flag,relres,iter] = pcg(Afctn,rhs,tolCG,maxIterCG);
%     %     otherwise
%     %         if isnumeric(H)
%     %             % if H is a matrix, solve the linear system using MATLAB's backslash
%     %             dy = H\rhs;
%     %         else
%     %             error(solver)
%     %         end
%   otherwise,
%     keyboard
%     error(1)
%     %
%   case {'PCG-hyperElastic'}
%     M =@(x) H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x));
%     Afctn = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
%     Ddiag = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr);
%     D = H.d2D.P(full(Ddiag))  +  H.d2S.diag(H.d2S.yc);
%     PC = @(x) D.\x; % Jacobi preconditioner
%     [dy,flag,relres,iter] = pcg(Afctn,rhs,tolCG,maxIterCG,PC);
%
%     %   case {'PCG-elastic','PCG-curvature'}
%     %     % operator
%     %     A         = @(x) ...
%     %       H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x)) ...
%     %       + H.d2S.d2S(x,H.omega,H.m);
%     %     % preconditioner
%     %     Ddiag = diag(H.d2D.dr'*H.d2D.d2psi*H.d2D.dr);
%     %     D = H.d2D.P(full(Ddiag))  +  H.d2S.diag(H.omega,H.m);
%     %     PC = @(x) D.\x; % Jacobi preconditioner
%     %     [dy,flag,relres,iter] = pcg(A,rhs,tolCG,maxIterCG,PC);
%
%
%     %     case {'Joint-CG-hyperElastic'}
%     %         Afctn = @(x) H.M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
%     %         [dy,flag,relres,iter] = pcg(Afctn,rhs,tolCG,maxIterCG);
%     %     case {'CG-hyperElastic'}
%     %         M =@(x) H.d2D.P((H.d2D.dr'*H.d2D.d2psi*H.d2D.dr)*H.d2D.P(x));
%     %         Afctn = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m,H.d2S.yc);
%     %         [dy,flag,relres,iter] = pcg(Afctn,rhs,tolCG,maxIterCG);
%     %     otherwise
%     %         if isnumeric(H)
%     %             % if H is a matrix, solve the linear system using MATLAB's backslash
%     %             dy = H\rhs;
%     %         else
%     %             error(solver)
%     %         end


text:
function z = choose(n, m)
% n choose m = n!/m!/(n-m)!
z = prod([n-m+1:n]) / factorial(m);

author: zouchuhang

repo: LayoutNet


text:
function model = vargplvmUpdateStats(model, X_u)
% VARGPLVMUPDATESTATS Update stats for VARGPLVM model.
%
% COPYRIGHT : Michalis K. Titsias, 2009-2011
% COPYRIGHT : Neil D. Lawrence, 2009-2011
% COPYRIGHT : Andreas C. Damianou, 2010-2011
%
% SEEALSO : vargplvmOptimise, vargplvmExpandParam

% VARGPLVM

jitter = 1e-6;
%model.jitter = 1e-6;

% %%% Precomputations for the KL term %%%
%%% Maybe we should add something like model.dynamics.X =
%%% model.dynamics.vardist.means (if visualisation requires that).
 if isfield(model, 'dynamics') && ~isempty(model.dynamics)
     % This actually only operates on model.dynamics. The only reason we
     % send the whole model is that we need to get model.vardist.means.
     model = vargplvmDynamicsUpdateStats(model);
 end

%%% Precomputations for (the likelihood term of) the bound %%%

model.K_uu = kernCompute(model.kern, X_u);

% Always add jitter (so that the inducing variables are "jitter" function variables)
% and the above value represents the minimum jitter value
% Putting jitter always ("if" in comments) is like having a second
% whiteVariance in the kernel which is constant.
%if (~isfield(model.kern, 'whiteVariance')) | model.kern.whiteVariance < jitter
   % There is no white noise term so add some jitter.
if ~strcmp(model.kern.type, 'rbfardjit')
       model.K_uu = model.K_uu ...
                         + sparseDiag(repmat(jitter, size(model.K_uu, 1), 1));
end
%end

model.Psi0 = kernVardistPsi0Compute(model.kern, model.vardist);
model.Psi1 = kernVardistPsi1Compute(model.kern, model.vardist, X_u);
[model.Psi2, AS] = kernVardistPsi2Compute(model.kern, model.vardist, X_u);

%K_uu_jit = model.K_uu + model.jitter*eye(model.k);
%model.Lm = chol(K_uu_jit, 'lower');

                                      % M is model.k
%model.Lm = chol(model.K_uu, 'lower');
model.Lm = jitChol(model.K_uu)';      % M x M: L_m (lower triangular)   ---- O(m^3)
model.invLm = model.Lm\eye(model.k);  % M x M: L_m^{-1}                 ---- O(m^3)
model.invLmT = model.invLm'; % L_m^{-T}
model.C = model.invLm * model.Psi2 * model.invLmT;
model.TrC = sum(diag(model.C)); % Tr(C)
% Matrix At replaces the matrix A of the old implementation; At is more stable
% since it has a much smaller condition number than A=sigma^2 K_uu + Psi2
model.At = (1/model.beta) * eye(size(model.C,1)) + model.C; % At = beta^{-1} I + C

model.Lat = jitChol(model.At)'; % UNCOMMENT THIS IF YOU ARE NOT DOING THE DEBUG BELOW

%try %%%% DEBUG1
% model.Lat = jitChol(model.At)';
%catch e %%% DEBUG
%    model_At = model.At;
%    save('modellAt', 'model_At');
%    e.throw %%% DEBUG
%end %%%% DEBUG

% DEBUG2
% warning('') %%% DEBUG
% model.Lat = jitChol(model.At)';
% tmpWarn = lastwarn; %%% DEBUG
% if length(tmpWarn) > 41 && strcmp(tmpWarn(1:42), 'Matrix is not positive definite in jitChol') %%% DEBUG
%     % do stuff...
%     %%%model.id
%     warning('') % This is a good place to put a stop for the debugger...
% end %%% DEBUG

model.invLat = model.Lat\eye(size(model.Lat,1));
model.invLatT = model.invLat';
model.logDetAt = 2*(sum(log(diag(model.Lat)))); % log |At|

model.P1 = model.invLat * model.invLm; % M x M

% First multiply the two last factors; so, the large N is only involved
% once in the calculations (P1: MxM, Psi1':MxN, Y: NxD)
model.P = model.P1 * (model.Psi1' * model.m);

% Needed for both, the bound's and the derivs. calculations.
model.TrPP = sum(sum(model.P .* model.P));

%%% Precomputations for the derivatives (of the likelihood term) of the bound %%%

%model.B = model.invLmT * model.invLatT * model.P; %next line is better
model.B = model.P1' * model.P;
model.invK_uu = model.invLmT * model.invLm;
Tb = (1/model.beta) * model.d * (model.P1' * model.P1);
	Tb = Tb + (model.B * model.B');
model.T1 = model.d * model.invK_uu - Tb;

model.X = model.vardist.means;

%{
%%%%-------------------------- TEMP (for analysing how ind. pts are optimised) !!!!!!!!!
d=dist2(model.X_u, model.vardist.means);
[C,I]=min(d,[],2);
dNorm=dist2(model.X_u, model.vardist.means)/norm(model.X_u)^2;
[Cnorm,Inorm]=min(dNorm,[],2);

if Inorm-I ~= 0
    fprintf(1, 'Something is wrong with the distances! (Update Stats)');
end

try
    load TEMPInducingMeansDist
catch exception
    TEMPInducingMeansDist=[];
end
try
    load TEMPInducingIndices
catch exception
    TEMPInducingIndices=[];
end
try
    load TEMPInducingMeansDistNorm
catch exception
    TEMPInducingMeansDistNorm=[];
end
TEMPInducingMeansDist = [TEMPInducingMeansDist sum(C)];
save 'TEMPInducingMeansDist.mat' 'TEMPInducingMeansDist';
%fprintf(1,'# Total X_u - vardistMeans = %d\n',sum(C));
TEMPInducingIndices = [TEMPInducingIndices I];
save 'TEMPInducingIndices.mat' 'TEMPInducingIndices';

TEMPInducingMeansDistNorm = [TEMPInducingMeansDistNorm sum(Cnorm)];
save 'TEMPInducingMeansDistNorm.mat' 'TEMPInducingMeansDistNorm';
%}

author: SheffieldML

repo: vargplvm


text:
function [coef, res, yhat] = ar_fit( y,Pmax,nosearch)
% y(t) = coef(1) + coef(2)*y(t-1) + coef(3)*y(t-2) ...
%
% y is a column vector with the time-series we are fitting
%
% Pmax is the max model order (not counting the constant)

if(nargin<3)
    nosearch=false;
end

n = length(y);
Pmax=max(min(Pmax,n-1),1);

Xf = lagmatrix(y, 1:Pmax);
Xb = lagmatrix(flipud(y), 1:Pmax);

X = [ones(2*n,1) [Xf; Xb]];

yy=[y; flipud(y)];
lstValid=~isnan(yy) & ~isnan(sum(X,2));

if(nosearch)
    [Q,R] = qr(X(lstValid,:),0); % note that the zero is very important for performance
    invR = pinv(R);
    coef = invR * Q'*yy(lstValid);
    res = yy(lstValid) - X(lstValid,:)*coef;
else
    [coef, res] = stepwise(X(lstValid,:), yy(lstValid));
end
res = res(1:n);
yhat = y - res;

end

author: BUNPC


text:
function sFileOut = out_fopen_spm(OutputFile, sFileIn, ChannelMat)
% OUT_FOPEN_SPM: Saves the header of a new empty SPM .mat/.dat file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2017-2022

% Initialize SPM12+CAT12
[isInstalled, errMsg] = bst_plugin('Install', 'spm12');
if ~isInstalled
    error(errMsg);
end

% Get the two output file names: .mat and .dat
[fPath, fBase, fExt] = bst_fileparts(OutputFile);
MatFile = bst_fullfile(fPath, [fBase, '.mat']);
DatFile = bst_fullfile(fPath, [fBase, '.dat']);

% Create .mat structure
D.type      = 'continuous';
D.Nsamples  = round((sFileIn.prop.times(2) - sFileIn.prop.times(1)) .* sFileIn.prop.sfreq) + 1;
D.Fsample   = sFileIn.prop.sfreq;
D.timeOnset = sFileIn.prop.times(1);
% Trials
D.trials.label  = 'Undefined';
D.trials.events = repmat(struct('type', [], 'time', [], 'value', [], 'offset', [], 'duration', []),0);
D.trials.onset  = sFileIn.prop.times(1);
D.trials.bad    = 0;
D.trials.tag    = [];
D.trials.repl   = 1;
% Events
for iEvt = 1:length(sFileIn.events)
    for iOcc = 1:size(sFileIn.events(iEvt).times,2)
        i = length(D.trials.events) + 1;
        D.trials.events(i).type = sFileIn.events(iEvt).label;
        if ~isempty(sFileIn.events(iEvt).notes) && ~isempty(sFileIn.events(iEvt).notes{iOcc})
            D.trials.events(i).type = [D.trials.events(i).type '-' sFileIn.events(iEvt).notes{iOcc}];
        end
        D.trials.events(i).time = sFileIn.events(iEvt).times(1,iOcc);
        D.trials.events(i).value = iEvt;
        D.trials.events(i).offset = 0;
        if (size(sFileIn.events(iEvt).times,1) == 2)
            D.trials.events(i).duration = sFileIn.events(iEvt).times(2,iOcc) - sFileIn.events(iEvt).times(1,iOcc);
        else
            D.trials.events(i).duration = [];
        end
    end
end
% Channels
for i = 1:length(ChannelMat.Channel)
    D.channels(i).bad      = (sFileIn.channelflag(i) == -1);
    D.channels(i).label    = ChannelMat.Channel(i).Name;
    D.channels(i).type     = ChannelMat.Channel(i).Type;
    D.channels(i).X_plot2D = [];
    D.channels(i).Y_plot2D = [];
    D.channels(i).units    = 'm';
end
% Data
D.data = file_array(DatFile, [length(ChannelMat.Channel), D.Nsamples], 'float32-le');
% File name
D.fname = [fBase, '.mat'];
D.path = fPath;

% Get sensor types
iEeg = channel_find(ChannelMat.Channel, 'EEG, SEEG, ECOG, NIRS');
iMeg = channel_find(ChannelMat.Channel, 'MEG, MEG REG');
if ~isempty(iMeg)
    error(['MEG sensors are currently not supported by this function.' 10 ...
           'Please contact us through the Brainstorm user forum to request this feature.']);
end
% If all the channels have other types: consider it's all EEG
if isempty(iEeg) && isempty(iMeg)
    iEeg = 1:length(ChannelMat.Channel);
end
% Sensors
for i = 1:length(iEeg)
    if ~isempty(ChannelMat.Channel(iEeg(i)).Loc) && ~all(ChannelMat.Channel(iEeg(i)).Loc(:) == 0)
        D.sensors.eeg.chanpos(i,:) = ChannelMat.Channel(iEeg(i)).Loc(:,1)';
        D.sensors.eeg.elecpos(i,:) = ChannelMat.Channel(iEeg(i)).Loc(:,1)';
    else
        D.sensors.eeg.chanpos(i,:) = [NaN NaN NaN];
        D.sensors.eeg.elecpos(i,:) = [NaN NaN NaN];
    end
    D.sensors.eeg.chantype{i} = lower(ChannelMat.Channel(iEeg(i)).Type);
    D.sensors.eeg.chanunit{i} = 'V';
    D.sensors.eeg.label{i}    = ChannelMat.Channel(iEeg(i)).Name;
    D.sensors.eeg.type        = 'ctf';
    D.sensors.eeg.unit        = 'm';
    D.sensors.eeg.balance.current = 'none';
end
D.sensors.chantype = 'eeg';
% Rest of the structure
D.fiducials    = repmat(struct(), 0);
D.transform.ID = 'time';
D.condlist     = {};
D.montage.M    = [];
D.montage.Mind = 0;
D.history      = repmat(struct(), 0);
D.other        = repmat(struct(), 0);

% Save file
save(MatFile, 'D');

% Create a new header structure
sFileOut = sFileIn;
sFileOut.filename  = MatFile;
sFileOut.condition = '';
sFileOut.format    = 'SPM-DAT';
sFileOut.byteorder = 'l';
sFileOut.comment   = fBase;
% Force the destination compensation level
sFileOut.prop.currCtfComp = sFileOut.prop.destCtfComp;
sFileOut.header.ctfcomp   = sFileOut.prop.destCtfComp;
% Save pointer to the SPM file in the header
sFileOut.header.file_array = D.data;
sFileOut.header.nchannels  = length(ChannelMat.Channel);

author: brainstorm-tools


text:
function [im, kspFull] = pocs( ksp, iter, watchProgress )
%Partial-Fourier Reconstruction with POCS
%
% [im, kspFull] = pocs( kspIn, iter, watchProgr )
%
% === Input ===
%
%   kspIn:      Reduced Cartesian MRI Data-Set
%               Any dimension may be reduced,
%               but only one reduction dim. is allowed due to Physics/Math.
%
%               Allowed shapes for kspIn are...
%                 ... Ny x Nx
%                 ... Nc x Ny x Nx
%                 ... Nc x Ny x Nx x Nz
%
%               With Nc == number of receive Channels / Coils.
%
%               kspIn can either be a zero-padded array, so the partial Fourier property is obvious.
%               Or kspIn can be the measured data only, then we try to find k-space centre automagically
%               and create a zero-padded array with the full size, first.
%               Errors are however more likely to occur in the latter case.
%
%
%   iter:       No. of iterations
%   (optional)  default: iter = 20
%               Try on your own if larger iter improves your results!
%
%   watchProgr: true/false; Whether the progress of the reconstruction should
%   (optional)  be monitored in an image window.
%               In 3D data, only the central partition will be shown.
%
%
% === Output ===
%
%   im:         Reconstructed Images (channels not combined)
%
%   kspFull:    Reconstructed full k-space data
%
%
%
% === About the code ===
%
%   (1) We find out whether input data is
%       a) already zero-filled or
%       b) the pure asymmetric dataset, only
%
%       If b) is true, we zero-fill the data ourselves, which means we have to
%       determine the dimension first, in which the partial Fourier reduction was done.
%       We therefor find the position of the max. intensity in k-space which should
%       be identical to k-space centre. If the k-space centre is different from the
%       centre of the matrix, we know the partial Fourier dimension.
%       We then enlarge the matrix to its desired full size and fill the new part
%       with zeros.
%       If a) was true, finding the partial Fourier dimension is easy:
%       It is the dimension with all the zeros. :-)
%
%   (2) We create one low resolution image per channel/coil:
%
%       We need a symmetrically sampled part around the central k-space. Think of a
%       small stripe of phase encoding lines in the central k-space.
%       We only use these symmetric data (setting the rest zero) to reconstruct
%       low-resolution images. In order to avoid Gibbs-Ringing, a Hamming-filter
%       with the width of the stripe is multiplied with the data.
%       Additionally, all the fully sampled dimensions get a Hamming filter, too,
%       since we increase SNR, reduce further Gibbs-ringing and do not lose much
%       resolution.
%
%   (3) The phase of the low-resolution images is saved
%
%       POCS uses the fact that k-space data of real objects (no imaginary part)
%       have a point symmetry:
%           S(-k)  =  S*(k)      with k = (kx, ky, kz)
%       Our MRI objects are always complex, but we assume that phase variations
%       are due to coil sensitivities and B0-inhomgeneities,
%       which are both slowly varying (no high res. required).
%       Small-scale phase pertubations will decrease the reconstruction quality.
%
%   (4) Reference phase is applied in image space
%
%       We...
%       ... transform our zero-filled data to image space (IFFT)
%       ... remove the phase --> abs(image)
%       ... set the phase of our reference phase map --> image .* exp(1i.*phase)
%       ... transform back to k-space (FFT)
%       ... re-insert the measured data (self-consistency!)
%       ... goto "We..."
%
%       Iterating through the above steps fills the missing k-space points
%       with reasonable values.
%       If the phase varies slowly and there is no aliasing, this works very well.
%
% Aliasing artifacts are very challenging for POCS.
% So try to prevent aliasing in the first place (sufficient Field of View).

% =========================================================================
% Original code by Martin Blaimer
% * changed by Uvo Hoelscher
% * changed by Michael Völker
%        -- auto-detect PF dimension
%        -- auto-find centre point/line/partition
%        -- accept zerofilled or "pure" data
%        -- for multichannel or plain 2D data (single-channel)
%        -- 2D and 3D
%        -- error handling
%        -- comments, comments, comments
%        -- added option to monitor progress
%        -- moved code to seperate functions
%        -- smooth transition between acquired signal and
%           reconstructed data
%
% Problems? Suggestions?
%  --> michael.voelker@mr-bavaria.de
% =========================================================================

    % ( ===================================================================
    % Input Handling
    %
        if ~exist( 'ksp', 'var' ) || isempty(ksp) || ~isnumeric(ksp)
            error('pocs:input', 'First input must be Cartesian k-space data.')
        end
        if ~exist('iter','var') || isempty(iter) || numel(iter) ~= 1 || ~isnumeric(iter)
            iter = 20;
        end
        if ~exist('watchProgress','var')  || isempty(watchProgress) ||  numel(watchProgress) ~= 1 || ~isfinite(watchProgress)
            watchProgress = false;
        else
            watchProgress = logical( watchProgress );
        end

        Ndim = ndims( ksp );

        if Ndim > 4 || Ndim < 2
            error('pocs:shape','First input ''kspace'' should have one of these shapes:\n\n\t... Ny x Nx\n\t... Nc x Ny x Nx\n\t... Nc x Ny x Nx x Nz')
        end
        if Ndim == 2    % Ny x Nx
            ksp = reshape( ksp, [1 size(ksp)] );    %  1 x Ny x Nx  --> now we have one channel...
            wasAddedCoilDim = true;
            Ndim = 3;
        else
            wasAddedCoilDim = false;
        end

        % read the properties of the data
        sz   = size( ksp );
        sz   = sz(2:end);           % the (k-)spatial size of the array (i.e. without channels)
        prec = class( ksp );        % single or double precision?
    % ) ===================================================================

    % First: Check the sampling pattern (which parts of input are actually data?)
    smplPtrn = reshape( sum(abs(ksp),1) ~= 0, sz);        % Ny x Nx x Nz

    % ( ===================================================================
    % If input data is not yet zero-filled, do it here
    %
        if nnz(smplPtrn) == numel(smplPtrn)     % only the sampled data were passed / |N|umber of |N|on |Z|ero elements

            [ ksp, pfDim, isUpper, isLower, Nsmp ] = zerofillPFdim( ksp, wasAddedCoilDim );

            sz = size( ksp );
            sz = sz(2:end);     % ignore channels
        else
            [ pfDim, isUpper, isLower, Nsmp ] = detectPFdim( smplPtrn, wasAddedCoilDim );
        end
        clear  smplPtrn
    % ) ===================================================================

    if numel(sz) < 3
        sz(3) = 1;
    end
    Ny = sz(1);
    Nx = sz(2);
    Nz = sz(3);

    % ( ===================================================================
    % Handle ugly problems.
    %
        if ~isUpper && ~isLower
            error('pocs:UnknownErrorFound', 'I thought we are partial Fourier, but things seem to make no sense... :-(')
        end
    % ) ===================================================================

    % =====================================================================
    %
    %        We can now be sure to operate with zero-padded data.
    %
    % =====================================================================

    % initialize a cell of subscripts
    subs = { ':', ':', ':', ':' };      % all channels / all Ny / all Nx / all Nz

    % If the first entries are zero-filled (instead of the trailing ones),
    % flip the entries so we can treat them as if we pf'ed the first half of kspace.
    if isLower
        subs{pfDim+1} = sz(pfDim):-1:1;     % ...esreveR
        ksp           = ksp(subs{:});       % !ecaps-k si sihT
        subs{pfDim+1} = 1:sz(pfDim);        % lalala, we didn't do anything...
    end

    % Find out which point is in the centre and which indices belong to the
    % symmetrically sampled part of k-space.
    [ centreLine, idxSym ] = findSymSampled( ksp, pfDim, Nsmp );

    szSym = numel( idxSym );                % 2 * (Nsmp - centreLine) + 1

    if isUpper
        fprintf('Using %g points around point %g\n', szSym,    centreLine    );
    else
        fprintf('Using %g points around point %g\n', szSym, sz(pfDim)-centreLine+1 );
    end

    % ( ===================================================================
    % build up a symmetric low-pass filter
    %
    filter = cast( 1, prec );
    for d = 1:Ndim-1

        reshRule = ones(1,Ndim);    % how the filter will be reshaped

        if d ~= pfDim   % Each standard dimension gets a simple low-pass filter

            filt1D = hamming( sz(d), 'periodic'  );

        else            % our partial Fourier dimension gets an extra nice filter

            % create a narrow filter and remove everything else
            filt1D          = zeros(sz(d), 1, prec);            % full-size filter
            tmp             = hann( szSym + 2, 'symmetric' );   % a very narrow window
            filt1D(idxSym)  = tmp(2:end-1);                     % cut out the zeros at the edges (we have data there!)

            % take a look:
            %figure, plot(filt1D)
        end

        % reshape the filter according to the dimension it represents
        reshRule(d+1) = sz(d);

        filt1D = reshape( filt1D, reshRule );
        filter = bsxfun( @times, filter, filt1D );      % iteratively build up a multidimensional filter
    end
    % ) ===================================================================

    % Apply the low-pass filter
    kspLowRes = bsxfun( @times, filter, ksp);
    clear  filt1D  filter  reshRule  idxSym

    % ( ===================================================================
    % prerequisites prior to the iteration loop
    %
    %  Set everything up here, do computations that you don't have
    %  to do in the loop, remove no longer needed variables...
    %
        % fftshift everything once before and after for-looping
        %  => less overhead during iteration
        ksp       = cmshiftnd(       ksp, [0  sz/2] );
        kspLowRes = cmshiftnd( kspLowRes, [0  sz/2] );

        % reorder arrays such that the fft-dimensions come first
        % => faster memory access
        ksp       = permute(       ksp, [2 3 4 1] );      % Ny x Nx x Nz x Nc
        kspLowRes = permute( kspLowRes, [2 3 4 1] );      %
        subs      = { subs{2}, subs{3}, subs{4}, subs{1} };

        % calc. initial image and the reference phase map
        im        =  fft( fft( fft( conj(ksp), [], 1), [], 2), [], 3);  % im's phase is wrong now, but we only want it's abs() to be correct
        phase     = ifft(ifft(ifft( kspLowRes, [], 1), [], 2), [], 3);
        phase     = exp(1i * angle(phase));

        % We use a trick in the loop to avoid using ifft (fft is faster).
        % We only need to calculate the factor 1/N ourselves, with N = prod(sz)
        phase = phase ./ prod(sz);      % 1/N is absorbed inside the phase array, once

        % create image with calculated phasemap from low res image
        im = abs(im) .* phase;

        % In the loop, we want to know where we have to copy the
        % measured data to, so we set the subscript of the pf dimension
        % accordingly.
        % We have to do this due to the ifftshift'ing above.
        tmp         = false( 1, sz(pfDim));
        tmp(1:Nsmp) = true;
        subs{pfDim} = find(ifftshift(tmp));

        % release RAM
        clear  tmp  kspLowRes

        % only keep the acquired data in memory
        ksp = ksp(subs{:});
    % ) ===================================================================

    % Helpers for pretty-printing:
    % Such a mess for such beautiful output!
    b = repmat('=',1,80);
    progress_str = 'starting POCS loop...';
    fprintf( '%s\n%s\n%s   %s', b, b(1), b(1), progress_str )
    edging = sprintf( '\n%s\n%s', b(1), b );
    fprintf( edging )

    % ( ===================================================================
    % iterative reconstruction POCS
    %
    tic
    for ii = double(~watchProgress) : iter

        if ii > 0

            % Fourier transform the image to k-space
            im = fft(fft(fft(  im  ,[],1),[],2),[],3);      % "im" is a really bad variable name now
                                                            % but we save a lot of RAM with this
            % Data Consistency:
            % insert original data where we have them
            im(subs{:}) = ksp;                              % "im" is still our reconstructed k-space signal

            % Fourier transform into image domain
            im = conj( im );
            im = fft(fft(fft(  im  ,[],1),[],2),[],3);      % Now, "im" is an image again.

            % create image with calculated phasemap from low res image
            im = abs(im) .* phase;

            prevLength = numel(progress_str) + numel(edging);
            t = toc;
            ETA = (t./ii) * iter  - t;
            progress_str = sprintf( 'Iteration %g/%g, in %g s,  ETA: %g s...', ii, iter, t, ETA );

            fprintf([repmat('\b',1,prevLength) '%s' '%s'], progress_str, edging );

        end % if ii > 0

        % a rough way to monitor the progress
        %
        if watchProgress
            tmp = ifftshift(sqrt(sum(abs(im(:,:,1,:).^2),4)));      % due to fftshift(), the 1st partition is the central one
            maxRange = sort( tmp(:), 'descend' );
            maxRange = maxRange( ceil(0.05 * numel(maxRange)) );    % ignore the "hottest" 5%

            if ~exist('pic','var')
                pic = [tmp tmp zeros(size(tmp),prec)];
                diffScale = 1;
            else
                delta = abs( pic(:,Nx+(1:Nx)) - tmp );
                diffScale = 0.5 * maxRange / median( delta(:) );
                pic(:,  Nx+(1:Nx)) = tmp;
                pic(:,2*Nx+(1:Nx)) = diffScale * delta;
                clear delta
            end

            figure(999)
            imagesc( pic, [0    maxRange ] )
            title(sprintf('\\bfiteration %g\ninitial     |    current     |     abs(previous - current) × %g', ii, diffScale ))
            axis image
            colormap(gray(256))
            drawnow
            clear tmp

            %if Nz == 1          % little pause for 2D (too fast otherwise)
            %    pause(2 / iter)
            %end
        end

    end % for ii = 1:iter
    fprintf([repmat('\b',1, numel(progress_str) + numel(edging)) 'POCS done! (%g s)' '%s\n\n'], t, edging );
    % ) ===================================================================

    clear  phase  pic

    % ( ===================================================================
    % The main part is over. Time for some thoughts.
    %
    % We began with a dataset that had fewer data samples than would be
    % necessary for an unambiguous image reconstruction. As a consequence,
    % an infinite number of images corresponds to the acquired data.
    % The above iteration picks that single image whose abs() fits the data
    % AND whose phase corresponds to the low-resolution phase, obtained
    % using the symmetric part of the data.
    %
    % Viewed in k-space, there is almost always a severe edge at the border
    % between acquired and interpolated data, which is due to imperfections
    % in the assumptions made.
    % Namely, phase often has some high frequency components which cannot be
    % accounted for in the low-resolution map. Additionally, there is noise
    % and we may have changing contrast or trajectory errors in our MRI
    % sequence.
    %
    %         ^
    %         | A A A A A A A A   \
    %         | A A A A A A A A
    %         | A A A A A A A A     acquired signal
    %      k2 | A A A A A A A A
    %         | A A A A A A A A   /
    %         | I I I I I I I I   \
    %         | I I I I I I I I     interpolated data
    %         | I I I I I I I I   /
    %         ----------------->
    %                 k1
    %
    % Empirically, it should be wise to create a smoother transition from
    % the acquired part of the signal to the interpolated data.
    %
        Ntrans = floor( (szSym-1)/3 );      % width of the transition zone

        % Create subscripts where we intend to keep the measured data, only.
        tmp                 = false( 1, sz(pfDim));
        tmp(1:Nsmp-Ntrans)  = true;
        subsPure            = subs;
        subsPure{pfDim}     = find(ifftshift(tmp));

        % Create subscripts where we want to have a smooth transition between
        % measured and phase-corrected data.
        subsTrans           = subs;
        subsTrans{pfDim}    = setdiff( subs{pfDim}, subsPure{pfDim} );

        % build a filter for the transition:
        tmp         = hann( 2*Ntrans+3, 'symmetric');
        filterTrans = tmp( Ntrans+3 : end-1 );
        filterTrans = reshape( filterTrans, [ ones(1,pfDim-1)  Ntrans  1] );

        % Seperate data in unfiltered part and transition zone.
        tmp = zeros( size(im), prec );
        tmp(subs{:}) = ksp;
        kspPure  = tmp(subsPure{:});
        kspTrans = tmp(subsTrans{:});
        clear  tmp  ksp

        im = fft(fft(fft(  im  ,[],1),[],2),[],3);      % "im" becomes k-space signal, again

        im(subsPure{:}) = kspPure;                      % strict data consistency for Nsmp-Ntrans samples
        im(subsTrans{:}) =  bsxfun( @times,   filterTrans,  kspTrans         )     ...
                          + bsxfun( @times, 1-filterTrans,  im(subsTrans{:}) );

        clear  subsPure  subsTrans  filterTrans  kspPure  kspTrans

        if nargout > 1
            kspFull = im;
        else
            kspFull = double.empty([sz 0]);     % kspFull exists, but no memory required
        end

        im = ifft(ifft(ifft(  im  ,[],1),[],2),[],3);   % "im" is an image, again
    % ) ===================================================================

    % ( ===================================================================
    % Undo the prerequisites (--> postrequisites???)
    %
        % undo the permutations
        im      = permute(      im, [4 1 2 3] );
        kspFull = permute( kspFull, [4 1 2 3] );
        subs    = { subs{4}, subs{1}, subs{2}, subs{3} };

        % undo the fftshifts
        im      = cmshiftnd(      im, [0  sz/2] );
        kspFull = cmshiftnd( kspFull, [0  sz/2] );

        % undo flipping
        if isLower
            subs{pfDim+1} = sz(pfDim):-1:1;
            im            = im(subs{:});
            kspFull       = kspFull(subs{:});
        end
    % ) ===================================================================

    if wasAddedCoilDim                              % we initially reshaped a simple 2D raw data matrix to be of size 1 x Ny x Nx
        im      = reshape(      im, Ny, Nx, [] );
        kspFull = reshape( kspFull, Ny, Nx, [] );
    end

end     % of pocs()

% =========================================================================
%                                                                         =
%                      SWAPPED CODE                                       =
%                                                                         =
% =========================================================================

function [ ksp, pfDim, isUpper, isLower, Nsmp ] = zerofillPFdim( ksp, wasAddedCoilDim )
    % Only the acquired data were passed and we have to find the asymmetric
    % dimension. Then we increase the size along this dimension and pad with 0.

    Ndim    = ndims( ksp ) - 1;     % one dimension was for the channels
    sz      = size(  ksp );
    sz      = sz(  2:end );         % ignore channel dimension
    Nc      = size(  ksp, 1 );
    prec    = class( ksp );

    % init some helper variables
    pfDim    = 0;               % partial Fourier reduction dimension
    isUpper  = false;
    isLower  = false;
    isPartialFourier = false(Ndim,1);

    % ( ===============================================================
    % autodetect the Partial Fourier dimension
    %
    for d = 1:Ndim

        centre = floor( sz(d)/2 ) + 1;

        tmp = squeeze( sum(abs(ksp),1) );
        for d2 = 1:Ndim
            if d2 ~= d
                tmp = max(tmp,[],d2);       % keep only the maximum of non-partial data points
            end
        end
        [ dummy, maxPos(d) ] = max( tmp(:) );   %#ok <-- don't use "~", for compatibility

        if abs(maxPos(d) - centre) >= 2      % significant asymmetry ==> partial Fourier acquisition

            isPartialFourier(d) = true;
            pfDim = d;
            Nsmp = sz(d);

            isUpper = maxPos(d) > centre;   % Did we sample the upper matrix part, so the lower part is missing...
            isLower = maxPos(d) < centre;   % ... or are the first data points missing (e.g. asymmetric echo)?
        end
    end % for d = 1:Ndim
    %
    % ) ===== (PF dim detection) ======================================

    switch nnz(isPartialFourier)    % |N|umber of |N|on |Z|ero elements
        case 0
            error( 'pocs:NoPfDim', 'No partial Fourier dimension found.' )
        case 1
            fprintf( 'Found partial Fourier along array dimension %d\n', pfDim + ~wasAddedCoilDim )
        otherwise
            error( 'pocs:TooManyPfDims', 'Partial Fourier only allowed in 1 dimension, but %g were found!', nnz(isPartialFourier) )
    end

    if pfDim == 0   % our init value above
        error('zerofillPF:NoPF','No partial Fourier property found!')
    end

    % initialize a cell of subscripts
    subs = { ':', ':', ':', ':' };      % all channels / all Ny / all Nx / all Nz

    c = maxPos(pfDim);

    if isUpper

        sz(pfDim) = 2 * (c - mod(c,2));         % determine the blown-up size we want to achieve
        subs{pfDim+1} = 1:Nsmp;

    elseif  isLower

        sz(pfDim) = 2 * (Nsmp - c + 1);
        c = floor( sz(pfDim)/2 ) + 1;
        sz(pfDim) = sz(pfDim) + 2*~mod(c,2);    % A hack for Stefan's data... keep an eye on this!
        subs{pfDim+1} = (1:Nsmp) + (sz(pfDim)-Nsmp);

    else
        error( 'zerofillPF:PFdimNotClassified', 'Could not tell how partial Fourier was implemented.' )
    end

    % do the zerofilling
    tmp = zeros( [Nc sz], prec );
    tmp(subs{:}) = ksp;
    ksp = tmp;

end % of zerofillPFdim()

function [ pfDim, isUpper, isLower, Nsmp ] = detectPFdim( smplPtrn, wasAddedCoilDim )
    % User passed already zero-padded data. This was nice, now it's easy
    % to find the partial Fourier dimension!

    Ndim = ndims( smplPtrn );
    sz = size( smplPtrn );

    % init some helper variables
    pfDim    = 0;               % partial Fourier reduction dimension
    isUpper  = false;
    isLower  = false;
    isPartialFourier = false( Ndim, 1 );

    % ( ===============================================================
    % Determine if this is a zerofilled partial Fourier measurement
    % and along which dimension the data is reduced.
    %
    % smplPtrn in Partial Fourier looks like this:
    %
    %     ^
    %     | 1 1 1 1 1 1 1 1     --->  sampling pattern is the same
    %     | 1 1 1 1 1 1 1 1           for all k1 points
    %     | 1 1 1 1 1 1 1 1
    %  k2 | 1 1 1 1 1 1 1 1             i.e. for programming:
    %     | 1 1 1 1 1 1 1 1             smplPtrn == repmat( smplPtrn(:,1,1), [1 Nx Nz] )
    %     | O O O O O 0 0 0
    %     | O O O O O 0 0 0
    %     | O O O O O 0 0 0
    %      ----------------->
    %           k1
    %
    for d = 1:Ndim

        subs = { ones(1,sz(d)),     ... % initialize a cell of subscripts we might be interested in
                 ones(1,sz(d)),     ...
                 ones(1,sz(d))  };
        subs{d} = 1:sz(d);              % we ask for all entries in the d'th dimension

        idx_d = sub2ind( sz, subs{:} ); % convert to linear array indices

        oneCol = smplPtrn( idx_d );     % one column of the d'th dimension

        % create a rule how to reshape oneCol
        reshRule = ones(1,Ndim);
        reshRule(d) = sz(d);                    % e.g. reshRule = [   1   1 128 ]
        oneCol = reshape( oneCol, reshRule);

        % create a rule how to replicate oneCol
        repRule = sz;
        repRule(d) = 1;                         % e.g. repRule  = [ 256 256   1 ]

        % Check if we get the sampling pattern again
        % just by replicating oneCol along the other dimensions
        isPartialFourier(d) = isequal( smplPtrn, repmat( oneCol, repRule ) );

        if isPartialFourier(d)
            pfDim = d;
            Nsmp  = nnz( oneCol );      % how many fully sampled lines do we have?

            % Sampled upper or lower part of k-space matrix?
            isUpper = isequal( oneCol(:).', [ true( 1,Nsmp)          false(1,sz(d)-Nsmp)    ]);
            isLower = isequal( oneCol(:).', [ false(1,sz(d)-Nsmp)    true( 1,Nsmp)          ]);
        end
    end
    % ) ===============================================================

    switch nnz(isPartialFourier)    % |N|umber of |N|on |Z|ero elements
        case 0
            error( 'pocs:NoPfDim', 'No partial Fourier dimension found.' )
        case 1
            fprintf( 'Found partial Fourier along array dimension %d\n', pfDim + ~wasAddedCoilDim )
        otherwise
            error( 'pocs:TooManyPfDims', 'Partial Fourier only allowed in 1 dimension!' )
    end

end % of detectPFdim()

function [ centreLine, idxSym ] = findSymSampled( ksp, pfDim, Nsmp )

    Ndim = ndims( ksp ) - 1;    % one for channels
    sz = size( ksp );
    sz = sz(2:end);

    % autodetect the central k-space line
    %if ~exist('centreLine', 'var') || isempty(centreLine)
        tmp = squeeze( sum(abs(ksp),1) );
        for d = 1:Ndim
           if d ~= pfDim
               tmp = max(tmp,[],d);     % keep only the maximum of non-partial data points
           end
        end
        [ dummy, centreLine] = max( tmp(:) );   %#ok the central line has the max intensity
    %end

    % calculate the size of the symmetric part and the full dataset
    startSym = centreLine - (Nsmp - centreLine);    % start of our symmetric sampling
    endSym   = centreLine + (Nsmp - centreLine);    % end of symmetric part
    idxSym   = startSym : endSym;

    if any(idxSym < 1)    ||   any(idxSym > sz(pfDim))
       error( 'pocs:BadDataProperty' , 'Symmetric part of k-space out of bounds.\nThe maximum k-space intensity is at index %g whereas it should be centred => near %g.\nThe way, zerofilling was done is probably wrong.\nCheck your input k-space.', centreLine, round(sz(pfDim)/2) )
    end

end % of findSymmetricSampled()

function x = cmshiftnd( x, shifts)
%Function to circularly shift N-D arrays

    if nargin < 2 || all(shifts(:) == 0)
       return                       % no shift
    end

    sz      = size( x );
    numDims = ndims(x);             % number of dimensions
    idx = cell(1, numDims);         % creates cell array of empty matrices,
                                    % one cell for each dimension

    for k = 1:numDims

        m = sz(k);
        p = ceil(shifts(k));

        if p < 0
            p = m + p;
        end

        idx{k} = [p+1:m  1:p];
    end

    % Use comma-separated list syntax for N-D indexing.
    x = x(idx{:});

end % of cmshiftnd()

% Avoid the need for the signal toolbox and implement
% hamming() and hann() manually:
%
function w = hamming( N, symFlag )
%Hamming window
%
% w = hamming(L) returns an L-point symmetric Hamming window in the column vector w.
% L should be a positive integer.
%
%  The coefficients of a Hamming window are computed from the following equation:
%
%       w(n) = 0.54  +  0.46 * cos(2*pi*n/N),   0 <= n <= N
%
%
% w = hamming( L, 'symFlag') returns an L-point Hamming window using the window sampling
% specified by 'symFlag', which can be either 'periodic'  or 'symmetric' (the default).
% The 'periodic' flag is useful for DFT/FFT purposes, such as in spectral analysis.
% The DFT/FFT contains an implicit periodic extension and the periodic flag enables a signal
% windowed with a periodic window to have perfect periodic extension.
% When 'periodic' is specified, hamming computes a length L+1 window and returns the first L points.
% When using windows for filter design, the 'symmetric' flag should be used.
%
% --> http://www.mathworks.de/de/help/signal/ref/hamming.html
% --> https://de.wikipedia.org/wiki/Hamming-Fenster

% implemented by Michael.Voelker@mr-bavaria.de, 2012

    if ~exist( 'N', 'var' ) || isempty(N) || numel(N) ~= 1 || ~isnumeric(N)  || ~isfinite(N) || N < 1 || floor(N) ~= N
        error( 'hamming:badSize', 'Window lenght must be a positive integer.' )
    end
    if ~exist( 'symFlag', 'var' ) || isempty(symFlag)
        symFlag = 'symmetric';
    end

    if N == 1
        w = 1;
        return
    end

    switch symFlag
        case 'symmetric'
            L = N-1;
        case 'periodic'
            L = N;
        otherwise
            error('hamming:symFlag', 'Unknown symmetry flag. Try ''symmetric'' (default) or ''periodic''.')
    end

    w = (0:N-1) - L/2;
    w = 0.54  +  0.46 * cos(2*pi * w(:)./L);

end % of hamming()

function w = hann( N, symFlag )
%von-Hann (Hanning) window
%
% w = hann(L) returns an L-point symmetric Hann window in the column vector w.
% L must be a positive integer.
%
% The coefficients of a Hann window are computed from the following equation:
%
%      w(n) = 0.5 * (1 + cos(2*pi*n/N)),   0 <= n <= N
%
% The window length is L = N+1.
%
% w = hann(L,'sflag') returns an L-point Hann window using the window sampling specified by 'sflag',
% which can be either 'periodic' or 'symmetric' (the default). The 'periodic' flag is useful for DFT/FFT purposes,
% such as in spectral analysis.
% The DFT/FFT contains an implicit periodic extension and the periodic flag enables a signal windowed
% with a periodic window to have perfect periodic extension.
% When 'periodic' is specified, hann computes a length L+1 window and returns the first L points.
% When using windows for filter design, the 'symmetric' flag should be used.
%
% --> http://www.mathworks.de/de/help/signal/ref/hann.html
% --> https://de.wikipedia.org/wiki/Hann-Fenster

% implemented by Michael.Voelker@mr-bavaria.de, 2012

    if ~exist( 'N', 'var' ) || isempty(N) || numel(N) ~= 1 || ~isnumeric(N)  || ~isfinite(N) || N < 1 || floor(N) ~= N
        error( 'hann:badSize', 'Window lenght must be a positive integer.' )
    end
    if ~exist( 'symFlag', 'var' ) || isempty(symFlag)
        symFlag = 'symmetric';
    end

    if N == 1
        w = 1;
        return
    end

    switch symFlag
        case 'symmetric'
            L = N-1;
        case 'periodic'
            L = N;
        otherwise
            error('hann:symFlag', 'Unknown symmetry flag. Try ''symmetric'' (default) or ''periodic''.')
    end

    w = (0:N-1) - L/2;
    w = 0.5 * ( 1 + cos(2*pi * w(:)./L) );

end % of hann()

author: Sable

repo: mcbench-benchmarks


text:
function [node,elem]=meshgrid6(varargin)
%
% [node,elem]=meshgrid6(v1,v2,v3,...)
%
% mesh an ND rectangular lattice by splitting
% each hypercube into 6 tetrahedra
%
% author: John D'Errico
% URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
% modified by Qianqian Fang (q.fang at neu.edu)
%
% input:
%    v1,v2,v3,... - numeric vectors defining the lattice in
%                   each dimension.
%                   Each vector must be of length >= 1
%
% output:
%    node - factorial lattice created from (v1,v2,v3,...)
%           Each row of this array is one node in the lattice
%    elem - integer array defining simplexes as references to
%           rows of "node".
%
% example:
%     [node,elem]=meshgrid6(0:5,0:6,0:4);
%     plotmesh(node,elem);
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

% dimension of the lattice
n = length(varargin);

% create a single n-d hypercube
% list of node of the cube itself
vhc=('1'==dec2bin(0:(2^n-1)));
% permutations of the integers 1:n
p=perms(1:n);
nt=factorial(n);
thc=zeros(nt,n+1);
for i=1:nt
  thc(i,:)=find(all(diff(vhc(:,p(i,:)),[],2)>=0,2))';
end

% build the complete lattice
nodecount = cellfun('length',varargin);
if any(nodecount<2)
  error 'Each dimension must be of size 2 or more.'
end
node = lattice(varargin{:});

% unrolled index into each hyper-rectangle in the lattice
ind = cell(1,n);
for i=1:n
ind{i} = 0:(nodecount(i)-2);
end
ind = lattice(ind{:});
k = cumprod([1,nodecount(1:(end-1))]);
ind = 1+ind*k';
nind = length(ind);

offset=vhc*k';
elem=zeros(nt*nind,n+1);
L=(1:nind)';
for i=1:nt
  elem(L,:)=repmat(ind,1,n+1)+repmat(offset(thc(i,:))',nind,1);
  L=L+nind;
end

if(n==2 || n==3)
  elem=meshreorient(node,elem);
end

% ======== subfunction ========
function g = lattice(varargin)
% generate a factorial lattice in n variables
n=nargin;
sizes = cellfun('length',varargin);
c=cell(1,n);
[c{1:n}]=ndgrid(varargin{:});
g=zeros(prod(sizes),n);
for i=1:n
g(:,i)=c{i}(:);
end

author: fangq


text:
function [kern, noise, ivmInfo] = ivmDeconstruct(model)

% IVMDECONSTRUCT break IVM in pieces for saving.
% FORMAT
% DESC takes an IVM model structure and breaks it into component
% parts for saving.
% ARG model : the model that needs to be saved.
% RETURN kern : the kernel component of the IVM model.
% RETURN noise : the noise component of the IVM model.
% RETURN ivmInfo : a structure containing the other information
% from the IVM: what the active set is, what the inactive set is
% and what the site parameters are.
%
% SEEALSO : ivmReconstruct
%
% COPYRIGHT : Neil D. Lawrence, 2005

% IVM

kern = rmfield(model.kern, 'Kstore');
kern = rmfield(kern, 'diagK');
noise = model.noise;
ivmInfo.I = model.I;
ivmInfo.J = model.J;
ivmInfo.m = model.m;
ivmInfo.beta = model.beta;
ivmInfo.type = 'ivm';

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function  [Erange,eta,Eshort,fs] = erange(CIJ)
%ERANGE     Shortcuts
%
%   [Erange,eta,Eshort,fs] = erange(CIJ);
%
%   Shorcuts are central edges which significantly reduce the
%   characteristic path length in the network.
%
%   Input:      CIJ,        binary directed connection matrix
%
%   Outputs:    Erange,     range for each edge, i.e. the length of the
%                           shortest path from i to j for edge c(i,j) AFTER
%                           the edge has been removed from the graph.
%               eta         average range for entire graph.
%               Eshort      entries are ones for shortcut edges.
%               fs          fraction of shortcuts in the graph.
%
%   Follows the treatment of 'shortcuts' by Duncan Watts
%
%
%   Olaf Sporns, Indiana University, 2002/2007/2008

N = size(CIJ,1);
K = length(nonzeros(CIJ));
Erange = zeros(N,N);
[i,j] = find(CIJ==1);

for c=1:length(i)
   CIJcut = CIJ;
   CIJcut(i(c),j(c)) = 0;
   [R,D] = reachdist(CIJcut);                   %#ok<ASGLU>
   Erange(i(c),j(c)) = D(i(c),j(c));
end;

% average range (ignore Inf)
eta = sum(Erange((Erange>0)&(Erange<Inf)))/length(Erange((Erange>0)&(Erange<Inf)));

% Original entries of D are ones, thus entries of Erange
% must be two or greater.
% If Erange(i,j) > 2, then the edge is a shortcut.
% 'fshort' is the fraction of shortcuts over the entire graph.

Eshort = Erange>2;
fs = length(nonzeros(Eshort))/K;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ECGtask_heartbeat_classifier < ECGtask

% ECGtask for ECGwrapper (for Matlab)
% ---------------------------------
%
% Description:
%
% Abstract class for defining ECGtask interface
%
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 18/2/2013
% Last update: 18/2/2013

    properties(GetAccess = public, Constant)
        name = 'ECG_heartbeat_classifier';
        target_units = 'ADCu';
        doPayload = true;
    end

    properties( GetAccess = public, SetAccess = private)
        % if user = memory;
        % memory_constant is the fraction respect to user.MaxPossibleArrayBytes
        % which determines the maximum input data size.
        memory_constant = 0.3;

        started = false;

    end

    properties( Access = private, Constant)

        fig_hdl = 1;
        cKnownModesOfOperation = {'auto' 'slightly-assisted', 'assisted'};

    end

    properties( Access = private )

    end

    properties

        progress_handle
        payload
        mode
        tmp_path

        % minimum amount of heartbeats required
        min_heartbeats_required = 500;
        % maximum amount of heartbeats to process each iteration
        max_heartbeats_per_iter = 4000;

    end

    methods

        function obj = ECGtask_heartbeat_classifier(obj)

            if( isempty(obj.mode) )
                obj.mode = obj.cKnownModesOfOperation{1};
            end

        end

        function Start(obj, ECG_header, ECG_annotations)

            if( ECG_header.nsig == 1 )
                disp_string_framed(2, 'The classifier was not developed for single lead signals.');
                return
            end

            lead_idx = get_ECG_idx_from_header(ECG_header);

            if( length(lead_idx) <= 1 )
                disp_string_framed(2, 'The classifier was not developed for single lead signals.');
                return
            end

            if(( isfield(ECG_annotations, 'time') && ~isempty(ECG_annotations.time)) )
                obj.started = true;
            else

                % payload property is used in this task to input an external QRS
                % detector, or manually corrected detections.
                aux_val = get_QRS_from_payload(obj.payload);

                if(~isempty(aux_val) && length(aux_val.time) > obj.min_heartbeats_required )
                    obj.started = true;
                else
                    warning('ECGtask_HB_classifier:BadPayload', 'Not enough heartbeats to perform classification. Try decreasing "min_heartbeats_required" or increasing the amount of heartbeats.\n');
                end
            end

        end

        function payload_out = Process(obj, ECG, ECG_start_offset, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx )

            payload_out = [];

            if( ~obj.started )
                obj.Start(ECG_header);
                if( ~obj.started )
                    cprintf('*[1,0.5,0]', 'Task %s unable to be started for %s.\n', obj.name, ECG_header.recname);
                    return
                end
            end

            % payload property is used in this task to input an external QRS
            % detector, or manually corrected detections.
            aux_QRS_loc = get_QRS_from_payload(obj.payload, ECG_start_offset, ECG_sample_start_end_idx);

            if(isempty(aux_QRS_loc))
                annotations_used = ECG_annotations;
            else
                annotations_used = aux_QRS_loc;
            end

%             aux_val = ECG_sample_start_end_idx + ECG_start_offset - 1;
%             disp_string_title(1, sprintf( 'Classifying from %d to %d', aux_val(1), aux_val(2) ) );

            [lead_idx, ECG_header ] = get_ECG_idx_from_header(ECG_header);

            [aux_val, ~, lablist ] = a2hbc('ECG', ECG(:,lead_idx), ...
                            'ECG_header', ECG_header, ...
                            'ECG_annotations', annotations_used, ...
                            'op_mode', obj.mode );

            bAux = annotations_used.time > ECG_sample_start_end_idx(1) & annotations_used.time < ECG_sample_start_end_idx(2);
            payload_out.anntyp = lablist(colvec(aux_val(bAux)),1);
            payload_out.time = annotations_used.time(bAux) + ECG_start_offset - 1;
            payload_out.lablist = lablist;

        end

        function payload = Finish(obj, payload, ECG_header)

        end

        function payload = Concatenate(obj, plA, plB)

            if( isempty(plA) )

                payload = plB;

            else
                for fn = rowvec(fieldnames(plA))

                    if( isfield(plA, fn{1}) && isfield(plB, fn{1}) )
                        payload.(fn{1}) = [ plA.(fn{1}); plB.(fn{1}) ];
                    elseif( ~isfield(plA, fn{1}) && isfield(plB, fn{1}) )
                        payload.(fn{1}) = plB.(fn{1});
                    end

                end
            end

        end

        %% property restriction functions

        function set.mode(obj,x)
            if( ischar(x) && ~isempty(strcmpi(x, obj.cKnownModesOfOperation)) )
                obj.mode = x;
            else
                warning('ECGtask_HB_classifier:BadArg', 'Invalid mode.');
            end
        end

    end

    methods ( Access = private )

    end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [s,xform] = relaxAlignAll(s, refImg, outMm, bb, interp)
%
% [s,xform] = relaxAlignAll(s, refImg, [outMm=refImg.pixdim], [bb=true], [interp=7])
%
% Aligns all the series in s to the NIFTI file (or struct) refImg.
% If refImg is empty, all series are aligned to the first series.
%
% If bb==true, then the bounding box will be taken from the reference
% image. If bb==false, it will be taken from the first image in the series.
% If bb is a valid bounding box, then that will be used.
%
%
% HISTORY:
% 2008.01.01 RFD: wrote it.
% 2008.10.06 RFD: fixed bug where s(1) was getting zeroed-out when refImg
% was empty.
% 2009.12.07 RFD: no longer set NANs to 0. NANs are useful for knowing what
% regions are outside the FOV. Also, changed it to use the actual transform
% from the dicom header.

if(~exist('refImg','var'))
  refImg = [];
end
if(ischar(s))
  s = dicomLoadAllSeries(s);
end
if(~exist('outMm','var'))
  outMm = [];
end
if(~exist('bb','var')||isempty(bb))
  bb = true;
end
if(~exist('interp','var')||isempty(interp))
  interp = 1;
end
if(numel(interp)<6)
   interp = [interp(1) interp(1) interp(1) 0 0 0];
end

% Don't try to process localizers...
locInds = ~cellfun('isempty',regexp({s.seriesDescription},'LOCALIZER'));
if(any(locInds))
    fprintf('Discarding series [%d] since they appear to be localizers.\n', find(locInds));
    s = s(~locInds);
end

nSeries = numel(s);

disp('Applying cannonical xform to all series...');
for(ii=1:nSeries)
    canXform = mrAnatComputeCannonicalXformFromDicomXform(s(ii).imToScanXform,size(s(ii).imData));
    [s(ii).imData,s(ii).mmPerVox] = applyCannonicalXform(s(ii).imData, canXform, s(ii).mmPerVox, false);
    % Now that the data have been flipped to our cannonical (axial)
    % orientation, we need to flip the imToScanXform to also reflect this
    % new orientation:
    s(ii).imToScanXform = inv(canXform*inv(s(ii).imToScanXform));
end

disp(['Aligning images to the first series "' s(1).seriesDescription '"...']);
s(1).xformToFirst = inv(s(1).imToScanXform);
Vfirst.uint8 = uint8(round(mrAnatHistogramClip(double(s(1).imData),0.3,0.99).*255));
Vfirst.mat = s(1).imToScanXform;

for(ii=2:nSeries)
    Vin.uint8 = uint8(round(mrAnatHistogramClip(double(s(ii).imData),0.3,0.99).*255));
    Vin.mat = s(ii).imToScanXform;
    evalc('transRot = spm_coreg(Vfirst, Vin);');
    fprintf('Series %d (%s) alignment: trans = [%0.4f %0.4f %0.4f], rot = [%0.4f %0.4f %0.4f].\n',ii,s(ii).seriesDescription,transRot);
    % Vin.mat\spm_matrix(transRot(end,:))*Vfirst.mat
    s(ii).xformToFirst = inv(Vin.mat)*spm_matrix(transRot(end,:));
end

if(~isempty(refImg))
    if(ischar(refImg))
        refImg = niftiRead(refImg);
    end
    disp(['Aligning first image to the template image in ' refImg.fname '...']);
    img = double(refImg.data);
    Vref.uint8 = uint8(round(mrAnatHistogramClip(img,0.4,0.98).*255));
    Vref.mat = refImg.qto_xyz;
    if(isempty(outMm))
        outMm = refImg.pixdim([1:3]);
    end

    if(numel(bb)~=6)
        if(bb)
            bb = mrAnatXformCoords(Vref.mat, [1 1 1; size(img)]);
        else
            bb = mrAnatXformCoords(Vref.mat, [1 1 1; size(s(1).imData)]);
        end
    end

    % Align to ref
    evalc('transRot = spm_coreg(Vref, Vfirst);');
    fprintf('Ref alignment: trans = [%0.4f %0.4f %0.4f], rot = [%0.4f %0.4f %0.4f].\n',transRot);
    % Vin.mat\spm_matrix(transRot(end,:))*Vfirst.mat
    disp('Resampling all series...');
    % Now resample
    for(ii=1:nSeries)
        fprintf('   resampling series %d...\n',ii);
        xf = s(ii).xformToFirst*spm_matrix(transRot(end,:));
        [s(ii).imData,s(ii).imToScanXform] = mrAnatResliceSpm(double(s(ii).imData), xf, bb, outMm, interp, 0);
    	s(ii).mmPerVox = outMm;
    end
else
    % No ref image- just resample everything to be aligned to the first
    % series and at outMm voxel size.
    if(numel(bb)~=6)
        bb = mrAnatXformCoords(Vfirst.mat, [1 1 1; size(s(1).imData)]);
    end
    if(isempty(outMm))
        outMm = s(1).mmPerVox;
    end
    disp('Resampling all series...');
    for(ii=1:nSeries)
        [s(ii).imData,s(ii).imToScanXform] = mrAnatResliceSpm(double(s(ii).imData), s(ii).xformToFirst, bb, outMm, interp, 0);
    	s(ii).mmPerVox = outMm;
    end
end
xform = s(1).imToScanXform;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function current_sensor()

# Opamp gain resistors
R = 68000.0; # 5400
r = 1000.0; # 1000

# Opamp input resistor network
r1 = 665.0; # 660
r2 = 562.0; # 560
r3 = 56000.0; # 4700

# Shunt resistance
shunt = 0.0025;

# Opamp gain
printf("Opamp gain: ");
G = (R + r)/r

# Relation between the input signal e and signal v calculation

A = -100.0:1.0:100.0;
shunt_v = shunt * A;

# k parameter
k = ( 1/r1 + 1/r2 + 1/r3);
U = 5.0;

# translation voltage
printf("Input translate: ");
trans = U/(k*r3)

printf("Oputput translate: ");
out_trans = trans * G

# opamp input
v = (trans + shunt_v/(k*r1));

# opamp output
S = v * G;

# maxi/min value
printf("Max: ");
max(S)
printf("Min: ");
min(S)

plot(A, shunt_v, "1", A, v, "2", A, S, "3")
legend('e', 'v', 'S');

# calculate output voltage for 1A
shunt_v_1 = shunt * 1;
v_1 = (trans + shunt_v_1/(k*r1));

printf("Output voltage at 1A");
S_1 = v_1 * G

author: open-bldc

repo: open-bldc-hardware


text:
function ld = sn_ld(l, n_sensor, n_source, n_snapshot)
%SN_LD Sufficient statistic for source number detection in MDL/AIC.
%   This function is used internally.
%Syntax:
%   ld = SN_LD(l, n_sensor, n_source, n_snapshot);
%Inputs:
%   l - Eigenvalues of the covariance matrix in descending order.
%   n_sensor - Number of sensors used. Should match the length of l.
%   n_source - Number of sources.
%   n_snapshot - Number of snapshots.
%Outputs:
%   ld - Computed value.
%Reference:
%   H. L. Van Trees, Optimum array processing. New York: Wiley, 2002.

diff = n_sensor - n_source;
l = l(n_source+1:end);
ld = n_snapshot * diff * log(sum(l)/diff/(prod(l)^(1/diff)));

end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function varargout = projPointOnPolyline(point, poly, varargin)
%PROJPOINTONPOLYLINE Compute position of a point projected on a polyline.
%
%   POS = projPointOnPolyline(POINT, POLYLINE)
%   Compute the position of the orthogonal projection of a point on a
%   polyline.
%   POINT is a 1-by-2 row vector containing point coordinates
%   POLYLINE is a N-by-2 array containing coordinates of polyline vertices
%   POS is the position of the point on the polyline, between 0 and the
%   number of vertices of the polyline. POS can be a non-integer value, in
%   this case, the integer part corresponds to the polyline edge index
%   (between 0 and Nv-1), and the floating-point part corresponds to the
%   relative position on i-th edge (between 0 and 1, 0: edge start, 1: edge
%   end).
%
%   When POINT is an array of points, returns a column vector with as many
%   rows as the number of points.
%
%   POS = projPointOnPolyline(POINT, POLYLINE, CLOSED)
%   Specifies if the polyline is closed or not. CLOSED can be one of:
%     'closed' -> the polyline is closed
%     'open' -> the polyline is open
%     a column vector of logical with the same number of elements as the
%       number of points -> specify individually if each polyline is
%       closed (true=closed).
%
%   [POS, DIST] = projPointOnPolyline(...)
%   Also returns the distance between POINT and POLYLINE.
%
%   Example
%     poly = [10 10; 20 10;20 20;10 20];
%     projPointOnPolyline([15 0], poly)
%     ans =
%         0.5000
%     projPointOnPolyline([0 16], poly)
%     ans =
%         3.0000
%
%   See also
%   points2d, polygons2d, polylinePoint, projPointOnPolygon
%   distancePointPolyline
%

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2009-04-30, using Matlab 7.7.0.471 (R2008b)
% Copyright 2009-2022 INRA - Cepia Software Platform

% check if input polyline is closed or not
closed = false;
if ~isempty(varargin)
    var = varargin{1};
    if strcmp('closed', var)
        closed = true;
    elseif strcmp('open', var)
        closed = false;
    elseif islogical(var)
        closed = var;
    end
end

% closes the polyline if necessary
if closed
    poly = [poly ; poly(1,:)];
end

% number of points
Np = size(point, 1);

% allocate memory results
pos     = zeros(Np, 1);
minDist = inf*ones(Np, 1);

% iterate on points
for p = 1:Np
    % build set of edges
    edges = [poly(1:end-1, :) poly(2:end, :)];

    % compute distance between current point and all edges
    [dist, edgePos] = distancePointEdge(point(p, :), edges);

    % update distance and position if necessary
    [minDist(p), edgeIndex] = min(dist);
    pos(p) = edgeIndex - 1 + edgePos(edgeIndex);
end

% process output arguments
if nargout <= 1
    varargout{1} = pos;
elseif nargout == 2
    varargout{1} = pos;
    varargout{2} = minDist;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [varargout] = maxDose(planC, varargin)
%Calculate a maxDose metric.
%Stand alone metric :       maxDose(planC, structNum, doseNum, 'Percent' or 'Absolute');
%Request m object   :       maxDose(planC, 'getnewmetric');
%Evaluate m object  :       maxDose(planC, m, 'evaluate');
%LM:  4 Jan 06, JOD: changed interface.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
optS = planC{indexS.CERROptions};

if(nargin > 1)
    call = varargin{end};
else
    warning('Incorrect Usage, try: maxDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');');
    return;
end

switch upper(call)
    case 'GETNEWMETRIC'
        m.name = 'maxDose';
        m.valueV = [];
        m.description = 'Returns the max dose in specified structure';
        m.functionName = @maxDose;
        m.note = '';
        m.params(1).name = 'Structure';
        m.params(1).type = 'DropDown';
        m.params(1).list = {planC{indexS.structures}.structureName};
        m.params(1).value = 1;
        m.params(2).name = 'IsTarget?';
        m.params(2).type = 'DropDown';
        m.params(2).list ={'Target', 'Not Target'};
        m.params(2).value = 2;
        m.units = [];
        m.range = [0 inf];
        m.doseSets = [1];
        varargout = {planC,m};
    return;

    case 'EVALUATE'
        m = varargin{1};
        structName = planC{indexS.structures}(m.params(1).value).structureName;
        m.note = structName;
        m.valueV = [];

        maxD = -inf;
		for i=1:length(m.doseSets)
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, m.params(1).value, m.doseSets(i));
            ans = calc_maxDose(doseBinsV, volsHistV, m.params(2).value);
            m.valueV = [m.valueV ans];
            doseArray = getDoseArray(planC{indexS.dose}(m.doseSets(i)));
            maxD = max(maxD, max(doseArray(:)));
        end
        if m.params(2).value == 1
            %Is target, high dose is best.
            m.range = [0 max(m.valueV)];
        elseif m.params(2).value == 2
            m.range = [maxD min(m.valueV)];
        end
        m.units = getDoseUnitsStr(i,planC);
        varargout = {planC,m};
        return;

    otherwise  %No specific call has been made, assume user just wants raw metric.
        if(nargin > 3)
            structNum = varargin{1};
            doseNum = varargin{2};
            volumeTypeString = varargin{3};
            if(strcmp(upper(volumeTypeString), 'ABSOLUTE'))
                volumeType = 2;
            else
                volumeType = 1;
            end
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, structNum, doseNum);
            varargout = {calc_maxDose(doseBinsV, volsHistV, volumeType)};
        else
            error('Not enough parameters, try: maxDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');')
            return;
        end
    end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function S=tria(A)
%%TRIA Square root matrix triangularization. Given a rectangular square
%      root matrix, obtain a lower-triangular square root matrix that is
%      square.
%
%INPUTS: A A numRowXnumCol matrix that is generally not square.
%
%OUTPUTS: S A lower-triangular matrix such that S*S'=A*A'. If
%           numCol>=numRow, then S is a square numRowXnumRow matrix.
%           Otherwise, S is a numRowXnumCol matrix.
%
%This is the tria function needed for various steps in the cubature Kalman
%filter and the square root Kalman filter. It is described in [1]. It has
%been slightly modified from the paper so that the diagonal elements remain
%positive.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%July 2012 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    [~,R]=qr(A',0);
    S=R';

    %Make the diagonal elements all positive.
    sel=diag(S)<0;
    S(:,sel)=-S(:,sel);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function v = ml_get_isolines(xtest,xtrain,alpha,kernel,kpar)
%ML_GET_ISOLINES: Plot the isolines of a kernel method.
%
%   input ---------------------------------------------------------------
%
%       o xtest       : (N x P), dataset of N samples of dimension P
%
%       o xtrain      : (M x P), original training data points
%
%       o alpha       : (M x 1), eigenvector of dimension M
%
%       o kernel      : string, the type of kernel
%                               [gauss,diag-gauss,poly,linear].
%
%       o kpar        : sig = kpar(1);  % variance for the gauss kernels.
%                       p   = kpar(1);	% polynome order
%                       c   = kpar(2);	% additive constant
%
%   output ----------------------------------------------------------------
%
%       o v           : value when projected onto the eigenvetors alpha
%
%   comment ---------------------------------------------------------------
%
%       Each data point of xtest is projected onto every alpha eigenvector.
%       Usually xtest will be (N x 2) since we are using this function for
%       plotting. The dim input is for the dimensions of the xtrain we want
%       to use in conjunction with xtest points. Basically both xtest
%       and xtrain should be of the same dimension.
%

N = size(xtest,1);
v = zeros(N,1);

if(strcmp(kernel,'poly') == true)
   kpar_ = [kpar(2) kpar(1)];
   kpar  = kpar_;
end

for i=1:N
    % Compute gramm matrix between test point i and all train points
    K_i  = km_kernel(xtest(i,:),xtrain,kernel,kpar);
    v(i) = alpha(:,1)' * K_i(:);

end

end

author: epfl-lasa


text:
function nodes2 = grAdjacentNodes(edges, node)
%GRADJACENTNODES Find list of nodes adjacent to a given node.
%
%   NEIGHS = grAdjacentNodes(EDGES, NODE)
%   EDGES: the complete edges list (containing indices of neighbor nodes)
%   NODE: index of the node
%   NEIGHS: the nodes adjacent to the given node.
%
%   NODE can also be a vector of node indices, in this case the result is
%   the set of neighbors of any input node, excluding the input nodes.
%
%   Example
%     % create a basic graph and display it
%     nodes = [10 10;20 10;10 20;20 20;27 15];
%     edges = [1 2;1 3;2 4;2 5;3 4;4 5];
%     figure; drawGraph(nodes, edges);
%     hold on; drawNodeLabels(nodes, 1:5)
%     axis equal; axis([0 40 0 30]);
%     % compute list of nodes adjacent to node with index 2
%     grAdjacentNodes(edges, 2)
%     ans =
%         1
%         4
%         5
%
%   See also
%     grAdjacentEdges

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-08-16
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

[i, j] = find(ismember(edges, node)); %#ok<ASGLU>
nodes2 = edges(i,1:2);
nodes2 = unique(nodes2(:));
nodes2 = sort(nodes2(~ismember(nodes2, node)));


text:
function [xs,ys] = snakeIterate(alpha,beta,gamma,x,y,NI,Fx,Fy)
%SNAKEITERATE Iterative solution of the snake equation.
%   [XS,YS] = SNAKEITERATE(ALPHA,BETA,GAMMA,X,Y,NI,Fx,Fy) computes the
%   [XS,YS] coordinates of a segmentation snake using the iterative
%   solution in Eq.(12-7) of DIPUM3E. Vectors X and Y are the initial
%   coordinates of the snake (provided in sequential order). These
%   vectors are updated during iteration. ALPHA, BETA, and GAMMA are
%   parameters in Eq. (12-7) and (12-8), and Fx, Fy are the 2D force
%   arrays obtained, for example, using DIPUM3E function snakeForce.
%
%   This function is normally run within an outer loop with snake-point
%   respacing after each execution of the loop. NI controls the number
%   of iterations of Eq. (12-7) before the snake points are respaced. A
%   common value of NI is 1, indicating one execution of point respacing
%   after each iteration of Eq. (12-7).
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% PRELIMINARIES.
K = numel(x);
% Multiply the forces by gamma.
Fx = gamma*Fx;
Fy = gamma*Fy;

% CONSTRUCT MATRIX A IN EQ. (12-8) FOR USE IN EQ. (12-7).
% First construct matrix D2 in Eq. (12-9).
a = -2*ones(K,1);
b = 1*ones(K-1,1);
D2 = diag(a) + diag(b,-1) + diag(b,1);
D2(1,K) = 1;
D2(K,1) = 1;
% Next construct D4 in Eq. (12-10).
a = 6*ones(K,1);
b = -4*ones(K-1,1);
c = 1*ones(K-2,1);
D4 = diag(a) + diag(b,-1) + diag(b,1) + diag(c,-2) + diag(c,2);
D4(1,K) = -4;
D4(K,1) = -4;
D4(1,K-1) = 1;
D4(K-1,1) = 1;
D4(2,K) = 1;
D4(K,2) = 1;
% Construct matrix A. The inverse operation is performed during
% iteration.
D = alpha*D2 - beta*D4;
A = eye(K) - D;

% ITERATIVE SOLUTION.
for I = 1:NI
	% Obtain the force vectors fx and fy in Eq. (12-7). These are
	% obtained from Fx and Fy, the 2-D array components of the force
	% field F, by interpolating values of Fx and Fy at the locations of
	% the current snake coordinates, x and y. The 0 in the following
	% function call avoids NaNs in areas where there are no points to
	% interpolate. The order y,x is because function interp2 works with
	% (col,row), as opposed to (row,col) as in the book.
	fx = interp2(Fx,y,x,'linear',0);
	fy = interp2(Fy,y,x,'linear',0);

	% Compute new values of x and y using Eq. (12-7). Note the use of x =
	% A\(x + fx), as opposed to x = inv(A)*(x + fx), and similarly for y.
	% The former method is faster and more accurate.
	x = A\(x + fx);
	y = A\(y + fy);
end

% FORM THE SNAKE.
xs = x;
ys = y;

author: dipum

repo: dipum-toolbox

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


text:
% child script of canlab_glm_publish
% runs scn_spm_design_check in publishable manner
% runs in directory containing subject level SPM analyses

z = '_________________________________________________';
diaryname = fullfile('model_review_work_log.txt');

if exist(fullfile(pwd,'SPM.mat'),'file')
    sublevs{1} = pwd;
else
    % figure out which directories are lower level analyses
    % criterion: they contain an SPM.mat file
    spmfiles = filenames(fullfile('*','SPM.mat'),'absolute');
    sublevs = cellfun(@fileparts,spmfiles,'UniformOutput',false);
end

for i = 1:numel(sublevs)
    diary(diaryname), fprintf('%s\n%s\n%s\n', z, sublevs{i}, z), diary off

    try
        scn_spm_design_check(sublevs{i}, 'events_only');
        snapnow
    catch exc
        diary(diaryname)
        disp(getReport(exc,'extended'))
%         fprintf('Design review FAILED for: %s.\n', sublevs{i})
        diary off
    end
end

close all

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% This is the training demo of FFDNet for denoising noisy grayscale images corrupted by
% AWGN.
%
% To run the code, you should install Matconvnet (http://www.vlfeat.org/matconvnet/) first.
%
% @article{zhang2018ffdnet,
%   title={FFDNet: Toward a Fast and Flexible Solution for {CNN} based Image Denoising},
%   author={Zhang, Kai and Zuo, Wangmeng and Zhang, Lei},
%   journal={IEEE Transactions on Image Processing},
%   volume={27},
%   number={9},
%   pages={4608-4622},
% }
%
% If you have any question, please feel free to contact with me.
% Kai Zhang (e-mail: cskaizhang@gmail.com)

%% ********************  Note  **********************************
% ** You should set the training images folders of "generatepatches.m" first. Then you can run "Demo_Train_FFDNet_gray.m" directly.
% **
% ** folders    = {'path_of_your_training_dataset'};% set this from "generatepatches.m" first!
% ** stride     = 10;                       % control the number of image patches, from "generatepatches.m"
% ** nimages  = round(length(filepaths));   % control the number of image patches, from "generatepatches.m"
% **
%% **************************************************************

format compact;
addpath('utilities');

%-------------------------------------------------------------------------
% Configuration
%-------------------------------------------------------------------------

opts.modelName        = 'FFDNet_gray'; % model name
opts.learningRate     = [logspace(-4,-4,100),logspace(-4,-4,100)/3,logspace(-4,-4,100)/(3^2),logspace(-4,-4,100)/(3^3),logspace(-4,-4,100)/(3^4)];% you can change the learning rate
opts.batchSize        = 64; % default
opts.gpus             = [1]; % this code can only support one GPU!
opts.numSubBatches    = 2;
opts.weightDecay      = 0.0005;
opts.expDir           = fullfile('data', opts.modelName);

%-------------------------------------------------------------------------
%  Initialize model
%-------------------------------------------------------------------------

net  = feval(['model_init_',opts.modelName]);

%-------------------------------------------------------------------------
%   Train
%-------------------------------------------------------------------------

[net, info] = model_train(net,  ...
    'expDir', opts.expDir, ...
    'learningRate',opts.learningRate, ...
    'numSubBatches',opts.numSubBatches, ...
    'weightDecay',opts.weightDecay, ...
    'batchSize', opts.batchSize, ...
    'modelname', opts.modelName, ...
    'gpus',opts.gpus) ;

author: cszn

repo: FFDNet

save_path: github-repos/MATLAB/cszn-FFDNet


text:
classdef ZeroLinePlotAdder < handle

    properties (Access = private)
        x
        y
        z
        figID
        meshBackground
        unfittedMesh
    end

    methods (Access = public)

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

        function addLine(obj)
            obj.createMeshBackground();
            obj.createUnfittedMesh();
            obj.plotLine();
        end
    end

    methods (Access = private)

        function init(obj,cParams)
            obj.x = cParams.x;
            obj.y = cParams.y;
            obj.z = cParams.z;
            obj.figID = cParams.figID;
        end

        function createMeshBackground(obj)
            [xp,yp] = meshgrid(obj.x,obj.y);
            xp = xp(:);
            yp = yp(:);
            tri = delaunay(xp,yp);
            s.coord = [xp yp];
            s.connec = tri;
            mesh = Mesh_Total(s);
            obj.meshBackground = mesh;
        end

        function createUnfittedMesh(obj)
            s.meshBackground = obj.meshBackground;
            s.backgroundMesh = obj.meshBackground.innerMeshOLD;
            s.boundaryMesh   = obj.meshBackground.boxFaceMeshes;
            s.unfittedType   = 'BOUNDARY';
            s.isInBoundary = false;
            cParams = SettingsMeshUnfitted(s);
            uMesh = UnfittedMesh(cParams);
            levelSet = obj.z;
            uMesh.compute(levelSet);
            obj.unfittedMesh = uMesh;
        end

        function plotLine(obj)
            figure(obj.figID)
            patch('vertices',obj.unfittedMesh.coord,...
                  'faces',obj.unfittedMesh.connec,...
                  'edgecolor',[0.5 0 0],...
                   'edgealpha',0.5,...
                   'edgelighting','flat',...
                   'facecolor',[1 0 0],...
                   'facelighting','flat',...
                   'LineWidth',6)
            axis('equal')
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [filename ok] = getROIfilename(vw, local)
% gets file names of ROIs to load
%
%     filename = getROIfilename(vw, [local])
%
% INPUTS:
%	vw: mrVista view.
%
%	local: flag indicating whether to use the local ROI directory (e.g.,
%	'Gray/ROIs/') or the shared ROI directory (e.g., '3DAnatomy/ROIs').
%	[default: 1, local]
%
% ras, 2009: now this calls a much-improved dialog in roiDialog.
if notDefined('vw'),    vw = getCurView; end
if notDefined('local')
    local = 1; %isequal(vw.viewType, 'Inplane');
else
    % if local is defined AND if local ~=1
    if ~local,
        if ispref('VISTA', 'verbose') && getpref('VISTA', 'verbose')==1
            disp('Loading from shared ROI directory');
        end
    end
end

switch vw.viewType
    case 'Inplane'
        initDir = 'Inplane/ROIs';
        dirList = {'Inplane/ROIs'};

    case 'Volume'
        sharedDir = fullfile(getAnatomyPath, 'ROIs');
        if local==1
            initDir = 'Volume/ROIs';
        else
            initDir = sharedDir;
        end
        dirList = {sharedDir 'Volume/ROIs' 'Gray/ROIs'};

    case 'Gray'
        if (ispref('VISTA','defaultROIPath'))

            sharedDir = fullfile(getAnatomyPath, getpref('VISTA','defaultROIPath'));
        else
            sharedDir = fullfile(getAnatomyPath, 'ROIs');
        end
        if local==1
            initDir = 'Gray/ROIs';
        else
            initDir = sharedDir;
        end
        dirList = {sharedDir 'Volume/ROIs' 'Gray/ROIs'};

    case 'Flat'
        initDir = fullfile(vw.subdir, 'ROIs');
        dirList = {initDir};

    otherwise
        initDir = '';
        dirList = {};
end

[filename ok] = roiDialog(initDir, dirList);
if ~ok, filename = ''; end

return