text:
%INITUNDISTORTRECTIFYMAP  Computes the undistortion and rectification transformation map
%
%     [map1, map2] = cv.initUndistortRectifyMap(cameraMatrix, distCoeffs, siz)
%     [...] = cv.initUndistortRectifyMap(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __cameraMatrix__ Input camera matrix `A = [f_x 0 c_x; 0 f_y c_y; 0 0 1]`
% * __distCoeffs__ Input vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12 or 14
%   elements. If the vector is empty, the zero distortion coefficients are
%   assumed.
% * __siz__ Undistorted image size `[w,h]`.
%
% ## Output
% * __map1__ The first output map. See `M1Type`.
% * __map2__ The second output map. See `M1Type`.
%
% ## Options
% * __R__ Optional rectification transformation in the object space (3x3
%   matrix). `R1` or `R2`, computed by cv.stereoRectify can be passed here. If
%   the matrix is empty, the identity transformation is assumed. default empty
% * __NewCameraMatrix__, __P__ New camera matrix (3x3)
%   `Ap = [fp_x 0 cp_x; 0 fp_y cp_y; 0 0 1]` or new projection matrix `P`
%   (3x4). If empty, uses the default new camera matrix from
%   cv.getDefaultNewCameraMatrix (with `CenterPrincipalPoint=true`).
%   default empty
% * __M1Type__ Type of the first output map, default -1 (equivalent to
%   `int16`). See cv.convertMaps. Accepted types are:
%   * __int16__ first output map is a MxNx2 `int16` array, second output map
%     MxNx1 `uint16` (fixed-point representation).
%   * __single1__ first output map is a MxNx1 `single` matrix, second output
%     map is MxNx1 `single` (separate floating-point representation).
%   * __single2__ first output map is a MxNx2 `single` matrix, second output
%     map is empty (combined floating-point representation).
%
% The function computes the joint undistortion and rectification
% transformation and represents the result in the form of maps for
% cv.remap. The undistorted image looks like original, as if it is captured
% with a camera using the camera matrix `NewCameraMatrix` and zero distortion.
% In case of a monocular camera, `NewCameraMatrix` is usually equal to
% `cameraMatrix`, or it can be computed by cv.getOptimalNewCameraMatrix for a
% better control over scaling. In case of a stereo camera, `NewCameraMatrix`
% is normally set to `P1` or `P2` computed by cv.stereoRectify.
%
% Also, this new camera is oriented differently in the coordinate space,
% according to `R`. That, for example, helps to align two heads of a stereo
% camera so that the epipolar lines on both images become horizontal and
% have the same y-coordinate (in case of a horizontally aligned stereo
% camera).
%
% The function actually builds the maps for the inverse mapping algorithm
% that is used by cv.remap. That is, for each pixel `(u,v)` in the
% destination (corrected and rectified) image, the function computes the
% corresponding coordinates in the source image (that is, in the original
% image from camera). The following process is applied:
%
%     x = (u - cp_x) / fp_x
%     y = (v - cp_y) / fp_y
%     [X Y Z]' = inv(R) * [x y 1]'
%     xp = X / W
%     yp = Y / W
%     r^2 = xp^2 + yp^2
%     xpp = xp*(1 + k1*r^2 + k2*r^4 + k3*r^6)/(1 + k4*r^2 + k5*r^4 + k6*r^6) +
%           2*p1*xp*yp + p2*(r^2 + 2*xp^2) + s1*r^2 + s2*r^4
%     ypp = yp*(1 + k1*r^2 + k2*r^4 + k3*r^6)/(1 + k4*r^2 + k5*r^4 + k6*r^6) +
%           p1*(r^2 + 2*yp^2) + 2*p2*xp*yp + s3*r^2 + s4*r^4
%       [xppp]   [R33(taux,tauy)              0 -R13(taux,tauy)]                  [xpp]
%     s*[yppp] = [             0 R33(taux,tauy) -R23(taux,tauy)] * R(taux,tauy) * [ypp]
%       [   1]   [             0              0               1]                  [  1]
%     map_x(u,v) = xppp * f_x + c_x
%     map_y(u,v) = yppp * f_y + c_y
%
% where `k1`, `k2`, `p1`, `p2`, `k3`, `k4`, `k5`, `k6`, `s1`, `s2`, `s3`,
% `s4`, `taux`, `tauy` are the distortion coefficients.
%
% In case of a stereo camera, this function is called twice: once for each
% camera head, after cv.stereoRectify, which in its turn is called after
% cv.stereoCalibrate. But if the stereo camera was not calibrated, it is still
% possible to compute the rectification transformations directly from the
% fundamental matrix using cv.stereoRectifyUncalibrated. For each camera, the
% function computes homography `H` as the rectification transformation in a
% pixel domain, not a rotation matrix `R` in 3D space. `R` can be computed
% from `H` as:
%
%     R = inv(cameraMatrix) * H * cameraMatrix
%
% where `cameraMatrix` can be chosen arbitrarily.
%
% See also: cv.remap, cv.convertMaps, cv.getOptimalNewCameraMatrix,
%  cv.stereoRectify, cv.stereoCalibrate, cv.stereoRectifyUncalibrated
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function f = polyVal2D(p,x,y,n,m)
%POLYVAL2D Evaluate a 2-D polynomial using Horner's method.
%   F = POLYVAL2D(P,X,Y) evaluates the 2-D polynomial P at the points
%   specified by X and Y, which must have the same dimensions. The output F
%   will have the same dimensions as X and  Y. N and M specify the order of X
%   and Y respectively. Polynomial coefficients are in the following order.
%
%   F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ...
%            P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ...
%            ...
%            P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)}
%
%
% See also: POLYFITN by John D'Errico on MathWorks MATLAB Central FEX
% http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn
%
% Mark Mikofski, http://poquitopicante.blogspot.com
% Version 1-0, 2013-03-13

%% LaTex
%
% $$f\left(x,y\right)=p_1 x^n y^m+p_2 x^{\left(n-1\right)} y^m+\ldots+p_{n+1} y^m+\ldots$$
%
% $$p_{n+2}x^ny^{\left(m-1\right)}+p_{n+3}x^{\left(n-1\right)}y^{\left(m-1\right)}+\ldots+p_{2\left(n+1\right)}y^{\left(m-1\right)}+\ldots$$
%
% $$\ldots$$
%
% $$p_{m\left(n+1\right)+1}*x^n+p_{m\left(n+1\right)+2}*x^{\left(n-1\right)}+\ldots+p_{\left(n+1\right)\left(m+1\right)}$$

%% check input args
validateattributes(p,{'numeric'},{'2d','nonempty','real','finite'}, ...
    'polyVal2D','p',1)
validateattributes(x,{'numeric'},{'nonempty','real','finite'}, ...
    'polyVal2D','x',2)
validateattributes(y,{'numeric'},{'nonempty','real','finite'}, ...
    'polyVal2D','y',3)
assert(all(size(x)==size(y)),'polyVal2D:sizeMismatch', ...
    'X and Y must be the same size.')
% use size of p to set n & m
pdims = size(p);
if nargin<4 && all(pdims>1)
    n = pdims(1)-1;
    m = pdims(2)-1;
end
validateattributes(n,{'numeric'},{'scalar','integer','positive','<',10}, ...
    'polyVal2D','n',4)
validateattributes(m,{'numeric'},{'scalar','integer','positive','<',10}, ...
    'polyVal2D','m',5)
if all(pdims>1)
    assert(pdims(1)==n+1,'polyVal2D:xOrderMismatch', ...
        'The number of x coefficients doesn''t match the order n.')
    assert(pdims(2)==m+1,'polyVal2D:yOrderMismatch', ...
        'The number of y coefficients doesn''t match the order m.')
end
%% evaluate polynomial P
f = p(1);
for ni = 1:n
    f = f.*x+p(1+ni);
end
for mi = 1:m
    mj = (n+1)*mi+1;
    g = p(mj);
    for ni = 1:n
        g = g.*x+p(mj+ni);
    end
    f = f.*y+g;
end

author: Sable

repo: mcbench-benchmarks


text:
clear all; close all; clc

mtime=21
L=40; n=512; x2=linspace(-L/2,L/2,n+1); x=x2(1:n); % spatial discretization
k=(2*pi/L)*[0:n/2-1 -n/2:-1].';                    % wavenumbers for FFT
t=linspace(0,2*pi,mtime);                      % time domain collection points

q=2*ones(mtime,1);
for j=2:2:mtime
  q(j)=4;
end
q(1)=1; q(mtime)=1;

N=2;
u=N*sech(x);
ut=fft(u);
[t,utsol]=ode45('ch_pod_sol_rhs',t,ut,[],k);
for j=1:length(t)
  usol(j,:)=ifft(utsol(j,:));
end

figure(1)
subplot(2,2,2)
waterfall(x,t,abs(usol)), shading interp, colormap([0 0 0])
%hold on
%pcolor(x,t,abs(usol2)), shading interp, %colormap(gray)
view(20,25)
set(gca,'Xlim',[-20 20],'Ylim',[0 2*pi],'Zlim',[0 4],'Xtick',[-20 -10 0 10 20],'Ytick',[0 3 6],'Ztick',[0 2 4], ...
    'Fontsize',[15]);
subplot(2,2,4)
waterfall(fftshift(k),t,abs(fftshift(utsol))), shading interp, colormap([0 0 0])
view(20,25)
set(gca,'Xlim',[-40 40],'Ylim',[0 2*pi],'Zlim',[0 80],'Xtick',[-40 -20 0 20 40],'Ytick',[0 3 6],'Ztick',[0 40 80], ...
    'Fontsize',[15]);

X=usol.';
for j=1:mtime
  X2(:,j)=X(:,j)*q(j);
end

[U,S,V]=svd(X);
[U2,S2,V2]=svd(X2);

%%

figure(2)
subplot(3,2,1)
semilogy(100*diag(S)/sum(diag(S)),'ko','Linewidth',[2]), grid on, hold on
semilogy(100*S(1,1)/sum(diag(S)),'bo','Linewidth',[2]),
semilogy(2,100*S(2,2)/sum(diag(S)),'go','Linewidth',[2]),
semilogy(3,100*S(3,3)/sum(diag(S)),'ro','Linewidth',[2]),

set(gca,'Xlim',[0 21],'Xtick',[1 6 11 16 21],'Ylim',[10^(-8) 10^3],'Ytick',[10^(-6) 10^(-3) 10^0 10^3],'Fontsize',[15])
xlabel('')

subplot(3,2,2)
semilogy(100*diag(S2)/sum(diag(S2)),'ko','Linewidth',[2]), grid on, hold on
semilogy(100*S2(1,1)/sum(diag(S2)),'bo','Linewidth',[2]),
semilogy(2,100*S2(2,2)/sum(diag(S2)),'go','Linewidth',[2]),
semilogy(3,100*S2(3,3)/sum(diag(S2)),'ro','Linewidth',[2]),

set(gca,'Xlim',[0 21],'Xtick',[1 6 11 16 21],'Ylim',[10^(-8) 10^3],'Ytick',[10^(-6) 10^(-3) 10^0 10^3],'Fontsize',[15])
xlabel('')

subplot(3,1,2)
sn=[-1 -1 1 -1 -1];
for j=1:5
  nrm=1;
  Up(:,j)=real(U(:,j))*sn(j)/nrm;
end
plot(x,real(Up(:,2)),'g','Linewidth',[2]), hold on
plot(x,real(Up(:,3)),'r','Linewidth',[2])
plot(x,real(Up(:,1)),'b','Linewidth',[2])
set(gca,'Xlim',[-10 10],'Xtick',[-10 -5 0 5 10],'Ylim',[-0.3 0.3],'Ytick',[-0.3 0 0.3],'Fontsize',[15])
subplot(3,1,3)
sn=[-1 -1 1 -1 -1];
for j=1:5
  nrm=sqrt(trapz(x,U2(:,j).^2));
  nrm=1;
  Up2(:,j)=real(U2(:,j))*sn(j)/nrm;
end
plot(x,real(Up2(:,3)),'r','Linewidth',[2]), hold on
plot(x,real(Up2(:,2)),'g','Linewidth',[2])
plot(x,real(Up2(:,1)),'b','Linewidth',[2])
set(gca,'Xlim',[-10 10],'Xtick',[-10 -5 0 5 10],'Ylim',[-0.3 0.3],'Ytick',[-0.3 0 0.3],'Fontsize',[15])
legend('mode 3','mode 2','mode 1')

figure(3)
plot(100*diag(S)/sum(diag(S)),'ro','Linewidth',[2])
hold on
plot(100*diag(S2)/sum(diag(S2)),'ko','Linewidth',[2])

figure(4)
sig1=100*diag(S)/sum(diag(S));
sig2=100*diag(S2)/sum(diag(S2));
diff=abs(sig1-sig2);
bar(diff)

%% low-rank reconstructions

X1=U(:,1:3)*S(1:3,1:3)*(V(:,1:3).');
X22=U2(:,1:3)*S2(1:3,1:3)*(V2(:,1:3).');

for j=1:mtime
  X222(:,j)=X22(:,j)/q(j);
end

figure(6)
subplot(2,2,1)
waterfall(x,t,abs(X1).'), shading interp, colormap([0 0 0])
subplot(2,2,2)
waterfall(x,t,abs(X222).'), shading interp, colormap([0 0 0])
E1=norm(X1-X);
E2=norm(X222-X);

author: dynamicslab


text:
function drawopt = drawtrackresult(drawopt, fno, frame, tmpl, param)
% function drawopt = drawtrackresult(drawopt, fno, frame, tmpl, param, pts)
%
%   drawopt : misc info for drawing, intitially []
%         [.showcoef] : shows coefficient
%         [.showcondens,thcondens] : show condensation candidates
%   fno : frame number
%   frame(fh,fw) : current frame
%   tmpl.mean(th,tw) : mean image
%       .basis(tN,nb) : basis
%   param.est : current estimate
%        .wimg : warped image
%       [.err,mask] : error, mask image
%       [.param,conf] : condensation
%
% uses: util/showimgs

% Copyright (C) 2005 Jongwoo Lim and David Ross.
% All rights reserved.

  figure(1); clf;
  set(gcf,'DoubleBuffer','on','MenuBar','none');
  colormap('gray');

sz = size(tmpl.mean);

axes('position', [0 0 1 1])
imagesc(frame, [0,1]); hold on;

text(5, 18, num2str(fno), 'Color','y', 'FontWeight','bold', 'FontSize',18);
drawbox(sz(1:2), param.est, 'Color','r', 'LineWidth',2.5);
axis off; hold off;

drawnow;

repo: FCNT


text:
function a = daub8_matrix ( n )

%*****************************************************************************80
%
%% DAUB8_MATRIX returns the DAUB8 matrix.
%
%  Discussion:
%
%    The DAUB8 matrix is the Daubechies wavelet transformation matrix
%    with 8 coefficients.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gilbert Strang, Truong Nguyen,
%    Wavelets and Filter Banks,
%    Wellesley-Cambridge Press, 1997,
%    ISBN: 0-9614088-7-1,
%    LC: TK7872.F5S79 / QA403.3.S87
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 8, and a multiple of 2.
%
%    Output, real A(N,N), the matrix.
%
  if ( n < 8 || mod ( n, 2 ) ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DAUB8_MATRIX - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 6 and a multiple of 2.\n' );
    error ( 'DAUB8_MATRIX - Fatal error!' );
  end

  a = zeros ( n, n );

  c = [ ...
    0.2303778133088964; ...
    0.7148465705529154; ...
    0.6308807679298587; ...
   -0.0279837694168599; ...
   -0.1870348117190931; ...
    0.0308413818355607; ...
    0.0328830116668852; ...
   -0.0105974017850690 ];

  for i = 1 : 2 : n - 1

    a(i,i)                  =   c(1);
    a(i,i+1)                =   c(2);
    a(i,i4_wrap(i+2,1,n))   =   c(3);
    a(i,i4_wrap(i+3,1,n))   =   c(4);
    a(i,i4_wrap(i+4,1,n))   =   c(5);
    a(i,i4_wrap(i+5,1,n))   =   c(6);
    a(i,i4_wrap(i+6,1,n))   =   c(7);
    a(i,i4_wrap(i+7,1,n))   =   c(8);

    a(i+1,i)                =   c(8);
    a(i+1,i+1)              = - c(7);
    a(i+1,i4_wrap(i+2,1,n)) =   c(6);
    a(i+1,i4_wrap(i+3,1,n)) = - c(5);
    a(i+1,i4_wrap(i+4,1,n)) =   c(4);
    a(i+1,i4_wrap(i+5,1,n)) = - c(3);
    a(i+1,i4_wrap(i+6,1,n)) =   c(2);
    a(i+1,i4_wrap(i+7,1,n)) = - c(1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% intlinear.m
%
%[A,B] = intlinear(x,y)
%
% Does de linear regression (least mean squares) for given (x,y).
% Returns A and B, with meaning y = A + B*x.

function    [A,B] = intlinear(x,y)

mx = mean(x);        my = mean(y);
mx2 = mean(x.^2);    my2 = mean(y.^2);
mxy = mean(x.*y);
A = (mx2*my-mx*mxy)/(mx2-mx^2);
B = (mxy - (mx*my))/(mx2-mx^2);


text:
function point = intersectEdges(edge1, edge2, varargin)
%INTERSECTEDGES Return all intersections between two set of edges.
%
%   P = intersectEdges(E1, E2);
%   returns the intersection point of edges E1 and E2.
%   E1 and E2 are 1-by-4 arrays, containing parametric representation of
%   each edge (in the form [x1 y1 x2 y2], see 'createEdge' for details).
%
%   In case of colinear edges, the result P contains [Inf Inf].
%   In case of parallel but not colinear edges, the result P contains
%   [NaN NaN].
%
%   If each input is N-by-4 array, the result is a N-by-2 array containing
%   the intersection of each couple of edges.
%   If one of the input has N rows and the other 1 row, the result is a
%   N-by-2 array.
%
%   P = intersectEdges(E1, E2, TOL);
%   Specifies a tolerance parameter to determine parallel and colinear
%   edges, and if a point belongs to an edge or not. The latter test is
%   performed on the relative position of the intersection point over the
%   edge, that should lie within [-TOL; 1+TOL].
%
%   See also
%   edges2d, intersectLines
%

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - Cepia Software Platform

% tolerance for precision
tol = 1e-14;

if nargin > 2
    tol = varargin{1};
end

%% Initialisations

% ensure input arrays are same size
N1  = size(edge1, 1);
N2  = size(edge2, 1);

% ensure input have same size
if N1 ~= N2
    if N1 == 1
        edge1 = repmat(edge1, [N2 1]);
        N1 = N2;
    elseif N2 == 1
        edge2 = repmat(edge2, [N1 1]);
    end
end

% initialize result array
x0  = zeros(N1, 1);
y0  = zeros(N1, 1);

%% Detect parallel and colinear cases

% indices of parallel edges
%par = abs(dx1.*dy2 - dx1.*dy2)<tol;
par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2));

% indices of colinear edges
% equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps
col = abs(  (edge2(:,1)-edge1(:,1)) .* (edge1(:,4)-edge1(:,2)) - ...
            (edge2(:,2)-edge1(:,2)) .* (edge1(:,3)-edge1(:,1)) ) < tol & par;

% Parallel edges have no intersection -> return [NaN NaN]
x0(par & ~col) = NaN;
y0(par & ~col) = NaN;

%% Process colinear edges

% colinear edges may have 0, 1 or infinite intersection
% Discrimnation based on position of edge2 vertices on edge1
if sum(col) > 0
    % array for storing results of colinear edges
    resCol = Inf * ones(size(col));

    % compute position of edge2 vertices wrt edge1
    t1 = edgePosition(edge2(col, 1:2), edge1(col, :));
    t2 = edgePosition(edge2(col, 3:4), edge1(col, :));

    % control location of vertices: we want t1<t2
    if t1 > t2
        tmp = t1;
        t1  = t2;
        t2  = tmp;
    end

    % edge totally before first vertex or totally after last vertex
    resCol(col(t2 < -tol))  = NaN;
    resCol(col(t1 > 1+tol)) = NaN;

    % set up result into point coordinate
    x0(col) = resCol(col);
    y0(col) = resCol(col);

    % touches on first point of first edge
    touch = col(abs(t2) < tol);
    x0(touch) = edge1(touch, 1);
    y0(touch) = edge1(touch, 2);

    % touches on second point of first edge
    touch = col(abs(t1-1) < tol);
    x0(touch) = edge1(touch, 3);
    y0(touch) = edge1(touch, 4);
end

%% Process non parallel cases

% process edges whose supporting lines intersect
i = find(~par);

% use a test to avoid process empty arrays
if sum(i) > 0
    % extract base parameters of supporting lines for non-parallel edges
    x1  = edge1(i,1);
    y1  = edge1(i,2);
    dx1 = edge1(i,3) - x1;
    dy1 = edge1(i,4) - y1;
    x2  = edge2(i,1);
    y2  = edge2(i,2);
    dx2 = edge2(i,3) - x2;
    dy2 = edge2(i,4) - y2;

    % compute intersection points of supporting lines
    delta = (dx2.*dy1 - dx1.*dy2);
    x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta;
    y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta;

    % compute position of intersection points on each edge
    % t1 is position on edge1, t2 is position on edge2
    t1  = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1);
    t2  = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2);

    % check position of points on edges.
    % it should be comprised between 0 and 1 for both t1 and t2.
    % if not, the edges do not intersect
    out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol;
    x0(i(out)) = NaN;
    y0(i(out)) = NaN;
end

%% format output arguments

point = [x0 y0];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function proj = projPointOnEllipse(point, elli)
%PROJPOINTONELLIPSE Project a point orthogonally onto an ellipse.
%
%   PROJ = projPointOnEllipse(PT, ELLI)
%   Computes the (orthogonal) projection of point PT onto the ellipse given
%   by ELLI.
%
%
%   Example
%     % create an ellipse and a point
%     elli = [50 50 40 20 30];
%     pt = [60 10];
%     % display reference figures
%     figure; hold on; drawEllipse(elli, 'b');
%     axis equal; axis([0 100 0 100]);
%     drawPoint(pt, 'bo');
%     % compute projected point
%     proj = projPointOnEllipse(pt, elli);
%     drawEdge([pt proj], 'b');
%     drawPoint(proj, 'k*');
%
%   See also
%     ellipses2d, distancePointEllipse, projPointOnLine, ellipsePoint
%     drawEllipse
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2022-07-17, using Matlab 9.12.0.1884302 (R2022a)
% Copyright 2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

% defaults
nMaxIters = 4;

% compute transform to centered axis-aligned ellipse
center = elli(1:2);
theta = deg2rad(elli(5));
transfo = createRotation(-theta) * createTranslation(-center);
pointT = transformPoint(point, transfo);

% retrieve ellipse semi axis lengths
a = elli(3);
b = elli(4);

% keep absolute values
px = abs(pointT(:,1));
py = abs(pointT(:,2));

% initial guess of solution
tx = 0.707;
ty = 0.707;

% iterate
for i = 1:nMaxIters
    x = a * tx;
    y = b * ty;

    ex = (a*a - b*b) * power(tx, 3) / a;
    ey = (b*b - a*a) * power(ty, 3) / b;

    rx = x - ex;
    ry = y - ey;

    qx = px - ex;
    qy = py - ey;

    r = hypot(ry, rx);
    q = hypot(qy, qx);

    tx = min(1, max(0, (qx .* r ./ q + ex) / a));
    ty = min(1, max(0, (qy .* r ./ q + ey) / b));
    t = hypot(ty, tx);
    tx = tx ./ t;
    ty = ty ./ t;

end

% computes coordinates of projection
projX = a * tx;
projY = b * ty;

% fix sign
projX(pointT(:,1) < 0) = -projX(pointT(:,1) < 0);
projY(pointT(:,2) < 0) = -projY(pointT(:,2) < 0);

% project pack to basis of original ellipse
proj = transformPoint([projX projY], inv(transfo));


text:
function [ n_data, x, fx ] = gud_values ( n_data )

%*****************************************************************************80
%
%% GUD_VALUES returns some values of the Gudermannian function.
%
%  Discussion:
%
%    The Gudermannian function relates the hyperbolic and trigonomentric
%    functions.  For any argument X, there is a corresponding value
%    GD so that
%
%      SINH(X) = TAN(GD).
%
%    This value GD is called the Gudermannian of X and symbolized
%    GD(X).  The inverse Gudermannian function is given as input a value
%    GD and computes the corresponding value X.
%
%    GD(X) = 2 * arctan ( exp ( X ) ) - PI / 2
%
%    In Mathematica, the function can be evaluated by:
%
%      2 * Atan[Exp[x]] - Pi/2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 13;

  fx_vec = [ ...
     -0.1301760336046015E+01, ...
     -0.8657694832396586E+00, ...
      0.0000000000000000E+00, ...
      0.9983374879348662E-01, ...
      0.1986798470079397E+00, ...
      0.4803810791337294E+00, ...
      0.8657694832396586E+00, ...
      0.1131728345250509E+01, ...
      0.1301760336046015E+01, ...
      0.1406993568936154E+01, ...
      0.1471304341117193E+01, ...
      0.1510419907545700E+01, ...
      0.1534169144334733E+01 ];

  x_vec = [ ...
     -2.0E+00, ...
     -1.0E+00, ...
      0.0E+00, ...
      0.1E+00, ...
      0.2E+00, ...
      0.5E+00, ...
      1.0E+00, ...
      1.5E+00, ...
      2.0E+00, ...
      2.5E+00, ...
      3.0E+00, ...
      3.5E+00, ...
      4.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function h = dmaxflat(N,d)

% returns 2-D diamond maxflat filters of order 'N'
% the filters are nonseparable and 'd' is the (0,0) coefficient, being 1 or 0 depending on use
% by Arthur L. da Cunha, University of Illinois Urbana-Champaign
% Aug 2004

if (N > 7 | N < 1), error('N must be in {1,2,3,4,5,6,7}');end;

switch N
 case 1
     h = [0 1 0; 1 0 1; 0 1 0]/4;
     h(2,2) = d;
 case 2
     h = [0 -1 0; -1 0 10; 0 10 0];
     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))]/32;
     h(3,3) = d;
 case 3
     h = [0  3   0  2;
          3  0  -27 0;
          0 -27  0 174;
          2  0  174  0];

     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))]/512;
     h(4,4) = d;
 case 4
     h = [0   -5     0   -3    0    ;
         -5    0    52    0    34   ;
          0    52    0   -276  0    ;
         -3    0   -276   0   1454  ;
          0    34    0   1454  0   ]/2^12;

     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))];
     h(5,5) = d;
 case 5
     h = [0    35    0    20    0    18;
          35   0   -425   0   -250    0;
          0   -425   0   2500   0    1610 ;
         20    0    2500  0  -10200   0;
          0   -250    0 -10200  0   47780;
         18    0    1610  0   47780   0]/2^17;

     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))];
     h(6,6) = d;
 case 6
     h = [0    -63    0    -35    0    -30    0   ;
         -63    0    882    0    495    0    444  ;
          0    882    0   -5910   0   -3420   0   ;
         -35    0   -5910   0   25875   0   16460 ;
          0    495    0    25875  0   -89730  0   ;
         -30    0   -3420   0  -89730   0   389112;
          0    44     0    16460  0   389112  0   ]/2^20;

     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))];
     h(7,7) = d;
 case 7
     h = [0    231     0     126     0      105      0      100   ;
         231    0    -3675    0    -2009     0     -1715     0    ;
          0   -3675    0    27930    0     15435     0     13804  ;
         126    0     27930   0   -136514    0    -77910     0    ;
          0   -2009    0   -136514   0     495145    0     311780 ;
         105    0     15435   0    495145    0    -1535709   0    ;
          0   -1715    0   -77910    0    -1535709   0    6305740 ;
         100    0    13804    0    311780    0    6305740    0    ]/2^24 ;
     h = [h fliplr(h(:,1:end-1))];
     h = [h ; flipud(h(1:end-1,:))];
     h(8,8) = d;
 end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function output = callvsdp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

% Convert from internal (sedumi-like) format to VSDPs sdpt3-like format
[blk,A,C,b]=sedumi2vsdp(F_struc(:,1),F_struc(:,2:end),c,K);

if options.savedebug
    ops = options.sdpt3;
    save sdpt3debug blk A C b ops -v6
end

% Solver to be used in VSDP
options.vsdp.model = interfacedata;

solvertime = tic;
[objt,Xt,yt,Zt,info] = mysdps_yalmip(blk,A',C,b,options);

% Compute rigorous lower bound (default behaviour)
if options.vsdp.verifiedlower
    [fL, Y, dL] = vsdplow_yalmip(blk,A',C,b,Xt,yt,Zt,[],options);
    if isnan(Y)
        info(1) = 11;
    end
    %[fL, Y, dL] = vsdplow(blk,A,C,b,Xt,yt,Zt)
else
    Y = [];
    fL = [];
    dL = [];
end

% Compute rigorous lower bound
if options.vsdp.verifiedupper
    %[fL, Y, dL] = vsdplow_yalmip(blk,A,C,b,[],[],[],[],options)
    % Now compute rigorous lower bound
    [fU, X, lb] = vsdpup_yalmip(blk,A',C,b,Xt,yt,Zt,[],options);
    %[fU, X, lb] = vsdpup(blk,A,C,b,Xt,yt,Zt);
else
    fU = [];
    X = [];
    lb = [];
end

solvertime = toc(solvertime);

Dual = [];
Slack = [];
top = 1;
if K.f>0
    Dual = [Dual;Xt{top}(:)];
    Slack = [Slack;Zt{top}(:)];
    top = top+1;
end
if K.l>0
    Dual = [Dual;[Xt{1:K.l}]'];
    Slack = [Slack;[Zt{1:K.l}]'];
    top = top + K.l;
end
if any(K.q)
    Dual = [Dual;Xt{top}(:)];
    Slack = [Slack;Zt{top}(:)];
    top = top + 1;
end
if any(K.s)
    for i = 1:length(K.s)
        Dual = [Dual;Xt{top+i-1}(:)];
        Slack = [Slack;Zt{top+i-1}(:)];
    end
end

if ~isempty(Y) & ~isnan(Y)
    Primal = Y;  % Primal variable in YALMIP
else
    Primal = yt;  % Primal variable in YALMIP
end

% Convert error code
switch info(1)
    case 0
        problem = 0; % No problems detected
    case {-1,-5}
        problem = 5; % Lack of progress
    case {-2,-3,-4,-7}
        problem = 4; % Numerical problems
    case -6
        problem = 3; % Maximum iterations exceeded
    case -10
        problem = 7; % YALMIP sent incorrect input to solver
    case 1
        problem = 2; % Dual feasibility
    case 2
        problem = 1; % Primal infeasibility
    case 11
        problem = 11;
    otherwise
        problem = -1; % Unknown error
end

% always save output
if options.savesolveroutput
    solveroutput.objt = objt;
    solveroutput.Xt = Xt;
    solveroutput.yt = yt;
    solveroutput.Zt = Zt;
    solveroutput.fL = fL;
    solveroutput.Y  = Y;
    solveroutput.dL  = dL;
    solveroutput.fU = fU;
    solveroutput.X  = X;
    solveroutput.lb  = lb;
    solveroutput.info = info;
 else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.blk = blk;
    solverinput.A   = A;
    solverinput.C   = C;
    solverinput.b   = b;
    solverinput.options   = options.sdpt3;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,Slack,problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

function [blk,A,C,b]=sedumi2vsdp(Cin,Ain,c,K);

C = {};
A = {};
b  = -c;
blk = {};
top = 1;
k = 1;

if any(K.f)
    C{k,1} = Cin(top:top+K.f-1);
%    A{k} = -Ain(top:top+K.f-1,:);
    for j = 1:length(c)
        A{j,k} = -Ain(top:top+K.f-1,j);
    end
    blk{k,1} = 'u';
    blk{k,2} = K.f;
    top = top + K.f;
    k = k + 1;
end

if any(K.l)
    K.s = [repmat(1,1,K.l) K.s];
    K.s(K.s == 0) = [];
%     C{k,1} = Cin(top:top+K.l-1);
% %    A{k} = -Ain(top:top+K.l-1,:);
%     for j = 1:length(c)
%         A{j,k} = -Ain(top:top+K.l-1,j);
%     end
%     blk{k,1} = 'l';
%     blk{k,2} = K.l;
%     top = top + K.l;
%     k = k + 1;
end

if any(K.s)
    for i = 1:length(K.s)
        C{k,1} = reshape(Cin(top:top+K.s(i)^2-1),K.s(i),K.s(i));
        for j = 1:length(c)
            A{j,k} = reshape(-Ain(top:top+K.s(i)^2-1,j),K.s(i),K.s(i));
        end
        blk{k,1} = 's';
        blk{k,2} = K.s(i);
        k = k + 1;
        top = top + K.s(i)^2;
    end
end

A = A';

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% FUNCTION solve_trace_norm_RMTL
%   trace norm projection
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Jianhui and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%

function L_hat = solve_trace_norm_RMTL(L, alpha)

[d1 d2] = size(L);

if (d1 > d2)

    [U S V] = svd(L, 0);

    thresholded_value = diag(S) - alpha / 2;

    diag_S = thresholded_value .* ( thresholded_value > 0 );

    L_hat = U * diag(diag_S) * V';

else

    new_L = L';

    [U S V] = svd(new_L, 0);

    thresholded_value = diag(S) - alpha / 2;

    diag_S = thresholded_value .* ( thresholded_value > 0 );

    L_hat = U * diag(diag_S) * V';

    L_hat = L_hat';

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
classdef CumPwrStorageStateOfChargeConstraint < AbstractConstraint
    %CumPwrStorageStateOfChargeConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = CumPwrStorageStateOfChargeConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, ~)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            [value, ~] = lvd_ElectricalPowerGlobalTasks(stateLogEntry, 'cumStorageSoC');

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                [valueStateComp, ~] = lvd_ElectricalPowerGlobalTasks(stateLogEntryStateComp, 'cumStorageSoC');
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Electrical Power Cumulative Storage State of Charge';
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'EC';
            lbLim = 0;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = CumPwrStorageStateOfChargeConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function polys = readPolygonSet(filename)
%READPOLYGONSET Read a set of simple polygons stored in a file
%
%   POLY = readPolygonSet(FILENAME);
%   Returns the polygon stored in the file FILENAME.
%   Polygons are assumed to be stored in text files, without headers, with
%   x and y coordinates packed in two separate lines:
%     X11 X12 X13 ... X1N
%     Y11 Y12 Y13 ... Y1N
%     X21 X22 X23 ... X2N
%     Y21 Y22 Y23 ... Y2N
%
%   Each polygon may have a different number of vertices. The result is a
%   cell array of polygon, each cell containing a N-by-2 array representing
%   the vertex coordinates.
%
%   See also:
%   polygons2d
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 11/04/2004.
%

% the set of polygons (no pre-allocation, as we do not know how many
% polygons are stored)
polys = {};

% index of polygon
p = 0;

% open file for reading
fid = fopen(filename, 'rt');

% use an infinite loop, terminated in case of EOF
while true
    % set of X, and Y coordinates
    line1 = fgetl(fid);
    line2 = fgetl(fid);

    % break loop if end of file is reached
    if line1 == -1
        break;
    end

    % create a new polygon by concatenating vertex coordinates
    p = p + 1;
    polys{p} = [str2num(line1)' str2num(line2)']; %#ok<AGROW,ST2NM>
end

% close file
fclose(fid);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ report ] = testReport( obj )
%testReport Analyze the sequence and return a text report
%   Currently no parameters are required. In future versions it may be
%   possible to (de)select some tests
%
% maxim.zaitsev@uniklinik-freiburg.de

% find the RF pulses and list flip angles
flipAnglesDeg=[];
for k=obj.rfLibrary.keys
    libData=obj.rfLibrary.data(k).array;
    if length(obj.rfLibrary.type)>=k
        rf=obj.rfFromLibData(libData,obj.rfLibrary.type(k));
    else
        rf=obj.rfFromLibData(libData);
    end
    flipAnglesDeg=[flipAnglesDeg abs(sum(rf.signal(1:end-1).*(rf.t(2:end)-rf.t(1:end-1))))*360]; %we use rfex.t(1) in place of opt.system.rfRasterTime
end
flipAnglesDeg=unique(flipAnglesDeg);

% calculate TE and TR

[duration, numBlocks, eventCount]=obj.duration();

%[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = obj.calculateKspace();
%[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = obj.calculateKspacePP();
[ktraj_adc, t_adc, ~, ~, t_excitation, ~] = obj.calculateKspacePP();

% trajectory calculation will fail for spin-echoes if seq is loaded from a
% file for the current file format revision (1.2.0) because we do not store
% the use of the RF pulses. Read function has an option 'detectRFuse' which
% may help...

%
kabs_adc=sum(ktraj_adc.^2,1).^0.5;
[kabs_echo, index_echo]=min(kabs_adc);
t_echo=t_adc(index_echo); % just a first estimate, see if we can improve it
if kabs_echo>eps
    i2check=[];
    % check if adc kspace trajectory has elements left and right to index_echo
    if index_echo > 1
        i2check=[i2check (index_echo-1)];
    end
    if index_echo < length(kabs_adc)
        i2check=[i2check (index_echo+1)];
    end
    for a=1:numel(i2check)
        v_i_to_0=-ktraj_adc(:,index_echo);
        v_i_to_t=ktraj_adc(:,i2check(a))-ktraj_adc(:,index_echo);
        % project v_i_to_0 to v_o_to_t
        p_vit=v_i_to_0'*v_i_to_t/(vecnorm(v_i_to_t)^2);
        if p_vit>0
            % we have forund a bracket for the echo and the proportionality
            % coefficient is p_vit
            t_echo=t_adc(index_echo)*(1-p_vit) + t_adc(i2check(a))*p_vit;
            break;
        end
    end
end

if ~isempty(t_excitation)
    t_ex_tmp=t_excitation(t_excitation<t_echo);
    TE=t_echo-t_ex_tmp(end);
    % TODO detect multiple TEs
else
    TE=NaN;
end

if (length(t_excitation)<2)
    TR=duration; % best estimate for now
else
    t_ex_tmp1=t_excitation(t_excitation>t_echo);
    if isempty(t_ex_tmp1)
        TR=t_ex_tmp(end)-t_ex_tmp(end-1);
    else
        TR=t_ex_tmp1(1)-t_ex_tmp(end);
    end
    % TODO check frequency offset to detect multiple slices
end

% check sequence dimensionality and spatial resolution
k_extent=max(abs(ktraj_adc),[],2);
k_scale=max(k_extent);
if (k_scale~=0)
    k_bins=4e6; % this defines our ability to separate k-space samples.
                % lower values give us imunity to rounding errors in k-space calculations
                % current code below (2nd pass) however merges neighboring cells (+-1)
    k_threshold=k_scale/k_bins;

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

    % bin the k-space trajectory to detect repetitions / slices
    k_len=size(ktraj_adc,2);
    k_repeat=zeros(1,k_len);
    k_storage=zeros(1,k_len);
    k_storage_next=1;
    % containers.Map only supports string as a key... (in older matlabs)
    %kmap = containers.Map('KeyType', 'char', 'ValueType', 'int32');
    kmap = java.util.HashMap;
    for i=1:k_len
        key_string = sprintf('%d ', int32(k_bins+round(ktraj_adc(:,i)/k_threshold)));
        % containers.Map does not have a proper find function so we use direct
        % access and catch the possible error
        k_storage_ind = kmap.get(key_string);
        if isempty(k_storage_ind)
            k_storage_ind=k_storage_next;
            kmap.put(key_string,k_storage_ind);
            k_storage_next=k_storage_next+1;
        end
        k_storage(k_storage_ind)=k_storage(k_storage_ind)+1;
        k_repeat(i) = k_storage(k_storage_ind);
    end
    % at this point k_storage(1:(k_storage_next-1)) is our visit frequency map
    Repeats_max=max(k_storage(1:(k_storage_next-1)));
    Repeats_min=min(k_storage(1:(k_storage_next-1)));
    Repeats_median=median(k_storage(1:(k_storage_next-1)));
    Repeats_unique=unique(k_storage(1:(k_storage_next-1)));
    Counts_unique=zeros(size(Repeats_unique));
    for i=1:numel(Repeats_unique)
        Counts_unique(i)=sum(Repeats_unique(i)==k_storage(1:(k_storage_next-1)));
    end

    ktraj_rep1=ktraj_adc(:,k_repeat==1);
    % TODO: think of something clever, e.g. detecting maximum delta-k
    % if length(k_extent)==2
    %     dt = delaunayTriangulation(ktraj_rep1');
    %     k = convexHull(dt);
    %     figure; plot(dt.Points(:,1),dt.Points(:,2), '.', 'markersize',10); hold on;
    %     plot(dt.Points(k,1),dt.Points(k,2), 'r'); hold off;
    %     %[V,R] = voronoiDiagram(dt);
    %     figure; voronoi(dt);
    % end
    % try to detect k-space lines or columns.
    k_counters=zeros(size(ktraj_rep1));
    dims=size(ktraj_rep1,1);
    %ordering=cell(1,dims);
    kmap = java.util.HashMap;
    for j=1:dims
        %kmap = containers.Map('KeyType', 'int32', 'ValueType', 'int32');
        k_storage=zeros(1,k_len);
        k_storage_next=1;
        kmap.clear();
        for i=1:size(ktraj_rep1,2)
            key=int32(round(ktraj_rep1(j,i)/k_threshold));
            k_storage_ind = kmap.get(key);
            if isempty(k_storage_ind) % attempt to account for rounding errors
                k_storage_ind = kmap.get(key+1);
            end
            if isempty(k_storage_ind) % attempt to account for rounding errors
                k_storage_ind = kmap.get(key-1);
            end
            if isempty(k_storage_ind)
                k_storage_ind=k_storage_next;
                kmap.put(key,k_storage_ind);
                k_storage_next=k_storage_next+1;
                k_storage(k_storage_ind)=ktraj_rep1(j,i);
                %fprintf('%d:%d(%g) ',k_storage_ind,key,ktraj_rep1(j,i));
            end
            %assert(k_storage_ind==k_storage(k_storage_ind));
            k_counters(j,i) = k_storage_ind;
        end
        %ordering{j}=cell2mat(kmap.values);
        %fprintf('\n');
    end
    unique_kpositions=max(k_counters,[],2);
    isCartesian=(prod(unique_kpositions)==size(ktraj_rep1,2));
else
    unique_kpositions=1;
end
% check gradient amplitudes and slew rates

% gradient waveform
gw_data=obj.waveforms_and_times(); % FIXME: avoid this second call for generating gradient shapes (1st one was inside of the k-space calculation routine)
gws=cell(size(gw_data));
ga=zeros(length(gw_data),1);
gs=zeros(length(gw_data),1);
% to calculate max absolute gradients and slew rates we have to play
% tricks... we namely have to interpolate the data to the common time axis
dim1ind = @(x, n) x(n,:);
common_time=unique(dim1ind([gw_data{:}],1));
gw_ct=zeros(length(gw_data),length(common_time));
gs_ct=zeros(length(gw_data),length(common_time)-1);
for gc=1:length(gw_data)
    if size(gw_data{gc},2)>0
        gws{gc}=(gw_data{gc}(2,2:end)-gw_data{gc}(2,1:end-1))./(gw_data{gc}(1,2:end)-gw_data{gc}(1,1:end-1)); % slew
        % interpolate to the common time
        gw_ct(gc,:)=interp1(gw_data{gc}(1,:),gw_data{gc}(2,:),common_time,'linear',0);
        gs_ct(gc,:)=(gw_ct(gc,2:end)-gw_ct(gc,1:end-1))./(common_time(2:end)-common_time(1:end-1));
        % max grad/slew per channel
        ga(gc)=max(abs(gw_data{gc}(2,:)));
        gs(gc)=max(abs(gws{gc}));
        % TODO: calculate grad RMS values (this is an interesting task in the piece-wise-linear domain)
    end
end

%figure; plot(common_time, gw_ct');
%figure; plot(common_time, sum(gw_ct.^2,1).^0.5);
%figure; plot(0.5*(common_time(1:end-1)+common_time(2:end)), sum(gs_ct.^2,1).^0.5);

% max absolute value grad/slew -- check for a worst case upon rotation
ga_abs=max(sum(gw_ct.^2,1).^0.5);
gs_abs=max(sum(gs_ct.^2,1).^0.5);

% check timing of blocks and delays (raster alignment)
[timing_ok, timing_error_report] = obj.checkTiming();

report = { sprintf('Number of blocks: %d\n',numBlocks),...
           [ sprintf('Number of events:\n'),...
             sprintf('   RF:   %6d\n',eventCount(2)),...
             sprintf('   Gx:   %6d\n',eventCount(3)),...
             sprintf('   Gy:   %6d\n',eventCount(4)),...
             sprintf('   Gz:   %6d\n',eventCount(5)),...
             sprintf('   ADC:  %6d\n',eventCount(6))],...
           [ sprintf('Sequence duration: %.6fs\n',duration),...
             sprintf('TE: %.6fs\n',TE),...
             sprintf('TR: %.6fs\n',TR) ],...
             sprintf('Flip angle: %.02f°\n', flipAnglesDeg),...
             sprintf('Unique k-space positions (a.k.a. columns, rows, etc): %d\n', unique_kpositions)};
if any(unique_kpositions>1)
    report = { report{:},...
           [ sprintf('Dimensions: %d\n', length(k_extent)),...
             sprintf('   Spatial resolution: %.02f mm\n', 0.5./k_extent*1e3) ],...
             sprintf('Repetitions/slices/contrasts: %.d  range: [%.d %.d]\n', Repeats_median, Repeats_min, Repeats_max) };
    report = { report{:},...
        sprintf('   %d k-space position(s) repeated %d times\n', [Counts_unique;Repeats_unique])};

    if isCartesian
       report = { report{:}, sprintf('Cartesian encoding trajectory detected\n') };
    else
       report = { report{:}, sprintf('Non-Cartesian/irregular encoding trajectory detected (e.g. EPI, spiral, radial, etc)\n') };
    end
end
if (timing_ok)
    report = { report{:}, sprintf('Block timing check passed successfully\n') };
else
    report = { report{:}, [ sprintf('Block timing check failed! Error listing follows:\n'),...
                            sprintf([timing_error_report{:}]) ] };
end
report = { report{:},...
    sprintf('Max. Gradient: %.0f Hz/m == %.02f mT/m\n', [ga mr.convert(ga,'Hz/m','mT/m')]'),...
    sprintf('Max. Slew Rate: %g Hz/m/s == %.02f T/m/s\n', [gs mr.convert(gs,'Hz/m/s','T/m/s')]') };
report = { report{:},...
    sprintf('Max. Absolute Gradient: %.0f Hz/m == %.02f mT/m\n', [ga_abs mr.convert(ga_abs,'Hz/m','mT/m')]'),...
    sprintf('Max. Absolute Slew Rate: %g Hz/m/s == %.02f T/m/s\n', [gs_abs mr.convert(gs_abs,'Hz/m/s','T/m/s')]') };

end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function H = cconvmtx2(h)

[block_size, num_blocks] = size(h);
num_elem = block_size*num_blocks;

H1 = spalloc(num_elem, block_size, block_size*nnz(h));
H = spalloc(num_elem, num_elem, num_elem*nnz(h));

% create the first n columns
for col = 1:block_size
    H1(:,col) = reshape(circshift(h, [col-1 0]), num_elem, 1);
end

% construct all blocks in H
for block = 1:num_blocks
    H(:,block_size*(block-1)+1:block_size*block) = circshift(H1, [(block-1)*block_size 0]);
end
end

author: flyers

repo: drone-tracking


text:
function simplex_monte_carlo_test02 ( )

%*****************************************************************************80
%
%% SIMPLEX_MONTE_CARLO_TEST02 estimates integrals in 6D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  e_test = [ ...
    0, 0, 0, 0, 0, 0; ...
    1, 0, 0, 0, 0, 0; ...
    0, 2, 0, 0, 0, 0; ...
    0, 2, 2, 0, 0, 0; ...
    0, 0, 0, 4, 0, 0; ...
    2, 0, 0, 0, 2, 2; ...
    0, 0, 0, 0, 0, 6 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Use SIMPLEX01_SAMPLE for a Monte Carlo estimate of an\n' );
  fprintf ( 1, '  integral over the interior of the unit simplex in 6D.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N' );
  fprintf ( 1, '        1      ' );
  fprintf ( 1, '        U      ' );
  fprintf ( 1, '         V^2   ' );
  fprintf ( 1, '         V^2W^2' );
  fprintf ( 1, '         X^4   ' );
  fprintf ( 1, '         Y^2Z^2' );
  fprintf ( 1, '         Z^6\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = simplex01_sample ( m, n, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:m) = e_test(1:m,j);

      value = monomial_value ( m, n, e, x );

      result = simplex01_volume ( m ) * sum ( value(1:n) ) / n;
      fprintf ( 1, '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Exact' );
  for j = 1 : 7

    e(1:m) = e_test(1:m,j);

    result = simplex01_monomial_integral ( m, e );
    fprintf ( 1, '  %14.6g', result );

  end

  fprintf ( 1, '\n' );

  return
end


text:
function spm_eeg_inv_vbecd_disp(action,varargin)
% Display the dipoles as obtained from VB-ECD
%
% FORMAT spm_eeg_inv_vbecd_disp('Init',D)
% Display the latest VB-ECD solution saved in the .inv{} field of the
% data structure D.
%
% FORMAT spm_eeg_inv_vbecd_disp('Init',D, ind)
% Display the ind^th .inv{} cell element, if it is actually a VB-ECD
% solution.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id$

% Note:
% unfortunately I cannot see how to ensure that when zooming in the image
% the dipole location stays in place...

global st

Fig     = spm_figure('GetWin','Graphics');
colors  = {'y','b','g','r','c','m'};              % 6 possible colors
marker  = {'o','x','+','*','s','d','v','p','h'};  % 9 possible markers
Ncolors = length(colors);
Nmarker = length(marker);

if nargin == 0, action = 'Init'; end;

switch lower(action),

%==========================================================================
case 'init'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('init',D,ind)
% Initialise the variables with GUI
%--------------------------------------------------------------------------

if nargin<2
    D = spm_eeg_load;
else
    D = varargin{1};
end
if nargin<3
    % find the latest inverse produced with vbecd
    Ninv = length(D.inv);
    lind = [];
    for ii=1:Ninv
        if isfield(D.inv{ii},'method') && ...
                strcmp(D.inv{ii}.method,'vbecd')
            lind = [lind ii];
        end
    end
    ind = max(lind);
    if ~ind,
        spm('alert*','No VB-ECD solution found with this data file!',...
                'VB-ECD display')
        return
    end
else
    ind = varargin{3};
end

% Stash dipole(s) information in sdip structure
sdip = D.inv{ind}.inverse;

% if the exit flag is not in the structure, assume everything went ok.
if ~isfield(sdip,'exitflag')
    sdip.exitflag = ones(1,sdip.n_seeds);
end

try
    Pimg = spm_vol(D.inv{ind}.mesh.sMRI);
catch
    Pimg = spm_vol(fullfile(spm('dir'), 'canonical', 'single_subj_T1.nii'));
end

spm_orthviews('Reset');
spm_orthviews('Image', Pimg, [0.0 0.45 1 0.55]);
spm_orthviews('MaxBB');
spm_orthviews('AddContext')
st.callback = 'spm_image(''shopos'');';
% remove clicking in image
for ii=1:3,
    set(st.vols{1}.ax{ii}.ax,'ButtonDownFcn',';');
end
WS = spm('WinScale');

% Build GUI
%==========================================================================
% Location:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[60 25 200 325].*WS, ...
    'DeleteFcn','spm_image(''reset'');');
uicontrol(Fig,'Style','Frame','Position',[70 250 180 90].*WS);
uicontrol(Fig,'Style','Text', 'Position',[75 320 170 016].*WS, ...
    'String','Current Position');
uicontrol(Fig,'Style','Text', 'Position',[75 295 35 020].*WS,'String','mm:');
uicontrol(Fig,'Style','Text', 'Position',[75 275 35 020].*WS,'String','vx:');
uicontrol(Fig,'Style','Text', 'Position',[75 255 75 020].*WS,'String','Img Intens.:');

st.mp = uicontrol(Fig,'Style','Text', 'Position',[110 295 135 020].*WS,'String','');
st.vp = uicontrol(Fig,'Style','Text', 'Position',[110 275 135 020].*WS,'String','');
st.in = uicontrol(Fig,'Style','Text', 'Position',[150 255  85 020].*WS,'String','');

c = 'if get(gco,''Value'')==1, spm_orthviews(''Xhairs'',''off''), else, spm_orthviews(''Xhairs'',''on''); end;';
uicontrol(Fig,'Style','togglebutton','Position',[95 220 125 20].*WS,...
    'String','Hide Crosshairs','Callback',c,'ToolTipString','show/hide crosshairs');

% Dipoles/seeds selection:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[300 25 180 325].*WS);
sdip.hdl.hcl = uicontrol(Fig,'Style','pushbutton','Position',[310 320 100 20].*WS, ...
        'String','Clear all','CallBack','spm_eeg_inv_vbecd_disp(''ClearAll'')');

sdip.hdl.hseed=zeros(sdip.n_seeds,1);
for ii=1:sdip.n_seeds
    if sdip.exitflag(ii)==1
        sdip.hdl.hseed(ii) = uicontrol(Fig,'Style','togglebutton','String',num2str(ii),...
            'Position',[310+rem(ii-1,8)*20 295-fix((ii-1)/8)*20 20 20].*WS,...
            'CallBack','spm_eeg_inv_vbecd_disp(''ChgSeed'')');
    else
        sdip.hdl.hseed(ii) = uicontrol(Fig,'Style','Text','String',num2str(ii), ...
            'Position',[310+rem(ii-1,8)*20 293-fix((ii-1)/8)*20 20 20].*WS) ;
    end
end

uicontrol(Fig,'Style','text','String','Select dipole # :', ...
        'Position',[310 255-fix((sdip.n_seeds-1)/8)*20 110 20].*WS);

txt_box = cell(sdip.n_dip,1);
for ii=1:sdip.n_dip, txt_box{ii} = num2str(ii); end
txt_box{sdip.n_dip+1} = 'all';
sdip.hdl.hdip = uicontrol(Fig,'Style','popup','String',txt_box, ...
        'Position',[420 258-fix((sdip.n_seeds-1)/8)*20 40 20].*WS, ...
        'Callback','spm_eeg_inv_vbecd_disp(''ChgDip'')');

% Dipoles orientation and strength:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[70 120 180 90].*WS);
uicontrol(Fig,'Style','Text', 'Position',[75 190 170 016].*WS, ...
                'String','Dipole orientation & strength');
uicontrol(Fig,'Style','Text', 'Position',[75 165 65 020].*WS, ...
                'String','x-y-z or.:');
uicontrol(Fig,'Style','Text', 'Position',[75 145 75 020].*WS, ...
                'String','theta-phi or.:');
uicontrol(Fig,'Style','Text', 'Position',[75 125 75 020].*WS, ...
                'String','Dip. intens.:');

sdip.hdl.hor1 = uicontrol(Fig,'Style','Text', 'Position', ...
                            [140 165  105 020].*WS,'String','a');
sdip.hdl.hor2 = uicontrol(Fig,'Style','Text', 'Position', ...
                            [150 145  85 020].*WS,'String','b');
sdip.hdl.int  = uicontrol(Fig,'Style','Text', 'Position', ...
                            [150 125  85 020].*WS,'String','c');

st.vols{1}.sdip = sdip;

% First plot = all the seeds that converged !
l_conv = find(sdip.exitflag==1);
if isempty(l_conv)
    error('No seed converged towards a stable solution, nothing to be displayed !')
else
    spm_eeg_inv_vbecd_disp('DrawDip',l_conv,1)
    set(sdip.hdl.hseed(l_conv),'Value',1); % toggle all buttons
end

%==========================================================================
case 'drawdip'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip,sdip)
% e.g. spm_eeg_inv_vbecd_disp('DrawDip',1,1,sdip)
% e.g. spm_eeg_inv_vbecd_disp('DrawDip',[1:5],1,sdip)
%--------------------------------------------------------------------------

if nargin < 2
    i_seed = 1;
else
    i_seed = varargin{1};
end
if nargin<3
    i_dip = 1;
else
    i_dip = varargin{2};
end

if nargin<4
    if isfield(st.vols{1},'sdip')
        sdip = st.vols{1}.sdip;
    else
        error('I can''t find sdip structure');
    end
else
    sdip = varargin{3};
    st.vols{1}.sdip = sdip;
end

if any(i_seed>sdip.n_seeds) || i_dip>(sdip.n_dip+1)
    error('Wrong i_seed or i_dip index in spm_eeg_inv_vbecd_disp');
end

% Note if i_dip==(sdip.n_dip+1) all dipoles are displayed simultaneously,
% The 3D cut will then be at the mean location of all sources !!!
if i_dip == (sdip.n_dip+1)
    i_dip = 1:sdip.n_dip;
end

% if seed indexes passed is wrong (no convergence) remove the wrong ones
i_seed(sdip.exitflag(i_seed)~=1) = [];
if isempty(i_seed)
    error('You passed the wrong seed indexes...')
end
if size(i_seed,2)==1, i_seed=i_seed'; end

% Display business
%--------------------------------------------------------------------------
loc_mm = sdip.loc{i_seed(1)}(:,i_dip);
if length(i_seed)>1
%     unit = ones(1,sdip.n_dip);
    for ii = i_seed(2:end)
        loc_mm = loc_mm + sdip.loc{ii}(:,i_dip);
    end
    loc_mm = loc_mm/length(i_seed);
end
if length(i_dip)>1
    loc_mm = mean(loc_mm,2);
end

% Place the underlying image at right cuts
spm_orthviews('Reposition',loc_mm);

if length(i_dip)>1
    tabl_seed_dip = [kron(ones(length(i_dip),1),i_seed') ...
                        kron(i_dip',ones(length(i_seed),1))];
else
    tabl_seed_dip = [i_seed' ones(length(i_seed),1)*i_dip];
end

% Scaling, according to all dipoles in the selected seed sets.
% The time displayed is the one corresponding to the maximum EEG power !
Mn_j = -1;
l3 = -2:0;
for ii = 1:length(i_seed)
    for jj = 1:sdip.n_dip
        Mn_j = max([Mn_j sqrt(sum(sdip.j{ii}(jj*3+l3,sdip.Mtb).^2))]);
    end
end
st.vols{1}.sdip.tabl_seed_dip = tabl_seed_dip;

% Display all dipoles, the 1st one + the ones close enough.
% Run through the 6 colors and 9 markers to differentiate the dipoles.
% NOTA: 2 dipoles coming from the same set will have same colour/marker
ind = 1 ;
dip_h = zeros(9,size(tabl_seed_dip,1),1);
    % each dipole displayed has 9 handles:
    %   3 per view (2*3): for the line, for the circle & for the error
js_m = zeros(3,1);

% Deal with case of multiple i_seed and i_dip displayed.
% make sure dipole from same i_seed have same colour but different marker.

pi_dip = find(diff(tabl_seed_dip(:,2)));
if isempty(pi_dip)
    % i.e. only one dip displayed per seed, use old fashion
    for ii=1:size(tabl_seed_dip,1)
        if ii>1
            if tabl_seed_dip(ii,1)~=tabl_seed_dip(ii-1,1)
                ind = ind+1;
            end
        end
        ic = mod(ind-1,Ncolors)+1;
        im = fix(ind/Ncolors)+1;

        loc_pl = sdip.loc{tabl_seed_dip(ii,1)}(:,tabl_seed_dip(ii,2));
        js = sdip.j{tabl_seed_dip(ii,1)}(tabl_seed_dip(ii,2)*3+l3,sdip.Mtb);
        vloc = sdip.cov_loc{tabl_seed_dip(ii,1)}(tabl_seed_dip(ii,2)*3+l3,tabl_seed_dip(ii,2)*3+l3);
        dip_h(:,ii) = add1dip(loc_pl,js/Mn_j*20,vloc, ...
                            marker{im},colors{ic},st.vols{1}.ax,Fig,st.bb);
        js_m = js_m+js;
    end
else
    for ii=1:pi_dip(1)
        if ii>1
            if tabl_seed_dip(ii,1)~=tabl_seed_dip(ii-1,1)
                ind = ind+1;
            end
        end
        ic = mod(ind-1,Ncolors)+1;
        for jj=1:sdip.n_dip
            im = mod(jj-1,Nmarker)+1;

            loc_pl = sdip.loc{tabl_seed_dip(ii,1)}(:,jj);
            js = sdip.j{tabl_seed_dip(ii,1)}(jj*3+l3,sdip.Mtb);
            vloc = sdip.cov_loc{tabl_seed_dip(ii,1)}(jj*3+l3,jj*3+l3);
            js_m = js_m+js;
            dip_h(:,ii+(jj-1)*pi_dip(1)) = ...
                add1dip(loc_pl,js/Mn_j*20,vloc, ...
                        marker{im},colors{ic},st.vols{1}.ax,Fig,st.bb);
        end
    end
end
st.vols{1}.sdip.ax = dip_h;

% Display dipoles orientation and strength
js_m = js_m/size(tabl_seed_dip,1);
[th,phi,Ijs_m] = cart2sph(js_m(1),js_m(2),js_m(3));
Njs_m = round(js_m'/Ijs_m*100)/100;
Angle = round([th phi]*1800/pi)/10;

set(sdip.hdl.hor1,'String',[num2str(Njs_m(1)),' ',num2str(Njs_m(2)), ...
        ' ',num2str(Njs_m(3))]);
set(sdip.hdl.hor2,'String',[num2str(Angle(1)),'� ',num2str(Angle(2)),'�']);
set(sdip.hdl.int,'String',Ijs_m);

% Change the colour of toggle button of dipoles actually displayed
for ii=tabl_seed_dip(:,1)
    set(sdip.hdl.hseed(ii),'BackgroundColor',[.7 1 .7]);
end

%==========================================================================
case 'clearall'
%==========================================================================
% Clears all dipoles, and reset the toggle buttons
%--------------------------------------------------------------------------

if isfield(st.vols{1},'sdip')
    sdip = st.vols{1}.sdip;
else
    error('I can''t find sdip structure');
end

disp('Clears all dipoles')
spm_eeg_inv_vbecd_disp('ClearDip');
for ii=1:st.vols{1}.sdip.n_seeds
    if sdip.exitflag(ii)==1
        set(st.vols{1}.sdip.hdl.hseed(ii),'Value',0);
    end
end
set(st.vols{1}.sdip.hdl.hdip,'Value',1);

%==========================================================================
case 'chgseed'
%==========================================================================
% Changes the seeds displayed
%--------------------------------------------------------------------------

% disp('Change seed')

sdip = st.vols{1}.sdip;
if isfield(sdip,'tabl_seed_dip')
    prev_seeds = p_seed(sdip.tabl_seed_dip);
else
    prev_seeds = [];
end

l_seed = zeros(sdip.n_seeds,1);
for ii=1:sdip.n_seeds
    if sdip.exitflag(ii)==1
        l_seed(ii) = get(sdip.hdl.hseed(ii),'Value');
    end
end
l_seed = find(l_seed);

% Modify the list of seeds displayed
if isempty(l_seed)
    % Nothing left displayed
    i_seed=[];
elseif isempty(prev_seeds)
    % Just one dipole added, nothing before
    i_seed=l_seed;
elseif length(prev_seeds)>length(l_seed)
    % One seed removed
    i_seed = prev_seeds;
    for ii=1:length(l_seed)
        p = find(prev_seeds==l_seed(ii));
        if ~isempty(p)
            prev_seeds(p) = [];
        end % prev_seeds is left with the index of the one removed
    end
    i_seed(i_seed==prev_seeds) = [];
    % Remove the dipole & change the button colour
    spm_eeg_inv_vbecd_disp('ClearDip',prev_seeds);
    set(sdip.hdl.hseed(prev_seeds),'BackgroundColor',[.7 .7 .7]);
else
    % One dipole added
    i_seed = prev_seeds;
    for ii=1:length(prev_seeds)
        p = find(prev_seeds(ii)==l_seed);
        if ~isempty(p)
            l_seed(p) = [];
        end % l_seed is left with the index of the one added
    end
    i_seed = [i_seed ; l_seed];
end

i_dip = get(sdip.hdl.hdip,'Value');
spm_eeg_inv_vbecd_disp('ClearDip');
if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
case 'chgdip'
%==========================================================================
% Changes the dipole index for the first seed displayed
%--------------------------------------------------------------------------

disp('Change dipole')
sdip = st.vols{1}.sdip;

i_dip = get(sdip.hdl.hdip,'Value');
if isfield(sdip,'tabl_seed_dip')
    i_seed = p_seed(sdip.tabl_seed_dip);
else
    i_seed = [];
end

if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('ClearDip')
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
case 'cleardip'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('ClearDip',seed_i)
% e.g. spm_eeg_inv_vbecd_disp('ClearDip')
%       clears all displayed dipoles
% e.g. spm_eeg_inv_vbecd_disp('ClearDip',1)
%       clears the first dipole displayed
%--------------------------------------------------------------------------

if nargin>2
    seed_i = varargin{1};
else
    seed_i = 0;
end
if isfield(st.vols{1},'sdip')
        sdip = st.vols{1}.sdip;
    else
        return; % I don't do anything, as I can't find sdip strucure
end

if isfield(sdip,'ax')
    Nax = size(sdip.ax,2);
    else
        return; % I don't do anything, as I can't find axes info
end

if seed_i==0 % removes everything
    for ii=1:Nax
        for jj=1:9
            delete(sdip.ax(jj,ii));
        end
    end
    for ii=sdip.tabl_seed_dip(:,1)
        set(sdip.hdl.hseed(ii),'BackgroundColor',[.7 .7 .7]);
    end
    sdip = rmfield(sdip,'tabl_seed_dip');
    sdip = rmfield(sdip,'ax');
elseif seed_i<=Nax % remove one seed only
    l_seed = find(sdip.tabl_seed_dip(:,1)==seed_i);
    for ii=l_seed
        for jj=1:9
            delete(sdip.ax(jj,ii));
        end
    end
    sdip.ax(:,l_seed) = [];
    sdip.tabl_seed_dip(l_seed,:) = [];
else
    error('Trying to clear unspecified dipole');
end
st.vols{1}.sdip = sdip;

%==========================================================================
case 'redrawdip'
%==========================================================================
% spm_eeg_inv_vbecd_disp('RedrawDip')
% redraw everything, useful when zooming into image
%--------------------------------------------------------------------------

% spm_eeg_inv_vbecd_disp('ClearDip')
% spm_eeg_inv_vbecd_disp('ChgDip')

% disp('Change dipole')
sdip = st.vols{1}.sdip;

i_dip = get(sdip.hdl.hdip,'Value');
if isfield(sdip,'tabl_seed_dip')
    i_seed = p_seed(sdip.tabl_seed_dip);
else
    i_seed = [];
end

if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('ClearDip')
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
otherwise
%==========================================================================
    warning('Unknown action string');
end
% warning(sw);
return

%==========================================================================
% dh = add1dip(loc,js,vloc,mark,col,ax,Fig,bb)
%==========================================================================
function dh = add1dip(loc,js,vloc,mark,col,ax,Fig,bb)
% Plots the dipoles on the 3 views, with an error ellipse for location
% Then returns the handle to the plots

global st
is = inv(st.Space);
loc = is(1:3,1:3)*loc(:) + is(1:3,4);
% taking into account the zooming/scaling only for the location
% NOT for the dipole's amplitude.
% Amplitude plotting is quite arbitrary anyway and up to some scaling
% defined for better viewing...

loc(1,:) = loc(1,:) - bb(1,1)+1;
loc(2,:) = loc(2,:) - bb(1,2)+1;
loc(3,:) = loc(3,:) - bb(1,3)+1;
% +1 added to be like John's orthview code

% prepare error ellipse
vloc = is(1:3,1:3)*vloc*is(1:3,1:3);
[V,E] = eig(vloc);
VE = V*diag(sqrt(diag(E))); % use std
% VE = V*E;   % or use variance ???

dh = zeros(9,1);
figure(Fig)
% Transverse slice, # 1
%----------------------
set(Fig,'CurrentAxes',ax{1}.ax)
set(ax{1}.ax,'NextPlot','add')
dh(1) = plot(loc(1),loc(2),[mark,col],'LineWidth',1);
dh(2) = plot(loc(1)+[0 js(1)],loc(2)+[0 js(2)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([1 2],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi)+loc(1);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(2);
dh(3) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{1}.ax,'NextPlot','replace')

% Coronal slice, # 2
%----------------------
set(Fig,'CurrentAxes',ax{2}.ax)
set(ax{2}.ax,'NextPlot','add')
dh(4) = plot(loc(1),loc(3),[mark,col],'LineWidth',1);
dh(5) = plot(loc(1)+[0 js(1)],loc(3)+[0 js(3)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([1 3],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi)+loc(1);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(3);
dh(6) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{2}.ax,'NextPlot','replace')

% Sagital slice, # 3
%----------------------
set(Fig,'CurrentAxes',ax{3}.ax)
set(ax{3}.ax,'NextPlot','add')
% dh(5) = plot(dim(2)-loc(2),loc(3),[mark,col],'LineWidth',2);
% dh(6) = plot(dim(2)-loc(2)+[0 -js(2)],loc(3)+[0 js(3)],col,'LineWidth',2);
dh(7) = plot(bb(2,2)-bb(1,2)-loc(2),loc(3),[mark,col],'LineWidth',1);
dh(8) = plot(bb(2,2)-bb(1,2)-loc(2)+[0 -js(2)],loc(3)+[0 js(3)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([2 3],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = -(e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi))+bb(2,2)-bb(1,2)-loc(2);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(3);
dh(9) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{3}.ax,'NextPlot','replace')

return

%==========================================================================
% pr_seed = p_seed(tabl_seed_dip)
%==========================================================================
function pr_seed = p_seed(tabl_seed_dip)
% Gets the list of seeds used in the previous display

ls = sort(tabl_seed_dip(:,1));
if length(ls)==1
    pr_seed = ls;
else
    pr_seed = ls([find(diff(ls)) ; length(ls)]);
end

%

% OLD STUFF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use it with arguments or not:
% - spm_eeg_inv_vbecd_disp('Init')
%        The routine asks for the dipoles file and image to display
% - spm_eeg_inv_vbecd_disp('Init',sdip)
%        The routine will use the avg152T1 canonical image
% - spm_eeg_inv_vbecd_disp('Init',sdip,P)
%        The routines dispays the dipoles on image P.
%
% If multiple seeds have been used, you can select the seeds to display
% by pressing their index.
% Given that the sources could have different locations, the slices
% displayed will be the 3D view at the *average* or *mean* locations of
% selected sources.
% If more than 1 dipole was fitted at a time, then selection of source 1
% to N is possible through the pull-down selector.
%
% The location of the source/cut is displayed in mm and voxel, as well as
% the underlying image intensity at that location.
% The cross hair position can be hidden by clicking on its button.
%
% Nota_1: If the cross hair is manually moved by clicking in the image or
%       changing its coordinates, the dipole displayed will NOT be at
%       the right displayed location. That's something that needs to be improved...
%
% Nota_2: Some seeds may have not converged within the limits fixed,
%       these dipoles are not displayed...
%
% Fields needed in sdip structure to plot on an image:
%       + n_seeds: nr of seeds set used, i.e. nr of solutions calculated
%       + n_dip: nr of fitted dipoles on the EEG time series
%       + loc: location of fitted dipoles, cell{1,n_seeds}(3 x n_dip)
%               remember that loc is fixed over the time window.
%       + j: sources amplitude over the time window,
%            cell{1,n_seeds}(3*n_dip x Ntimebins)
%       + Mtb: index of maximum power in EEG time series used

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % First point to consider
% loc_mm = sdip.loc{i_seed(1)}(:,i_dip);
%
% % PLace the underlying image at right cuts
% spm_orthviews('Reposition',loc_mm);
% % spm_orthviews('Reposition',loc_vx);
% % spm_orthviews('Xhairs','off')
%
% % if i_seed = set, Are there other dipoles close enough ?
% tabl_seed_dip=[i_seed(1) i_dip]; % table summarising which set & dip to use.
% if length(i_seed)>1
%   unit = ones(1,sdip.n_dip);
%   for ii = i_seed(2:end)'
%         d2 = sqrt(sum((sdip.loc{ii}-loc_mm*unit).^2));
%         l_cl = find(d2<=lim_cl);
%         if ~isempty(l_cl)
%             for jj=l_cl
%                 tabl_seed_dip = [tabl_seed_dip ; [ii jj]];
%             end
%         end
%   end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get(sdip.hdl.hseed(1),'Value')
% for ii=1:sdip.n_seeds, delete(hseed(ii)); end
% h1 = uicontrol(Fig,'Style','togglebutton','Position',[600 25 10 10].*WS)
% h2 = uicontrol(Fig,'Style','togglebutton','Position',[620 100 20 20].*WS,'String','1')
% h2 = uicontrol(Fig,'Style','checkbox','Position',[600 100 10 10].*WS)
% h3 = uicontrol(Fig,'Style','radiobutton','Position',[600 150 20 20].*WS)
% h4 = uicontrol(Fig,'Style','radiobutton','Position',[700 150 20 20].*WS)
% delete(h2),delete(h3),delete(h4),
% delete(hdip)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% test_dijkstra.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2012, 2015 University of Oxford
% Version: 0.3.3
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Basic Dijkstra

% distance matrix to describe the graph in
% http://en.wikipedia.org/wiki/Dijkstra's_algorithm
% (saved to matlab/test/Dijkstra_Animation_wikipedia.gif)

d = [...
    %  1   2   3   4   5   6
       0   7   9   0   0  14; ...
       7   0  10  15   0   0; ...
       9  10   0  11   0   2; ...
       0  15  11   0   6   0; ...
       0   0   0   6   0   9; ...
      14   0   2   0   9   0];

% run disjkstra's algorithm on first node
[dd, p] = dijkstra(sparse(d'), 1)

% dd' =
%
%      0     7     9    20    20    11
%
% p' =
%
%      0     1     1     3     6     3

% query several source nodes at the same time
[dd, p] = dijkstra(sparse(d'), [1 3 6])

% dd' =
%
%     0     7     9    20    20    11
%     9    10     0    11    11     2
%    11    12     2    13     9     0
%
% p' =
%
%      0     1     1     3     6     3
%      3     3     0     3     6     3
%      3     3     6     3     6     0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Dijkstra with limited number of targets

% query a source node, but we are only interested in the path to a single
% target node
[dd, p] = dijkstra(sparse(d'), 1, 3)

% dd' =
%
%     0     7     9    20   Inf    11
%
% p' =
%
%      0     1     1     3   NaN     3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Dijkstra with piggyback

% secondary distance matrix. The algorithm will ignore this distances in
% terms of searching for shortest-paths, but will create another distance
% output adding the values in D2 while following the path given by D

d2 = [...
    %  1   2   3   4   5   6
       0  10   7   0   0   6; ...
      10   0   2  11   0   0; ...
       7   2   0   9   0   9; ...
       0  11   9   0  14   0; ...
       0   0   0  14   0  15; ...
       6   0   9   0  15   0];

[dd, p, dd2] = dijkstra(sparse(d'), 1:6, 1:6, sparse(d2))

% dd2 =
%
%      0    10     7    16    31    16
%     10     0     2    11    25    11
%      7     2     0     9    24     9
%     16    11     9     0    14    18
%     31    26    24    14     0    15
%     16    11     9    18    15     0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Dijkstra with non-symmetric input matrix

% we start with a symmetric matrix
d = [...
    %  1   2   3   4   5   6
       0   7   9   0   0  14; ...
       7   0  10  15   0   0; ...
       9  10   0  11   0   2; ...
       0  15  11   0   6   0; ...
       0   0   0   6   0   9; ...
      14   0   2   0   9   0];

% the path from 1 to 5 is 1, 3, 6, 5, and has cost 20
[dd, p] = dijkstra(sparse(d'), 1);
dd(5)
graphpred2path(p', 5)

% in the opposite direction, same path, 5, 6, 3, 1
[dd, p] = dijkstra(sparse(d'), 5);
dd(1)
graphpred2path(p', 1)

% change the distance from 6 to 5 to be very expensive, so that the best
% path now is 1, 3, 4, 5, with cost 26
d(6, 5) = 20;
[dd, p] = dijkstra(sparse(d'), 1);
dd(5)
graphpred2path(p', 5)

% note that in the opposite direction nothing has changed, because the
% matrix is asymmetric
[dd, p] = dijkstra(sparse(d'), 5);
dd(1)
graphpred2path(p', 1)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%function calib_gui(mode)

% calib_gui(mode)
%
% Runs the Camera Calibration Toolbox.
% Set mode to 1 to run the memory efficient version.
% Any other value for mode will run the normal version (see documentation)
%
% INFORMATION ABOUT THE MEMORY EFFICIENT MODE FOR THE CAMERA CALIBRATION TOOLBOX:
%
% If your calibration images are large, or if you calibrate using a lot of images, you may have experienced memory problems
% in Matlab when using the calibration toolbox (OUT OF MEMORY errors). If this is the case, you can now run the
% new memory efficient version of the toolbox that loads every image one by one without storing them all in memory.
% If you choose to run the standard version of the toolbox now, you can always switch to the other memory efficient mode
% later in case the OUT OF MEMORY error message is encountered. The two modes of operation are totally compatible.

cell_list = {};

fig_number = 1;

title_figure = 'Camera Calibration Toolbox - Select mode of operation:';

cell_list{1,1} = {'Standard (all the images are stored in memory)','calib_gui_normal;'};
cell_list{2,1} = {'Memory efficient (the images are loaded one by one)','calib_gui_no_read;'};
cell_list{3,1} = {'Exit',['disp(''Bye. To run again, type calib_gui.''); close(' num2str(fig_number) ');']};

show_window(cell_list,fig_number,title_figure,290,18,0,'clean',12);

repo: ekfmonoslam


text:
function clenshaw_curtis_set_test ( )

%*****************************************************************************80
%
%% CLENSHAW_CURTIS_SET_TEST tests CLENSHAW_CURTIS_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CLENSHAW_CURTIS_SET_TEST\n' );
  fprintf ( 1, '  CLENSHAW_CURTIS_SET sets up a Clenshaw-Curtis rule;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the integral of sqrt(abs(x)) over [-1,+1].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N           Estimate             Error\n' );
  fprintf ( 1, '\n' );

  exact = 4.0 / 3.0;

  for n = 1 : 10

    [ x, w ] = clenshaw_curtis_set ( n );

    v(1:n,1) = sqrt ( abs ( x(1:n,1) ) );

    q = w' * v;
    e = abs ( q - exact );

    fprintf ( 1, '  %2d  %24.16g  %14.6e\n', n, q, e );

  end

  return
end


text:
function labels = grLabel(nodes, edges)
%GRLABEL Associate a label to each connected component of the graph
%
%   LABELS = grLabel(NODES, EDGES)
%   Returns an array with as many rows as the array NODES, containing index
%   number of each connected component of the graph. If the graph is
%   totally connected, returns an array of 1.
%
%   Example
%       nodes = rand(6, 2);
%       edges = [1 2;1 3;4 6];
%       labels = grLabel(nodes, edges);
%   labels =
%       1
%       1
%       1
%       2
%       3
%       2
%
%   See also
%   getNeighborNodes
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-08-14,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% init
Nn = size(nodes, 1);
labels = (1:Nn)';

% iteration until stability
modif = true;
while modif
    modif = false;

    % compute the minimum label in the neighborhood of each node
    for i = 1:Nn
        neigh = grAdjacentNodes(edges, i);
        neighLabels = labels([i;neigh]);

        % check for a modification
        if length(unique(neighLabels)) > 1
            modif = true;
        end

        % put new labels
        labels(ismember(labels, neighLabels)) = min(neighLabels);
    end
end

% renumbering to have fewer labels
labels2 = unique(labels);
for i = 1:length(labels2)
    labels(labels == labels2(i)) = i;
end


text:
filename='Bridge_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = VBA_spm_print(varargin)
% Print figure
% FORMAT spm_print(fname,F,opts)
% fname  - output filename [Default: 'spm_<date>']
% F      - figure handle or tag [Default: 'Graphics']
% opts   - structure containing printing options
%          [Default: defaults.ui.print from spm_defaults.m]
%
% FORMAT spm_print(job)
% Run a batch print job (see spm_cfg_print)
%__________________________________________________________________________
% Copyright (C) 1994-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_print.m 5969 2014-05-01 14:37:22Z guillaume $

%-Shortcut for list of graphics file formats available
%--------------------------------------------------------------------------
if nargin > 0 && isequal(lower(varargin{1}),'format')
    [varargout{1:nargout}] = print_format(varargin{2:end});
    return;
end

%-Conversion layer for spm.util.print batch
%--------------------------------------------------------------------------
if nargin == 1 && isstruct(varargin{1})
    job = varargin{1};
    VBA_spm_print(job.fname, job.fig.(char(fieldnames(job.fig))), job.opts);
    return;
end

%-Input parameters
%--------------------------------------------------------------------------
if nargin<1, fname = ''; else    fname = varargin{1}; end

if nargin<2, F = NaN;    else    F     = varargin{2}; end
if isnumeric(F) && ~isfinite(F), F     = [];          end
if isempty(F),                   F     = 'Graphics';  end

opts = spm_get_defaults('ui.print');
if ischar(opts)
    opts = VBA_spm_print('format',opts);
    if isempty(opts), disp('Print warning: unknown printing format.'); end
end
if nargin == 3
    o = varargin{3};
    if ischar(o), o = VBA_spm_print('format',o); end
    if isempty(o), disp('Print warning: unknown printing format.'); end
    try, opts.ext    = o.ext;    end
    try, opts.append = o.append; end
    try, opts.opt    = o.opt;    end
end
opts.opt = [opts.opt {'-noui','-painters'}];
if opts.append, opts.opt = [opts.opt {'-append'}]; end

%-Get output filename
%--------------------------------------------------------------------------
if isempty(fname)
    bname = ['spm_' datestr(now,'yyyymmmdd')];
    pname = pwd;
else
    bname = spm_file(fname,'basename');
    pname = spm_file(fname,'fpath');
end
fname = fullfile(pname,[bname opts.ext]);
if ~opts.append
    fname = spm_file(fname,'unique');
end

%-Get figure handle
%--------------------------------------------------------------------------
F = spm_figure('FindWin',F);
%if isempty(F), F = get(0,'CurrentFigure'); end
if isempty(F)
    disp('Print error: Figure not found.');
    return;
end

%-Print
%==========================================================================
try
    if ismember('-dfig',opts.opt)
        saveas(F, fname, 'fig');
    else
        if isdeployed
            deployprint(F, fname, opts.opt{:});
        else
            print(F, fname, opts.opt{:});
        end
    end
catch
    disp('Print error: nothing has been printed.');
    l = lasterror;
    disp(l.message);
    return;
end

%-Report
%--------------------------------------------------------------------------
str = '';
if ~isempty(get(F,'Tag')), str = sprintf('''%s'' ',get(F,'Tag')); end
fprintf('Printing %sfigure to:\n',str);
if ispc && ~ismember('-dfig',opts.opt), cmd = 'winopen(''%s'')';
else                                    cmd = 'open(''%s'')'; end
fprintf('  %s\n',spm_file(fname,'link',cmd));

%==========================================================================
% function pf = print_format(f)
%==========================================================================
function [pf, i] = print_format(f)

pf(1).name   = 'PostScript (PS)';
pf(1).opt    = {'-dpsc2'};
pf(1).append = true;
pf(1).ext    = '.ps';
pf(1).label  = {'ps'};

pf(2).name   = 'Encapsulated PostScript (EPS)';
pf(2).opt    = {'-depsc2'};
pf(2).append = false;
pf(2).ext    = '.eps';
pf(2).label  = {'eps'};

pf(3).name   = 'Portable Document Format (PDF)';
pf(3).opt    = {'-dpdf'};
pf(3).append = false;
pf(3).ext    = '.pdf';
pf(3).label  = {'pdf'};

pf(4).name   = 'JPEG image';
pf(4).opt    = {'-djpeg'};
pf(4).append = false;
pf(4).ext    = '.jpg';
pf(4).label  = {'jpg','jpeg'};

pf(5).name   = 'PNG image';
pf(5).opt    = {'-dpng'};
pf(5).append = false;
pf(5).ext    = '.png';
pf(5).label  = {'png'};

pf(6).name   = 'TIFF image';
pf(6).opt    = {'-dtiff'};
pf(6).append = false;
pf(6).ext    = '.tif';
pf(6).label  = {'tif','tiff'};

pf(7).name   = 'MATLAB figure';
pf(7).opt    = {'-dfig'}; % use saveas instead of print
pf(7).append = false;
pf(7).ext    = '.fig';
pf(7).label  = {'fig'};

if nargin
    for i=1:numel(pf)
        if ismember(lower(f),lower(pf(i).label))
            pf = pf(i); return;
        end
    end
    pf = []; i = [];
end

author: MBB-team

repo: VBA-toolbox


text:
function train_model = runTrainNRSFM(data,class)

% Run and cache NRSFM. Load cached NRSFM model if available
globals;
params = get_params();

% Path to cached NRSFM model
shapeModelNRSFMFile = fullfile(cachedir,class,'shapeModelNRSFM.mat');

fprintf('\n%%%%%%%%%%%% NRSFM %%%%%%%%%%%%\n');

% Run NRSFM if cached model not found
if(exist(shapeModelNRSFMFile,'file'))
    fprintf('Loading cached NRSFM model from \n%s\n',shapeModelNRSFMFile);
    load(shapeModelNRSFMFile,'train_model');
else
    % Setup training data struct (filter instances and normalize keypoints)
    [train,~] = prep_data(data.train,1,params.nrsfm.flip);

    % Train model
    fprintf('Train Non Rigid SFM model \n');
    train_model = train_nrsfm(params.class,train);

    % Caching trained NRSFM model
    fprintf('Caching NRSFM model in \n%s\n',shapeModelNRSFMFile);
    save(shapeModelNRSFMFile,'train_model');
end

end

repo: CategoryShapes


text:
function hf=sliceViewer(varargin)

%% Parse input

switch nargin
    case 1
        M=varargin{1};
        v=ones(1,3);
        viewerType=1;
    case 2
        M=varargin{1};
        v=varargin{2};
        viewerType=1;
    case 3
        M=varargin{1};
        v=varargin{2};
        viewerType=varargin{3};
end

%% Start viewer

switch viewerType
    case 1
        hf=sv2(M,v);
    case 2
        hf=sv3(M,v);
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [V,F] = readOBJfast(filename,varargin)
  % readOBJfast
  %
  % reads an OBJ file quickly, but OBJ file should be formated *simply*
  %
  % Input:
  %  filename  path to .obj file
  % Outputs:
  %  V  #V by 3 list of vertices
  %  F  #F by 3 list of triangle indices
  %
  % filename should be a text file containing first vertex positions where each
  % line is:
  % v x y z
  % then a list of faces where each line is:
  % f i j k
  % and exactly that. Any lines before vertices, or after vertices but before
  % faces are ignored. Comments throughout vetex and face lines will surely
  % break this.
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % For complete(r) support use readOBJ
  %
  % See also readOBJ
  %

  function file_wrangler()
    error_handler = @final_error;
    if isunix()
      % try to make the file easier to read
      tmpf = tempname();
      escaped_filename = filename;
      escaped_filename = strrep(escaped_filename,' ','\ ');
      escaped_filename = strrep(escaped_filename,'(','\(');
      escaped_filename = strrep(escaped_filename,')','\)');
      cmd = sprintf( ...
        ['grep "^v " %s >%s ; ' ...
          'grep "^f " %s | ' ...
          'sed -e "s/\\([0-9]*\\)\\(\\/[0-9]*\\)*/\\1/g" >>%s'], ...
        escaped_filename,tmpf,escaped_filename,tmpf);
      s = system(cmd);
      if s==0
        readOBJfast_helper(tmpf);
        if any(F<=0)
          final_error();
        end
        %delete(tmpf);
        return;
      end
    end
    final_error();
  end

  function final_error()
    error('Bad format. Try readOBJ...\n');
  end

  function readOBJfast_helper(filename)
    fp = fopen(filename);
    % read in vertices
    V = [];
    while(true)
      V = fscanf(fp,' v %lg %lg %lg ',inf);
      if(prod(size(V)) > 0)
        break;
      else
        line = fgets(fp);
        if(prod(size(line)) == 0)
          fclose(fp);
          error_handler();
          return;
        end
      end
    end
    V = reshape(V,3,size(V,1)/3)';
    % read in faces
    F = [];
    if quads
      format = ' f %d %d %d %d ';
    else
      format = ' f %d %d %d ';
    end
    while(true)
      F = fscanf(fp,format,inf);
      if(numel(F) > 0)
        break;
      else
        line = fgets(fp);
        if(numel(line) == 0 || (isnumeric(line) && line == -1))
          fclose(fp);
          error_handler();
          return;
        end
      end
    end
    if quads
      if mod(size(F,1),4) ~= 0
        fclose(fp);
        error_handler();
        return;
      end
      F = reshape(F,4,size(F,1)/4)';
    else
      if mod(size(F,1),3) ~= 0
        fclose(fp);
        error_handler();
        return;
      end
      F = reshape(F,3,size(F,1)/3)';
    end
    if ~feof(fp)
      fclose(fp);
      error_handler();
      return;
    end
    fclose(fp);
  end

  % default values
  quads = false;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Quads'}, ...
    {'quads'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  error_handler = @file_wrangler;
  readOBJfast_helper(filename)
  if all(F(:)<0)
    F = F+size(V,1)+1;
  elseif any(F(:)<0)
    final_error();
  end
end

author: alecjacobson

repo: gptoolbox


text:
function areaV = calculate_structure_cross_sectional_area(structIndV,planC)
% function areaV = calculate_structure_cross_sectional_area(structIndV,planC)
%
% APA, 08/29/2012

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

if ischar(structIndV)
    structIndV = {structIndV};
end

structureNamesC = {planC{indexS.structures}.structureName};
for i = 1:length(structIndV)
    if isnumeric(structIndV(i))
        structNum = structIndV(i);
    else
        structureName = structIndV{i};
        structNum = getMatchingIndex(structureName,structureNamesC,'regex');
        if length(structNum) > 1
            structNum = structNum(1);
        end
    end

    if isempty(structNum) || length(structNum) > 1
        areaV(i) = NaN;
        continue
    end

    [rasterSegments, planC, isError] = getRasterSegments(structNum,planC);
    if isempty(rasterSegments)
        areaV(i) = NaN;
        continue;
    end
    scanNum = getStructureAssociatedScan(structNum,planC);
    [~, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
    middleSlcNum = round(median(uniqueSlices));

    areaVal = 0;
    for segNum = 1:length(planC{indexS.structures}(structNum).contour(middleSlcNum).segments)
        areaVal = areaVal + polyarea(planC{indexS.structures}(structNum).contour(middleSlcNum).segments(segNum).points(:,1),planC{indexS.structures}(structNum).contour(middleSlcNum).segments(segNum).points(:,2));
    end

    areaV(i) = areaVal;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function ind = randompoint(prob, n)
%RANDOMNEW to generate n new point randomly from the mop problem given.

if (nargin==1)
    n=1;
end

randarray = rand(prob.pd, n);
lowend = prob.domain(:,1);
span = prob.domain(:,2)-lowend;
point = randarray.*(span(:,ones(1, n)))+ lowend(:,ones(1,n));
cellpoints = num2cell(point, 1);

indiv = struct('parameter',[],'objective',[], 'estimation', []);
ind = repmat(indiv, 1, n);
[ind.parameter] = cellpoints{:};

% estimation = struct('obj', NaN ,'std', NaN);
% [ind.estimation] = deal(repmat(estimation, prob.od, 1));
end

author: sfvsfv

repo: Mathematical-modeling


text:
function X = ttm(X,V,varargin)
%TTM Tensor times matrix for ktensor.
%
%   Y = TTM(X,A,N) computes the n-mode product of the ktensor X with a
%   matrix A; i.e., X x_N A.  The integer N specifies the dimension
%   (or mode) of X along which A should be multiplied.  If size(A) =
%   [J,I], then X must have size(X,N) = I.  The result will be a
%   ktensor of the same order and size as X except that size(Y,N) = J.
%
%   Y = TTM(X,{A,B,C,...}) computes the n-mode product of the ktensor
%   X with a sequence of matrices in the cell array.  The n-mode
%   products are computed sequentially along all dimensions (or modes)
%   of X. The cell array contains ndims(X) matrices.
%
%   Y = TTM(X,{A,B,C,...},DIMS) computes the sequence tensor-matrix
%   products along the dimensions specified by DIMS.
%
%   Y = TTM(...,'t') performs the same computations as above except
%   the matrices are transposed.
%
%   Examples
%   X = ktensor({rand(5,2),rand(3,2),rand(4,2),rand(2,2)});
%   A = rand(4,5); B = rand(4,3); C = rand(3,4); D = rand(3,2);
%   Y = ttm(X, A, 1)         %<-- computes X times A in mode-1
%   Y = ttm(X, {A,B,C,D}, 1) %<-- same as above
%   Y = ttm(X, A', 1, 't')   %<-- same as above
%   Y = ttm(X, {A,B,C,D}, [1 2 3 4]) %<-- 4-way multiply
%   Y = ttm(X, {D,C,B,A}, [4 3 2 1]) %<-- same as above
%   Y = ttm(X, {A,B,C,D})            %<-- same as above
%   Y = ttm(X, {A',B',C',D'}, 't')   %<-- same as above
%   Y = ttm(X, {C,D}, [3 4])     %<-- X times C in mode-3 & D in mode-4
%   Y = ttm(X, {A,B,C,D}, [3 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, [1 2 4])   %<-- 3-way multiply
%   Y = ttm(X, {A,B,C,D}, [1 2 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, -3)        %<-- same as above
%   Y = ttm(X, {A,B,C,D}, -3)      %<-- same as above
%
%   See also KTENSOR, KTENSOR/TTV, TENSOR/TTM.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

%%%%%%%%%%%%%%%%%%%%%%
%%% ERROR CHECKING %%%
%%%%%%%%%%%%%%%%%%%%%%

% Check the number of arguments
if (nargin < 2)
    error('TTM requires at least two arguments.');
end

% Check for transpose option
isTranspose = false;
if numel(varargin) > 0
  if isnumeric(varargin{1});
    dims = varargin{1};
  end
  isTranspose =  (ischar(varargin{end}) && (varargin{end} == 't'));
end

% Check for dims argument
if ~exist('dims','var')
    dims = [];
end

% Check that 2nd argument is cell array. If not, recall with V as a
% cell array with one element.
if ~iscell(V)
    X = ttm(X,{V},dims,varargin{end});
    return;
end

% Get sorted dims and index for multiplicands
[dims,vidx] = tt_dimscheck(dims,ndims(X),numel(V));

% Determine correct size index
if isTranspose
  j = 1;
else
  j = 2;
end

% Check that each multiplicand is the right size.
for i = 1:numel(dims)
    if (ndims(V) ~= 2) || (size(V{vidx(i)},j) ~= size(X,dims(i)))
disp(size(V{vidx(i)}))
disp(size(X))

        error('Multiplicand is wrong size');
    end
end

% Do the multiplications in the specified modes.
for i = 1:numel(dims)
  if isTranspose
    X.u{dims(i)} = V{vidx(i)}'* X.u{dims(i)};
  else
    X.u{dims(i)} = V{vidx(i)} * X.u{dims(i)};
  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%  INTERNAL FUNCTION: Creates a handle to an m-file not on the matlab search path
%
%  ::
%
%    hdl = func2fhandle(full_path/func)
%    hdl = func2fhandle(relative_path/func)
%    [hdl1,hdl2,...,hdln] = func2fhandle(p1,p2,...,pn)
%    [hdl1,hdl2,...,hdln] = func2fhandle({p1,p2,...,pn})
%
%  Args:
%
%     - **pth** [char|cellstr] : path names to create handles for. The path
%       could be relative to the current directory or absolute
%
%  Returns:
%     :
%
%     - **hdl** [function_handle] : function handle(s)
%
%

author: jmaih


function bool = hasTestFailed(status)
    % returns true when the test has failed

    if iscell(status) % allow for vectorization of the call
        bool = cellfun(@hasTestFailed, status, 'UniformOutput', true);
    else
        stages = getStagesFromStatus(status);
        bool = false;
        for jStage = 1:numel(stages)
            bool = bool || status.(stages{jStage}).error;
        end
    end
end