text:
function [ n_data_new, n, m, x, fx ] = legendre_associated_values ( n_data )

%*****************************************************************************80
%
%% LEGENDRE_ASSOCIATED_VALUES returns values of associated Legendre functions.
%
%  Discussion:
%
%    The function considered is the associated Legendre function P^M_N(X).
%
%    In Mathematica, the function
%
%      LegendreP [ n, m, x ]
%
%    returns the value of the associated Legendre function P^M_N(X).
%
%  Differential equation:
%
%    (1-X*X) * Y'' - 2 * X * Y + ( N (N+1) - (M*M/(1-X*X)) * Y = 0
%
%  First terms:
%
%    M = 0  ( = Legendre polynomials of first kind P(N)(X) )
%
%    P00 =    1
%    P10 =    1 X
%    P20 = (  3 X^2 -   1)/2
%    P30 = (  5 X^3 -   3 X)/2
%    P40 = ( 35 X^4 -  30 X^2 +   3)/8
%    P50 = ( 63 X^5 -  70 X^3 +  15 X)/8
%    P60 = (231 X^6 - 315 X^4 + 105 X^2 -  5)/16
%    P70 = (429 X^7 - 693 X^5 + 315 X^3 - 35 X)/16
%
%    M = 1
%
%    P01 =   0
%    P11 =   1 * SQRT(1-X*X)
%    P21 =   3 * SQRT(1-X*X) * X
%    P31 = 1.5 * SQRT(1-X*X) * (5*X*X-1)
%    P41 = 2.5 * SQRT(1-X*X) * (7*X*X*X-3*X)
%
%    M = 2
%
%    P02 =   0
%    P12 =   0
%    P22 =   3 * (1-X*X)
%    P32 =  15 * (1-X*X) * X
%    P42 = 7.5 * (1-X*X) * (7*X*X-1)
%
%    M = 3
%
%    P03 =   0
%    P13 =   0
%    P23 =   0
%    P33 =  15 * (1-X*X)^1.5
%    P43 = 105 * (1-X*X)^1.5 * X
%
%    M = 4
%
%    P04 =   0
%    P14 =   0
%    P24 =   0
%    P34 =   0
%    P44 = 105 * (1-X*X)^2
%
%  Recursion:
%
%    if N < M:
%      P(N,M) = 0
%    if N = M:
%      P(N,M) = (2*M-1)!! * (1-X*X)^(M/2) where N!! means the product of
%      all the odd integers less than or equal to N.
%    if N = M+1:
%      P(N,M) = X*(2*M+1)*P(M,M)
%    if M+1 < N:
%      P(N,M) = ( X*(2*N-1)*P(N-1,M) - (N+M-1)*P(N-2,M) )/(N-M)
%
%  Restrictions:
%
%    -1 <= X <= 1
%     0 <= M <= N
%
%  Special values:
%
%    P(N,0)(X) = P(N)(X), that is, for M=0, the associated Legendre
%    function of the first kind equals the Legendre polynomial of the
%    first kind.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, integer M, real X, the arguments of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 19;
  n_vec = [ ...
     1,  1,  1,  1, ...
     1,  2,  2,  2, ...
     3,  3,  3,  3, ...
     4,  5,  6,  7, ...
     8,  9, 10 ];
  m_vec = [ ...
     0,  0,  0,  0, ...
     1,  0,  1,  2, ...
     0,  1,  2,  3, ...
     2,  2,  3,  3, ...
     4,  4,  5 ];
  fx_vec = [ ...
     0.000000E+00,  0.500000E+00,  0.707107E+00,  1.000000E+00, ...
    -0.866025E+00, -0.125000E+00, -1.29904E+00,   2.25000E+00, ...
    -0.437500E+00, -0.324759E+00,  5.62500E+00,  -9.74278E+00, ...
     4.21875E+00,  -4.92187E+00,   12.7874E+00,   116.685E+00, ...
    -1050.67E+00,  -2078.49E+00,   30086.2E+00 ];
  x_vec = [ ...
    0.0E+00,       0.5E+00,       0.7071067E+00, 1.0E+00, ...
    0.5E+00,       0.5E+00,       0.5E+00,       0.5E+00, ...
    0.5E+00,       0.5E+00,       0.5E+00,       0.5E+00, ...
    0.5E+00,       0.5E+00,       0.5E+00,       0.5E+00, ...
    0.5E+00,       0.5E+00,       0.5E+00 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    m = 0;
    x = 0.0E+00;
    fx = 0.0E+00;
  else
    n = n_vec(n_data_new);
    m = m_vec(n_data_new);
    x = x_vec(n_data_new);
    fx = fx_vec(n_data_new);
  end

  return
end


text:
function li = lagrange_value ( nd, xd, ni, xi )

%*****************************************************************************80
%
%% LAGRANGE_VALUE evaluates the Lagrange basis polynomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(ND,1), the data nodes.
%
%    Input, integer NI, the number of evaluation points.
%
%    Input, real XI(NI,1), the evaluation points.
%
%    Output, real LI(NI,ND), the value, at the I-th point XI, of the
%    Jth basis function.
%
  li = zeros ( ni, nd );

  for i = 1 : ni
    for j = 1 : nd
      li(i,j) = prod ( ( xi(i) - xd(1:j-1)  ) ./ ( xd(j) - xd(1:j-1)  ) ) ...
              * prod ( ( xi(i) - xd(j+1:nd) ) ./ ( xd(j) - xd(j+1:nd) ) );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = niederreiter_mccurley4 ( x )

%*****************************************************************************80
%
%% NIEDERREITER_MCCURLEY4 evaluates the Niederreiter-McCurley function 4.
%
%  Discussion:
%
%    The spatial dimension is taken as 4.
%
%    The maximum is sought within the unit hypercube.
%
%    The reported maximum value occurs at approximately
%
%      X = [ 1, 1, 1, 1 ]
%
%    with
%
%      F(X) = 9.196986029286059
%
%    Note that FMINCON is designed to minimize functions, so we must negate
%    the output value here in order to satisfy FMINCON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, real X(N), the argument of the objection function.
%
%    Output, real F, the value of the objective function.
%
  f = - 100 * prod ( x ) * exp ( - x(4) ) ./ ( 1 + x(1) * x(2) * x(3) ).^2;

  return
end


text:
function [] = threed_to_vu(fname,x,e_conn,u,var)

%*****************************************************************************80
%
%  threed_to_vu.m - Writes out a 3D finite element solution to a format
%                readable by Vu
%                (http://www.invisu.ca/e/index.php)
%
%  Copyright (c) 2002, Jeff Borggaard, Virginia Tech
%  Version: 1.0
%
%  Usage:    [] = threed_tovu(fname,x,e_conn,u,var)
%
%  Variables:     fname
%                        Name of the output file (default="femvu.vu")
%                 x
%                        Nodal coordinates
%                 e_conn
%                        Element connectivity
%                 u
%                        Columns of scalar variables to write to Vu file
%                 var(n_var).names
%                        List of scalar variable names (default="var#")
%-----------------------------------------------------------------------
  [n_nodes,    n_dim      ] = size(x);
  [n_elements, nel_dof    ] = size(e_conn);
  [temp      , n_variables] = size(u);
  var(n_variables+1).name = '';

  if ( isempty(fname) )
    fname = 'femvu.vu';
  end
  fid = fopen(fname,'w');

  % Write out coordinate values
  fprintf(fid,'FIELD x_coord( ) =\n {\n');
  for i=1:n_nodes
    fprintf(fid,' %10.6e %10.6e %10.6e\n',x(i,1),x(i,2),x(i,3));
  end
  fprintf(fid,'};\n\n');

  % Write out element connectivity
  fprintf(fid,'FIELD<int> e_conn( ) =\n{\n');
  if ( nel_dof==4 ) % Test for linear tetrahedra
    for i=1:n_elements
      fprintf(fid,' %6i %6i %6i %6i\n',e_conn(i,:));
    end
  elseif ( nel_dof==10 ) % Test for quadratic tetrahedra
    for i=1:n_elements
      fprintf(fid,' %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i\n',...
              e_conn(i,[1 2 3 4 5 6 7 10 8 9]));
    end
  else
    error('threed_to_vu: element type not currently implemented\n')
    return
  end
  fprintf(fid,'};\n\n');

  % Write out variables
  for n_var=1:n_variables
    if ( isempty(var(n_var).name) )
      var(n_var).name = strcat('Var',int2str(n_var));
    end

    temp = strcat('FIELD@',var(n_var).name,'( ) = {\n');
    temp = strrep(temp,'D@','D ');
    fprintf(fid,temp);
    fprintf(fid,' %10.5e %10.5e %10.5e %10.5e %10.5e %10.5e \n',...
                  u(:,n_var) );
    fprintf(fid,'\n};\n');
  end
  fprintf(fid,'\n\n');

  % Write out zones and connectivity
  if ( nel_dof==4 )
    fprintf(fid,'MESH MyMesh( ) =\n{\n');
    fprintf(fid,'   ZONE Zone1( LagrTetra04, x_coord, e_conn );\n' );
    fprintf(fid,'};\n\n');
  elseif ( nel_dof==10 )
    fprintf(fid,'MESH MyMesh( ) =\n{\n');
    fprintf(fid,'   ZONE Zone1( LagrTetra10, x_coord, e_conn );\n' );
    fprintf(fid,'};\n\n');
  end

  fprintf(fid,'SOLUTION MySolution( ) =\n{\n');
  if ( nel_dof==4 )
    for n_var=1:n_variables
      temp = strcat('   VARIABLE@'   , var(n_var).name,...
                    '( LagrTetra04,@', var(n_var).name,...
                    ', e_conn, Zone1 );\n');
      temp = strrep(temp,'@',' ');
      fprintf(fid,temp);
    end
  else
    for n_var=1:n_variables
      temp = strcat('   VARIABLE@'   , var(n_var).name,...
                    '( LagrTetra10,@', var(n_var).name,...
                    ', e_conn, Zone1 );\n');
      temp = strrep(temp,'@',' ');
      fprintf(fid,temp);
    end
  end
  fprintf(fid,'};\n');

  fclose(fid);


text:
function d = NormL1NN_dual(x,weights)
% Dual of non-negative L1 gauge function

x(x < 0) = 0;
d = norm(x./weights,inf);

author: mpf


text:
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
J=fftshift(fft2(I));
[x, y]=meshgrid(-128:127, -128:127);
z=sqrt(x.^2+y.^2);
D1=10;  D2=30;
n=6;
H1=1./(1+(z/D1).^(2*n));
H2=1./(1+(z/D2).^(2*n));
K1=J.*H1;
K2=J.*H2;
L1=ifft2(ifftshift(K1));
L2=ifft2(ifftshift(K2));
figure;
subplot(131);
imshow(I);
subplot(132);
imshow(real(L1));
subplot(133);
imshow(real(L2))

author: HuangCongQing


text:
function [FV] = mesh_bet2matlab(fileprefix)

% mesh_bet2matlab - Read FSL BET tesselation (.coo/.dat)
%
% USEAGE: [FV] = mesh_bet2matlab(fileprefix)
%
% This function will load ascii files created by the FSL
% BET function, see http://www.fmrib.ox.ac.uk/fsl
%
% When used with the -x option, the FSL BET function can
% output the tesselation of the brain surface (ie, only
% a fairly smooth, non-detailed, surface), eg:
%
% bet T1 T1_brain -x
%
% The brain tesselation is output into two files, eg:
%
%   T1_brain.coo
%   T1_brain.dat
%
% The T1_brain.coo file contains the co-ordinates of the
% vertices, and the .dat file contains the details of
% the links between vertices (each line corresponds to
% each vertex in the .coo file, and the 1's correspond
% to which other vertices they are connected to).
%
% The returned FV struct contains the 'vertices' and
% 'faces' matrices, which can be input to the patch
% command, eg:
%
% Hpatch = patch('Vertices',FV.vertices,'Faces',FV.faces,...
%                'EdgeColor',[.8 .8 .8],...
%                'FaceColor',[0.9 0.9 0.9]);
%
% This will plot the mesh as a patch object.  See the patch
% command and matlab help for more information on coloring
% this object.
%
% See also: mesh_freesurfer2matlab
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  03/02 Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('fileprefix','var'),
    error('No input fileprefix');
end

if findstr('.coo',fileprefix),
    fprintf('MESH_BET2MATLAB: Removing .coo extension from ''%s''\n',fileprefix);
    fileprefix = strrep(fileprefix,'.coo','');
end
if findstr('.dat',fileprefix),
    fprintf('MESH_BET2MATLAB: Removing .dat extension from ''%s''\n',fileprefix);
    fileprefix = strrep(fileprefix,'.dat','');
end

FV.vertices = readCOO(fileprefix);

Nvertices   = size(FV.vertices,1);

edgematrix  = readDAT(fileprefix,Nvertices);

FV.faces    = edge2face(FV.vertices,edgematrix);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function faces = edge2face(vertices,edgematrix)

    Nvertices = size(vertices,1);

    plot3(vertices(:,1),vertices(:,2),vertices(:,3),'b.'), hold on

    for i = 1:Nvertices,

        % Get current vertex XYZ coordinates
        v  = vertices(i,:);

        plot3(v(:,1),v(:,2),v(:,3),'go');

        % Get neighbouring vertex XYZ coordinates
        vn = vertices(find(edgematrix(i,:)),:);

        plot3(vn(:,1),vn(:,2),vn(:,3),'ro');

        faces = tri([v;vn]);

        return

    end

    faces = [];

return

function faces = tri(V)

    % Must identify angle of rotation from
    % central vertex, V(1,:), to all other vertices
    % in a counterclockwise direction

    % First translate all vertices to an origin
    % given by the central vertex
    Vo = V - repmat(V(1,:),size(V,1),1);

    % Find direction cosines for line from centre to vertex
    d = zeros(size(Vo,1),1);
    l = zeros(size(Vo,1),1);
    m = zeros(size(Vo,1),1);
    n = zeros(size(Vo,1),1);

    for i = 1:length(Vo),

        x = Vo(i,1);  y = Vo(i,2);  z = Vo(i,3);

        d(i) = sqrt( (x)^2 + (y)^2 + (z)^2 );

        if d(i) > 0,
            l(i) = x/d(i); % cos alpha
            m(i) = y/d(i); % cos beta
            n(i) = z/d(i); % cos gamma
        end
    end

    faces = [];

    %hold off
    %plot3(Vo(:,1),Vo(:,2),Vo(:,3),'b.'), view(2), hold on

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vertices = readCOO(fileprefix)

	[path,name,ext] = fileparts(strcat(fileprefix,'.coo'));
	file = fullfile(path,[name ext]);

	fid = fopen(file,'r');

	if isequal(fid,-1),
        S=sprintf('Could not open file: "%s"',file);
        error(S);
	else
        fprintf('...Reading BET vertices...');
        tic;

        % Read vertices
        vertices = fscanf(fid,'%s%f%f%f',[4,inf]);
        fclose(fid);

        % remove last row (all zeros) and translate
        vertices = vertices(2:4,:)';

        t = toc;
        fprintf('...done (%6.2f sec).\n',t);
	end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function edgematrix = readDAT(fileprefix,Nvertices)

	[path,name,ext] = fileparts(strcat(fileprefix,'.dat'));
	file = fullfile(path,[name ext]);

	fid = fopen(file,'r');

	if isequal(fid,-1),
        S=sprintf('Could not open file: "%s"',file);
        error(S);
	else
        fprintf('...Reading BET edge matrix...');
        tic;

        % Read faces
        edgematrix = fscanf(fid,'%1d',[Nvertices,Nvertices]);
        fclose(fid);

        t = toc;
        fprintf('...done (%6.2f sec).\n',t);
	end
return

author: PatternRecognition

repo: OpenBMI


text:
function [ x, seed ] = burr_sample ( a, b, c, d, seed )

%*****************************************************************************80
%
%% BURR_SAMPLE samples the Burr PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, the parameters of the PDF.
%    0 < B,
%    0 < C.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = burr_cdf_inv ( cdf, a, b, c, d );

  return
end


text:
function ret = wrap_boundary_liu(img, img_size)
% wrap_boundary_liu.m
%
%   pad image boundaries such that image boundaries are circularly smooth
%
%     written by Sunghyun Cho (sodomau@postech.ac.kr)
%
% This is a variant of the method below:
%   Reducing boundary artifacts in image deconvolution
%     Renting Liu, Jiaya Jia
%     ICIP 2008
%
    [H, W, Ch] = size(img);
    H_w = img_size(1) - H;
    W_w = img_size(2) - W;

    ret = zeros(img_size(1), img_size(2), Ch);
    for ch = 1:Ch
        alpha = 1;
        HG = img(:,:,ch);

        r_A = zeros(alpha*2+H_w, W);
        r_A(1:alpha, :) = HG(end-alpha+1:end, :);
        r_A(end-alpha+1:end, :) = HG(1:alpha, :);
        a = ((1:H_w)-1)/(H_w-1);
        r_A(alpha+1:end-alpha, 1) = (1-a)*r_A(alpha,1) + a*r_A(end-alpha+1,1);
        r_A(alpha+1:end-alpha, end) = (1-a)*r_A(alpha,end) + a*r_A(end-alpha+1,end);

        A2 = solve_min_laplacian(r_A(alpha:end-alpha+1,:));
        r_A(alpha:end-alpha+1,:) = A2;
        A = r_A;

        r_B = zeros(H, alpha*2+W_w);
        r_B(:, 1:alpha) = HG(:, end-alpha+1:end);
        r_B(:, end-alpha+1:end) = HG(:, 1:alpha);
        a = ((1:W_w)-1)/(W_w-1);
        r_B(1, alpha+1:end-alpha) = (1-a)*r_B(1,alpha) + a*r_B(1,end-alpha+1);
        r_B(end, alpha+1:end-alpha) = (1-a)*r_B(end,alpha) + a*r_B(end,end-alpha+1);

        B2 = solve_min_laplacian(r_B(:, alpha:end-alpha+1));
        r_B(:,alpha:end-alpha+1,:) = B2;
        B = r_B;

        r_C = zeros(alpha*2+H_w, alpha*2+W_w);
        r_C(1:alpha, :) = B(end-alpha+1:end, :);
        r_C(end-alpha+1:end, :) = B(1:alpha, :);
        r_C(:, 1:alpha) = A(:, end-alpha+1:end);
        r_C(:, end-alpha+1:end) = A(:, 1:alpha);

        C2 = solve_min_laplacian(r_C(alpha:end-alpha+1, alpha:end-alpha+1));
        r_C(alpha:end-alpha+1, alpha:end-alpha+1) = C2;
        C = r_C;

        A = A(alpha:end-alpha-1, :);
        B = B(:, alpha+1:end-alpha);
        C = C(alpha+1:end-alpha, alpha+1:end-alpha);

        ret(:,:,ch) = [img(:,:,ch) B; A C];
    end

end

function [img_direct] = solve_min_laplacian(boundary_image)
% function [img_direct] = poisson_solver_function(gx,gy,boundary_image)
% Inputs; Gx and Gy -> Gradients
% Boundary Image -> Boundary image intensities
% Gx Gy and boundary image should be of same size
    [H,W] = size(boundary_image);

    % Laplacian
    f = zeros(H,W);                          clear j k

    % boundary image contains image intensities at boundaries
    boundary_image(2:end-1, 2:end-1) = 0;
    j = 2:H-1;      k = 2:W-1;      f_bp = zeros(H,W);
    f_bp(j,k) = -4*boundary_image(j,k) + boundary_image(j,k+1) + ...
        boundary_image(j,k-1) + boundary_image(j-1,k) + boundary_image(j+1,k);
    clear j k

    %f1 = f - reshape(f_bp,H,W); % subtract boundary points contribution
    f1 = f - f_bp; % subtract boundary points contribution
    clear f_bp f

    % DST Sine Transform algo starts here
    f2 = f1(2:end-1,2:end-1);                   clear f1
    % compute sine tranform
    tt = dst(f2);       f2sin = dst(tt')';      clear f2

    % compute Eigen Values
    [x,y] = meshgrid(1:W-2, 1:H-2);
    denom = (2*cos(pi*x/(W-1))-2) + (2*cos(pi*y/(H-1)) - 2);

    % divide
    f3 = f2sin./denom;                          clear f2sin x y

    % compute Inverse Sine Transform
    tt = idst(f3);      clear f3;       img_tt = idst(tt')';        clear tt

    % put solution in inner points; outer points obtained from boundary image
    img_direct = boundary_image;
    img_direct(2:end-1,2:end-1) = 0;
    img_direct(2:end-1,2:end-1) = img_tt;
end

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function basis = basis_col(A)
%
% Bases
%
% basis_col(A) produces a basis for the subspace of Eucldiean n-space
% spanned by the vectors {u1,u2,...}, where the matrix A is formed from
% these vectors as its columns. That is, the subspace is the column space
% of A. The columns of the matrix that is returned are the basis vectors
% for the subspace. These basis vectors will be a subset of the original
% vectors. An error is returned if a basis for the zero vector space is
% attempted to be produced.
%
% For example, if the vector space V = span{u1,u2,...}, where u1,u2,... are
% row vectors, then set A to be [u1' u2' ...].
%
% For example, if the vector space V = Col(B), where B is an m x n matrix,
% then set A to be equal to B.

matrix_size = size(A);

m = matrix_size(1,1);
n = matrix_size(1,2);

if A == zeros(m,n)
    error('There does not exist a basis for the zero vector space.');
elseif n == 1
    basis = A;
else
    flag = 0;

    if n == 2
        multiple = A(1,2)/A(1,1);
        count = 0;

        for i = 1:m
            if A(i,2)/A(i,1) == multiple
                count = count + 1;
            end
        end

        if count == m
            basis = A(1:m,1);
            flag = 1;
        end
    end

    if flag == 0
        [ref_A pivot_columns] = ref(A);

        for i = 1:size(pivot_columns,2)
            B(1:m,i) = A(1:m,pivot_columns(1,i));
        end

        basis = B;
    end
end

author: xiuyechen

repo: FishExplorer


text:
function varargout = volumization(varargin)
% GUI for visualizing 3D data
%
% To run type:
% >>volumization (data,'n')
% for example:voumization (data, '2')
% where data are of format [n x m x p]  and  'n' - which slice to display:
% '1'- all, '2' - every second, '3' - every third, etc.
%
% For users with PLS_Toolbox 3.5 can type
% >> volumization
% and then load data interactively using lddlgpls routine.
%
% Includes different types of visualization:
% 1. Orthogonal slices - each slice individually, as well as three orthogonal slices on one plot
% 2. Rendered volume through showing series of z-slices. Enter the number of slices to display:
% ('1' - all, '2' - every second, etc.)
% 3. Isosurfaces:
%	- one for a single value
%	- multiple for multiple values of interest on one figure (each with different color)
%
% Options to change for orthogonal slices (1) and volumes (2):
% - transparency;
% - color scheme.
%
% Options to changes for all three types of visualization (1-3):
% - top/bottom, left/right
% - aspect ratio.
%
% To display isosurface on top of rendered volume or orthogonal slices, use controls to
% display desired view and then use Multiple isosurface slider to display the isosurface
% of value of interest. (Single isosurface slider clears the figure on execution).
%
% Use Clear figure button to start over.
%
% created by K.Artyushkova
% November 2004
%
% Kateryna Artyushkova
% Postdoctoral Scientist
% Department of Chemical and Nuclear Engineering
% The University of New Mexico
% (505) 277-0750
% kartyush@unm.edu
%
% Last Modified by GUIDE v2.5 02-Feb-2006 11:38:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @volumization_OpeningFcn, ...
                   'gui_OutputFcn',  @volumization_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before volumization is made visible.
function volumization_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to volumization (see VARARGIN)

% Choose default command line output for volumization
handles.output = hObject;

% Update handles structure
if (nargin <4)
        H.Position=[181 477 332 275];
figure(H)
msgbox('Please load the data through the File Load menu')
    else
    data= varargin{1};
    display=varargin{2};
    data=double(data);
[m,n,p]=size(data);
set(handles.xmin,'string',1);
set(handles.xmax,'string',n);
set(handles.xi,'string',1);
set(handles.ymin,'string',1);
set(handles.ymax,'string',m);
set(handles.yi,'string',1);
set(handles.zmin,'string',1);
set(handles.zmax,'string',p);
set(handles.zi,'string',1);
s=data;
handles.s=s;
H.Position=[181 477 332 275];
figure(H)
step=str2double(display);
hz=slice(s,[],[],[1:step:p]);
alpha('color')
axis tight
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
alphamap('rampdown')
daspect([1 1 0.5])
handles.aspect=[1 1 0.5];
handles.count=0;
colormap(jet)

end
handles.v={'0'};
guidata(hObject, handles);

% UIWAIT makes volumization wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = volumization_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --------------------------------------------------------------------
function load_data_Callback(hObject, eventdata, handles)
% hObject    handle to load_data (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data =lddlgpls;
data=double(data);
[m,n,p]=size(data);
set(handles.xmin,'string',1);
set(handles.xmax,'string',n);
set(handles.xi,'string',1);
set(handles.ymin,'string',1);
set(handles.ymax,'string',m);
set(handles.yi,'string',1);
set(handles.zmin,'string',1);
set(handles.zmax,'string',p);
set(handles.zi,'string',1);
s=data;
handles.s=data;
figure(1)
hz=slice(s,[],[],[1:3:p]);
alpha('color')
axis tight
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
alphamap('rampdown')
daspect([1 1 0.5])
handles.aspect=[1 1 0.5];
handles.count=0;
colormap(jet)
guidata(hObject, handles);

% --- Executes on button press in alphaincrease.
function alphaincrease_Callback(hObject, eventdata, handles)
% hObject    handle to alphaincrease (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure(1)
alphamap('increase',.1)

% --- Executes on button press in alphadecrease.
function alphadecrease_Callback(hObject, eventdata, handles)
% hObject    handle to alphadecrease (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure(1)
alphamap('decrease',.1)

% --- Executes during object creation, after setting all properties.
function colormap_CreateFcn(hObject, eventdata, handles)
% hObject    handle to colormap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in colormap.
function colormap_Callback(hObject, eventdata, handles)
% hObject    handle to colormap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns colormap contents as cell array
%        contents{get(hObject,'Value')} returns selected item from colormap

val = get(hObject,'Value');
figure(1)
switch val
    case 1
   colormap(jet)
case 2
   colormap(hsv)
case 3
   colormap(hot)
case 4
   colormap(gray)
case 5
    colormap(bone)
case 6
    colormap(copper)
case 7
    colormap(pink)
case 8
    colormap(white)
case 9
    colormap(colorcube)
case 10
    colormap(vga)
case 11
    colormap(jet)
case 12
    colormap(prism)
case 13
    colormap(cool)
case 14
    colormap(autumn)
case 15
    colormap(spring)
case 16
    colormap(winter)
case 17
    colormap(summer)
end

% --- Executes on button press in aspect_ratio.
function aspect_ratio_Callback(hObject, eventdata, handles)
% hObject    handle to aspect_ratio (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

prompt={'X:','Y:','Z:'};
def={'1','1','0.5'};
dlgTitle='Input for Aspect ratio';
lineNo=1;
answer=inputdlg(prompt,dlgTitle,lineNo,def);
M=str2double(answer);
figure(1)
daspect([M(1) M(2) M(3)])
handles.aspect=[M(1) M(2) M(3)];
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function reverse_axes_CreateFcn(hObject, eventdata, handles)
% hObject    handle to reverse_axes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in reverse_axes.
function reverse_axes_Callback(hObject, eventdata, handles)
% hObject    handle to reverse_axes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns reverse_axes contents as cell array
%        contents{get(hObject,'Value')} returns selected item from reverse_axes

val = get(hObject,'Value');
switch val
case 1
    figure(1)
    set(gca,'Xdir','reverse')
case 2
    figure(1)
    set(gca,'Ydir','reverse')
case 3
    figure(1)
    set(gca,'Zdir','reverse')
end

% --- Executes during object creation, after setting all properties.
function display_CreateFcn(hObject, eventdata, handles)
% hObject    handle to display (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function display_Callback(hObject, eventdata, handles)
% hObject    handle to display (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of display as text
%        returns contents of display as a double

step=str2double(get(hObject,'String'));
s=handles.s;
figure(1)
[m,n,p]=size(s);
hz=slice(s,[],[],[1:step:p]);
alpha('color')
axis tight
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
alphamap('rampdown')
aspect=handles.aspect;
daspect(aspect)
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function sliderx_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sliderx (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function sliderx_Callback(hObject, eventdata, handles)
% hObject    handle to sliderx (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

s=handles.s;
[m,n,p]=size(s);
step=1/n;
slider_step(1)=step;
slider_step(2)=step;
set(handles.sliderx, 'SliderStep', slider_step, 'Max', n, 'Min',0)
ix=get(hObject,'Value');
ix=round(ix);
set(handles.xi,'string',ix);
set(handles.ix,'string',ix);
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hx=slice(x,y,z, s,ix,[],[]);
alpha('color')
set(hx,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
handles.IX=ix;
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function slidery_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slidery (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function slidery_Callback(hObject, eventdata, handles)
% hObject    handle to slidery (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
s=handles.s;
[m,n,p]=size(s);
step=1/m;
slider_step(1)=step;
slider_step(2)=step;
set(handles.slidery, 'SliderStep', slider_step, 'Max', m, 'Min',0)
iy=get(hObject,'Value');
iy=round(iy);
set(handles.yi,'string',iy);
set(handles.iy,'string',iy);
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hy=slice(x,y,z, s,[],iy,[]);
alpha('color')
set(hy,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
handles.IY=iy;
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function sliderz_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sliderz (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function sliderz_Callback(hObject, eventdata, handles)
% hObject    handle to sliderz (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

s=handles.s;
[m,n,p]=size(s);
step=1/p;
slider_step(1)=step;
slider_step(2)=step;
set(handles.sliderz, 'SliderStep', slider_step, 'Max', p, 'Min',0)
iz=get(hObject,'Value');
iz=round(iz);
set(handles.zi,'string',iz);
set(handles.iz,'string',iz);
figure(1)
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hz=slice(x,y,z, s,[],[],iz);
alpha('color')
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
handles.IZ=iz;
guidata(hObject, handles);

% --- Executes on button press in all_three.
function all_three_Callback(hObject, eventdata, handles)
% hObject    handle to all_three (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

s=handles.s;
ix=handles.IX;
iy=handles.IY;
iz=handles.IZ;
[m,n,p]=size(s);
figure(1)
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hz=slice(x,y,z, s,ix,iy,iz);
alpha('color')
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
aspect=handles.aspect;
daspect(aspect)
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function ix_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ix (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ix_Callback(hObject, eventdata, handles)
% hObject    handle to ix (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ix as text
%        str2double(get(hObject,'String')) returns contents of ix as a double
ix=str2double(get(hObject,'String')) ;
handles.IX=ix;
s=handles.s;
iy=handles.IY;
iz=handles.IZ;
[m,n,p]=size(s);
figure(1)
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hz=slice(x,y,z, s,ix,iy,iz);
alpha('color')
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
aspect=handles.aspect;
daspect(aspect)
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function iy_CreateFcn(hObject, eventdata, handles)
% hObject    handle to iy (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function iy_Callback(hObject, eventdata, handles)
% hObject    handle to iy (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of iy as text
%        str2double(get(hObject,'String')) returns contents of iy as a double

iy=str2double(get(hObject,'String')) ;
handles.IY=iy;
s=handles.s;
ix=handles.IX;
iz=handles.IZ;
[m,n,p]=size(s);
figure(1)
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hz=slice(x,y,z, s,ix,iy,iz);
alpha('color')
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
aspect=handles.aspect;
daspect(aspect)
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function iz_CreateFcn(hObject, eventdata, handles)
% hObject    handle to iz (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function iz_Callback(hObject, eventdata, handles)
% hObject    handle to iz (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of iz as text
%        str2double(get(hObject,'String')) returns contents of iz as a double

iz=str2double(get(hObject,'String')) ;
handles.IZ=iz;
s=handles.s;
ix=handles.IX;
iy=handles.IY;
[m,n,p]=size(s);
figure(1)
[x,y,z] = meshgrid([1:n],[1:m],[1:p]);
figure(1)
hz=slice(x,y,z, s,ix,iy,iz);
alpha('color')
set(hz,'EdgeColor','none','FaceColor','interp', 'FaceAlpha','interp')
axis tight
alphamap('rampdown')
aspect=handles.aspect;
daspect(aspect)
guidata(hObject, handles);

% --------------------------------------------------------------------
function smooth_Callback(hObject, eventdata, handles)
% hObject    handle to smooth (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

figure(1)
clf
s=handles.s;
size=inputdlg('Enter the size of smoothing kernel','Smoothing filter');
size=str2double(size);
s_sm = smooth3(s, 'gaussian',[size size size]);
handles.s=s_sm;
guidata(hObject, handles);

% --- Executes on slider movement.
function isosurface_Callback(hObject, eventdata, handles)
% hObject    handle to isosurface (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
s=handles.s;
aspect=handles.aspect;
Min=round(min(min(min(s))));
Max=round(max(max(max(s))));
v=[Min:1:Max];
K=Max-Min;
[m,n]=size(v);
step=1/n;
slider_step(1)=step;
slider_step(2)=step;
set(handles.isosurface, 'SliderStep', slider_step, 'Max', K, 'Min',0)
i=get(hObject,'Value');
i=round(i);
if i<=Min
    i=Min;
elseif i>=Max
    i=Max
else
end
i=i+Min;
handles.i=i;
set(handles.iso_value,'string',i);
[n,m,p]=size(s);
a=100*ones([n m p]);
s_inv=a-s;
figure(1)
clf
hiso=patch(isosurface(s,i),'FaceColor',[1,0.75,0.65], 'EdgeColor', 'none','FaceAlpha',0.7);
hcap=patch(isocaps(s_inv,(100-i)),'FaceColor',[1,0.75,0.65], 'Edgecolor','none','FaceAlpha',0.7);
view(45,30);
axis tight
daspect(aspect)
lightangle(45,30)
lighting phong
set(hcap,'AmbientStrength', 0.6)
set(hiso,'SpecularColorReflectance',0.3,'SpecularExponent', 50)
handles.count=0;
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function isosurface_CreateFcn(hObject, eventdata, handles)
% hObject    handle to isosurface (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes during object creation, after setting all properties.
function isosurf_multiple_CreateFcn(hObject, eventdata, handles)
% hObject    handle to isosurf_multiple (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function isosurf_multiple_Callback(hObject, eventdata, handles)
% hObject    handle to isosurf_multiple (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
count=handles.count;
s=handles.s;
value=handles.v;
aspect=handles.aspect;
Min=round(min(min(min(s))));
Max=round(max(max(max(s))));
v=[Min:1:Max];
K=Max-Min;
[m,n]=size(v);
step=1/n;
slider_step(1)=step;
slider_step(2)=step;
set(handles.isosurf_multiple, 'SliderStep', slider_step, 'Max', K, 'Min',0)
i=get(hObject,'Value');
i=round(i);
if i<=Min
    i=Min;
elseif i>=Max
    i=Max
else
    i=i;
end
i=i+Min;
handles.i=i;
set(handles.multiple_iso_value,'string',i);
[n,m,p]=size(s);
figure(1)
hold on
hiso=patch(isosurface(s,i),'FaceColor',[(0.85-count*0.1),(0.25+count*0.1),(0.35+count*0.1)], 'EdgeColor', 'none','FaceAlpha',0.7);
view(45,30);
axis tight
daspect(aspect)
set(hiso,'SpecularColorReflectance',0.3,'SpecularExponent', 50)
count=count+1;
handles.count=count;
value{count}=num2str(i);
outstring=textwrap(value,3);
set(handles.value,'string',outstring)
handles.v=value;
guidata(hObject, handles);

% --- Executes on button press in Clear.
function Clear_Callback(hObject, eventdata, handles)
% hObject    handle to Clear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure(1)
clf
handles.count=0;
handles.v={'0'};
set(handles.value,'string',{})
guidata(hObject, handles);

% --------------------------------------------------------------------
function about_Callback(hObject, eventdata, handles)
% hObject    handle to about (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
msgbox('This a GUI for visualizing 3D data. Created by K.Artyushkova. kartyush@unm.edu. December 2004','About GUI Volumization')

% --- Executes during object creation, after setting all properties.
function iso_value_multi_CreateFcn(hObject, eventdata, handles)
% hObject    handle to iso_value_multi (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

count=handles.count;
s=handles.s;
value=handles.v;
aspect=handles.aspect;
Min=round(min(min(min(s))));
Max=round(max(max(max(s))));
v=[Min:1:Max];
K=Max-Min;
[m,n]=size(v);
i=str2double(get(hObject,'String'));
handles.i=i;
[n,m,p]=size(s);
figure(1)
hold on
hiso=patch(isosurface(s,i),'FaceColor',[(0.85-count*0.1),(0.25+count*0.1),(0.35+count*0.1)], 'EdgeColor', 'none','FaceAlpha',0.7);
view(45,30);
axis tight
daspect(aspect)
set(hiso,'SpecularColorReflectance',0.3,'SpecularExponent', 50)
count=count+1;
handles.count=count;
value{count}=num2str(i);
outstring=textwrap(value,3);
set(handles.value,'string',outstring)
handles.v=value;
guidata(hObject, handles);

function multiple_iso_value_Callback(hObject, eventdata, handles)
% hObject    handle to multiple_iso_value (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of multiple_iso_value as text
%        str2double(get(hObject,'String')) returns contents of multiple_iso_value as a double
count=handles.count;
s=handles.s;
value=handles.v;
aspect=handles.aspect;
Min=round(min(min(min(s))));
Max=round(max(max(max(s))));
v=[Min:1:Max];
K=Max-Min;
[m,n]=size(v);
i=str2double(get(hObject,'String'));
handles.i=i;
[n,m,p]=size(s);
figure(1)
hold on
hiso=patch(isosurface(s,i),'FaceColor',[(0.85-count*0.1),(0.25+count*0.1),(0.35+count*0.1)], 'EdgeColor', 'none','FaceAlpha',0.7);
view(45,30);
axis tight
daspect(aspect)
set(hiso,'SpecularColorReflectance',0.3,'SpecularExponent', 50)
count=count+1;
handles.count=count;
value{count}=num2str(i);
outstring=textwrap(value,3);
set(handles.value,'string',outstring)
handles.v=value;
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function multiple_iso_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to multiple_iso_value (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Sable

repo: mcbench-benchmarks


text:
%%%  HeightControl_Sim
clear
path('./icon/',path);
Init;

% constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
% throttle when UAV is hovering
THR_HOVER = 0.609;
%% initial condition
ModelInit_PosE = [0, 0, -100];
ModelInit_VelB = [0, 0, 0];
ModelInit_AngEuler = [0, 0, 0];
ModelInit_RateB = [0, 0, 0];
ModelInit_Rads = 0;
%% control parameter
% attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
% position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
% integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

% max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
% max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
% maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
% throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
HeightControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function weighted_rmse = mohsst5_rmse_sparsepart(Yh, data)

Itest = load('ind20test.mat');
Itest = Itest.Itest;
if nargin < 2 || isempty(data)
  data = mohsst5_loaddata;
end

% Test set
Ytest = data.observations;
Ytest(~Itest) = nan;
missing = sum(isnan(Ytest(:))) / numel(Ytest);

% Take only data from EARLY parts
rows = size(data.observations,1);
rminds = data.time > datenum([1900 1 1 0 0 0]); % remove by time period
rminds = rowsum(~isnan(data.observations)) > 0.3 * rows; % remove by sparsity
rminds = data.time > datenum([9999 1 1 0 0 0]); % full data
Ytest(:,rminds) = nan;
missing = sum(isnan(Ytest(:))) / numel(Ytest);
testsize = sum(~isnan(Ytest(:)));

% $$$ % Take only data from sparse parts
% $$$ rows = size(data.observations,1);
% $$$ rminds = rowsum(~isnan(data.observations)) > 0.7 * rows;
% $$$ Ytest(:,rminds) = nan;
% $$$ missing = sum(isnan(Ytest(:))) / numel(Ytest)
% $$$ testsize = sum(~isnan(Ytest(:)))

weights = sqrt(cosd(data.coordinates(2,:)));
weights = weights(:);
Weights = bsxfun(@times, ~isnan(Ytest), weights(:)); % mask with missing values
Yhw = bsxfun(@times, Yh, weights(:));
Ytestw = bsxfun(@times, Ytest, weights(:));
err = Yhw(:) - Ytestw(:);
err(isnan(err)) = []; % remove missing values
weighted_rmse = sqrt(err'*err / sum(Weights(:).^2));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function rz = vecrotz(angle)
%
% Function Name:
%
%   vecrotz - Transformation matrix for a rotation around the z axis.
%
% Calling Sequence:
%
%   rz = vecrotz(angle);
%
% Parameters:
%
%   angle	: rotation angle defined in radians
%
%   rz		: (4x4) Transformation matrix.
%
%
% Description:
%
%   Return the (4x4) Transformation matrix for a rotation about the z axis
%   by the defined angle.
%
%   The matrix is:
%
%         [  cos(angle)  -sin(angle)       0          0]
%         [ -sin(angle)   cos(angle)       0          0]
%         [      0            0            1          0]
%         [      0            0            0          1]
%
% Examples:
%
%  Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees
%  around the z-axis
%
%    line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]);
%    trans = vecrotz(%pi/4);
%    rline = nrbtform(line, trans);
%
% See:
%
%    nrbtform

%  Dr D.M. Spink
%  Copyright (c) 2000.

sn = sin(angle);
cn = cos(angle);
rz = [cn -sn 0 0; sn cn 0 0; 0 0 1 0; 0 0 0 1];


text:
function []=panel4fevddisp(n,N,Units,endo,fevd_estimates,IRFperiods,pref)

if pref.plot
    % plot the figure
    fevd=figure('Tag','BEARresults');
    set(fevd,'Color',[0.9 0.9 0.9]);
    set(fevd,'name','forecast error variance decomposition');
    % initiate the count
    count=0;
    % loop over units
    for ii=1:N
        % loop over endogenous variables
        for jj=1:n
            % loop over shocks
            for kk=1:n
                % increment count
                count=count+1;
                % then plot
                subplot(N*n,n,count)
                temp=fevd_estimates{jj,kk,ii};
                hold on
                Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
                Ypatch=[temp(1,:) fliplr(temp(3,:))];
                FEVDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
                set(FEVDpatch,'facealpha',0.5);
                set(FEVDpatch,'edgecolor','none');
                plot(temp(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
                hold off
                set(gca,'XLim',[1 IRFperiods],'YLim',[0 1],'FontName','Times New Roman');
                % top labels
                if jj==1
                    title([Units{ii,1} '\_' endo{kk,1}],'FontWeight','normal');
                end
                % side labels
                if kk==1
                    ylabel([Units{ii,1} '\_' endo{jj,1}],'FontWeight','normal');
                end
            end
        end
    end
    % top supertitle
    ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
    set(get(ax,'Title'),'Visible','on')
    title('Contribution of shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
    % side supertitle
    ylabel('Variance of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
    set(get(ax,'Ylabel'),'Visible','on');
end

% save on Excel
% create the cell that will be saved on excel
fevdcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},IRFperiods+3,1);
horzspace=repmat({''},2,5*n);
% loop over units
for ii=1:N
    % loop over endogenous variables
    for jj=1:n
        % initiate the cell of results
        endocell={};
        % loop over shocks
        for kk=1:n
            % create a header
            header=[{['part of ' Units{ii,1} '_' endo{jj,1} ' fluctuation due to ' Units{ii,1} '_' endo{kk,1} ' shocks']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
            % complete the cell
            tempcell=[[header;num2cell((1:IRFperiods)') num2cell((fevd_estimates{jj,kk,ii})')] vertspace];
            % concatenate to the previous parts of unitcell
            endocell=[endocell tempcell];
        end
        % concatenate to the previous parts of irfcell
        fevdcell=[fevdcell;horzspace;endocell];
    end
    fevdcell=[fevdcell;horzspace];
end
% trim
fevdcell=fevdcell(3:end-2,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),fevdcell,'FEVD','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
clear all;
%% Reader config
%config.Filename = "D:\Dropbox\University of Liverpool\PhD\Workspace\Track Analytics\NelsonRadarAisTracking\data\ais_2017-01-23_2017-01-28.csv";
config.Filename = "C:\Users\sglvladi\Documents\GitHub\Stone-Soup-TA\examples\track_analytics\data\exact_earth\exactEarth_historical_data_2017-08-10.csv";
config.StateFields = {'Longitude', 'Latitude'};
config.TimeField = 'Time';
config.TimeFormat = 'yyyyMMdd_HHmmss';
%config.TimeFormat = 'yyyy-MM-dd''T''HH:mm:ssX';
config.TimeZone = 'UTC';

%% Initialise reader
reader = CsvMeasurementReaderX(config);

%% Read from reader
measurements = reader.read();

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX