text:
function vw = recenter3view(vw,orientation)
% view = recenter3view(view,orientation);
%
% recenters a 3-view window based on the current location of the
% mouse pointer with respect to the anatomy windows. This is set
% as a button-down callback function with respect to each anatomy view:
% so when the user presses at a particular location in a 3-view window,
% it will pick that location and redraw the 3-view centered there. (If the
% pointer happens to be outside the window or there is no 3-view window
% open, this function exits quietly).
%
% orientation is an integer from 1-3 that specifies which view
% the user has clicked on.
%       1 -- view is axial
%       2 -- view is coronal
%       3 -- view is sagittal
%
% See also: volume3view.m
%
% 3/06/03 by ras
% 01/04 ras: (chagrin) turns out the whole 1st half of this code, which was
% a convoluted way of figuring out where the user pressed, is completely
% superfluous if you just use the axes' 'CurrentPoint' property.

% % check if window exists
% winTag = ['3VolumeWindow: ',view.name];
% winExists = findobj('Tag',winTag);
%
% % if window doesn't exist, exit quietly
% if isempty(winExists)
%     return;
% end

pts = get(gca,'CurrentPoint');
locX = round(pts(1,1));
locY = round(pts(1,2));

% figure out R/L, A/P, S/I coordinates, based on view
switch orientation
    case 1, % axial view
        locRL = locX;
        locAP = locY;
        locSI = vw.loc(1);
    case 2, % coronal view
        locRL = locX;
        locAP = vw.loc(2);
        locSI = locY;
    case 3, % sagittal view
        locRL = vw.loc(3);
        locAP = locX;
        locSI = locY;
end

% might have imgs L/R flipped into radiological conventions
if (orientation<3) && isfield(vw.ui,'flipLR') && vw.ui.flipLR==1
    locRL = size(vw.anat,3) - locRL + 1;
end

% make the selected orientation the current one
setCurSliceOri(vw,orientation);

% mark the current location
loc = [locSI locAP locRL];

% if we're zoomed in, update the zoom bounds so that the crosshairs are
% centered:
maxZoom = [1 1 1; viewGet(vw,'Size')]';
if ~isequal(vw.ui.zoom, maxZoom)
	rng = diff(vw.ui.zoom, 1, 2) ./ 2;
	newZoom = [loc(:)-rng loc(:)+rng];

	vw.ui.zoom = newZoom;
end

% refresh volume 3-view
vw = volume3View(vw, loc);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function angle = triangle_angles_3d ( t )

%*****************************************************************************80
%
%% TRIANGLE_ANGLES_3D computes the angles of a triangle in 3D.
%
%  Discussion:
%
%    The law of cosines is used:
%
%      C * C = A * A + B * B - 2 * A * B * COS ( GAMMA )
%
%    where GAMMA is the angle opposite side C.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(DIM_NUM,3), the triangle vertices.
%
%    Output, real ANGLE(3,1), the angles opposite
%    sides P1-P2, P2-P3 and P3-P1, in radians.
%

%
%  Compute the length of each side.
%
  a = sqrt ( sum ( ( t(1:3,1) - t(1:3,2) ).^2 ) );
  b = sqrt ( sum ( ( t(1:3,2) - t(1:3,3) ).^2 ) );
  c = sqrt ( sum ( ( t(1:3,3) - t(1:3,1) ).^2 ) );

  angle = zeros ( 3, 1 );
%
%  Take care of a special case.
%
  if ( a == 0.0 & b == 0.0 & c == 0.0 )
    angle(1:3,1) = 2.0 * pi / 3.0;
    return
  end

  if ( c == 0.0 | a == 0.0 )
    angle(1,1) = pi;
  else
    angle(1,1) = r8_acos ( ( c * c + a * a - b * b ) / ( 2.0 * c * a ) );
  end

  if ( a == 0.0 | b == 0.0 )
    angle(2,1) = pi;
  else
    angle(2,1) = r8_acos ( ( a * a + b * b - c * c ) / ( 2.0 * a * b ) );
  end

  if ( b == 0.0 | c == 0.0 )
    angle(3,1) = pi;
  else
    angle(3,1) = r8_acos ( ( b * b + c * c - a * a ) / ( 2.0 * b * c ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function array = cell2array(cellarray, nD)

% The reason for not taking nD as ndims(cellarray{1}) is that an array
% intended to be 3D can look 2D because the size of the 3rd dimension is 1.
chkarg(istypesizeof(nD, 'int') && nD > 0, '"nD" should be positive integer.');
chkarg(istypesizeof(cellarray, 'complexcell', [1 0], zeros(1, nD)), ...
	'"cellarray" should be a row cell array whose each element is %dD array with complex elements.', nD);

ncell = length(cellarray);
if ncell == 0
	array = [];
else
	assert(ncell >= 1);
	dims = size(cellarray{1});
	chkarg(istypesizeof(cellarray, 'complexcell', [1 0], dims), 'elements of "cellarray" should be have the same size.');
	array = cat(nD+1, cellarray{:});
	array = permute(array, [nD+1, 1:nD]);
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%__________________________________________________________________________
% violin.m - Simple violin plot using matlab default kernel density estimation
% Last update: 10/2015
%__________________________________________________________________________
% This function creates violin plots based on kernel density estimation
% using ksdensity with default settings. Please be careful when comparing pdfs
% estimated with different bandwidth!
%
% Differently to other boxplot functions, you may specify the x-position.
% This is usefule when overlaying with other data / plots.
%__________________________________________________________________________
%
% Please cite this function as:
% Hoffmann H, 2015: violin.m - Simple violin plot using matlab default kernel
% density estimation. INRES (University of Bonn), Katzenburgweg 5, 53115 Germany.
% hhoffmann@uni-bonn.de
%
%__________________________________________________________________________
%
% INPUT
%
% Y:     Data to be plotted, being either
%        a) n x m matrix. A 'violin' is plotted for each column m, OR
%        b) 1 x m Cellarry with elements being numerical colums of nx1 length.
%
% varargin:
% xlabel:    xlabel. Set either [] or in the form {'txt1','txt2','txt3',...}
% facecolor: FaceColor. (default [1 0.5 0]); Specify abbrev. or m x 3 matrix (e.g. [1 0 0])
% edgecolor: LineColor. (default 'k'); Specify abbrev. (e.g. 'k' for black); set either [],'' or 'none' if the mean should not be plotted
% facealpha: Alpha value (transparency). default: 0.5
% mc:        Color of the bars indicating the mean. (default 'k'); set either [],'' or 'none' if the mean should not be plotted
% medc:      Color of the bars indicating the median. (default 'r'); set either [],'' or 'none' if the mean should not be plotted
% bw:        Kernel bandwidth. (default []); prescribe if wanted as follows:
%            a) if bw is a single number, bw will be applied to all
%            columns or cells
%            b) if bw is an array of 1xm or mx1, bw(i) will be applied to cell or column (i).
%            c) if bw is empty (default []), the optimal bandwidth for
%            gaussian kernel is used (see Matlab documentation for
%            ksdensity()
%
% OUTPUT
%
% h:     figure handle
% L:     Legend handle
% MX:    Means of groups
% MED:   Medians of groups
% bw:    bandwidth of kernel
%__________________________________________________________________________
%{
% Example1 (default):

disp('this example uses the statistical toolbox')
Y=[rand(1000,1),gamrnd(1,2,1000,1),normrnd(10,2,1000,1),gamrnd(10,0.1,1000,1)];
[h,L,MX,MED]=violin(Y);
ylabel('\Delta [yesno^{-2}]','FontSize',14)

%Example2 (specify facecolor, edgecolor, xlabel):

disp('this example uses the statistical toolbox')
Y=[rand(1000,1),gamrnd(1,2,1000,1),normrnd(10,2,1000,1),gamrnd(10,0.1,1000,1)];
violin(Y,'xlabel',{'a','b','c','d'},'facecolor',[1 1 0;0 1 0;.3 .3 .3;0 0.3 0.1],'edgecolor','b',...
'bw',0.3,...
'mc','k',...
'medc','r--')
ylabel('\Delta [yesno^{-2}]','FontSize',14)

%Example3 (specify x axis location):

disp('this example uses the statistical toolbox')
Y=[rand(1000,1),gamrnd(1,2,1000,1),normrnd(10,2,1000,1),gamrnd(10,0.1,1000,1)];
violin(Y,'x',[-1 .7 3.4 8.8],'facecolor',[1 1 0;0 1 0;.3 .3 .3;0 0.3 0.1],'edgecolor','none',...
'bw',0.3,'mc','k','medc','r-.')
axis([-2 10 -0.5 20])
ylabel('\Delta [yesno^{-2}]','FontSize',14)

%Example4 (Give data as cells with different n):

disp('this example uses the statistical toolbox')

Y{:,1}=rand(10,1);
Y{:,2}=rand(1000,1);
violin(Y,'facecolor',[1 1 0;0 1 0;.3 .3 .3;0 0.3 0.1],'edgecolor','none','bw',0.1,'mc','k','medc','r-.')
ylabel('\Delta [yesno^{-2}]','FontSize',14)
%}
%%
function[h,L,MX,MED,bw]=violin(Y,varargin)

%defaults:
%_____________________
xL=[];
fc=[1 0.5 0];
lc='k';
alp=0.5;
mc='k';
medc='r';
b=[]; %bandwidth
plotlegend=1;
plotmean=1;
plotmedian=1;
x = [];
%_____________________

%convert single columns to cells:
if iscell(Y)==0
    Y = num2cell(Y,1);
end

%get additional input parameters (varargin)
if isempty(find(strcmp(varargin,'xlabel')))==0
    xL = varargin{find(strcmp(varargin,'xlabel'))+1};
end
if isempty(find(strcmp(varargin,'facecolor')))==0
    fc = varargin{find(strcmp(varargin,'facecolor'))+1};
end
if isempty(find(strcmp(varargin,'edgecolor')))==0
    lc = varargin{find(strcmp(varargin,'edgecolor'))+1};
end
if isempty(find(strcmp(varargin,'facealpha')))==0
    alp = varargin{find(strcmp(varargin,'facealpha'))+1};
end
if isempty(find(strcmp(varargin,'mc')))==0
    if isempty(varargin{find(strcmp(varargin,'mc'))+1})==0
        mc = varargin{find(strcmp(varargin,'mc'))+1};
        plotmean = 1;
    else
        plotmean = 0;
    end
end
if isempty(find(strcmp(varargin,'medc')))==0
    if isempty(varargin{find(strcmp(varargin,'medc'))+1})==0
        medc = varargin{find(strcmp(varargin,'medc'))+1};
        plotmedian = 1;
    else
        plotmedian = 0;
    end
end
if isempty(find(strcmp(varargin,'bw')))==0
    b = varargin{find(strcmp(varargin,'bw'))+1}
    if length(b)==1
        disp(['same bandwidth bw = ',num2str(b),' used for all cols'])
        b=repmat(b,size(Y,2),1);
    elseif length(b)~=size(Y,2)
        warning('length(b)~=size(Y,2)')
        error('please provide only one bandwidth or an array of b with same length as columns in the data set')
    end
end
if isempty(find(strcmp(varargin,'plotlegend')))==0
    plotlegend = varargin{find(strcmp(varargin,'plotlegend'))+1};
end
if isempty(find(strcmp(varargin,'x')))==0
    x = varargin{find(strcmp(varargin,'x'))+1};
end
%%
if size(fc,1)==1
    fc=repmat(fc,size(Y,2),1);
end

%% Calculate the kernel density
i=1;
for i=1:size(Y,2)

    if isempty(b)==0
        [f, u, bb]=ksdensity(Y{i},'bandwidth',b(i));
    elseif isempty(b)
        [f, u, bb]=ksdensity(Y{i});
    end

    f=f/max(f)*0.3; %normalize
    F(:,i)=f;
    U(:,i)=u;
    MED(:,i)=nanmedian(Y{i});
    MX(:,i)=nanmean(Y{i});
    bw(:,i)=bb;

end
%%
%-------------------------------------------------------------------------
% Put the figure automatically on a second monitor
% mp = get(0, 'MonitorPositions');
% set(gcf,'Color','w','Position',[mp(end,1)+50 mp(end,2)+50 800 600])
%-------------------------------------------------------------------------
%Check x-value options
if isempty(x)
    x = zeros(size(Y,2));
    setX = 0;
else
    setX = 1;
    if isempty(xL)==0
        disp('_________________________________________________________________')
        warning('Function is not designed for x-axis specification with string label')
        warning('when providing x, xlabel can be set later anyway')
        error('please provide either x or xlabel. not both.')
    end
end

%% Plot the violins
i=1;
for i=i:size(Y,2)
    if isempty(lc) == 1
        if setX == 0
            h(i)=fill([F(:,i)+i;flipud(i-F(:,i))],[U(:,i);flipud(U(:,i))],fc(i,:),'FaceAlpha',alp,'EdgeColor','none');
        else
            h(i)=fill([F(:,i)+x(i);flipud(x(i)-F(:,i))],[U(:,i);flipud(U(:,i))],fc(i,:),'FaceAlpha',alp,'EdgeColor','none');
        end
    else
        if setX == 0
            h(i)=fill([F(:,i)+i;flipud(i-F(:,i))],[U(:,i);flipud(U(:,i))],fc(i,:),'FaceAlpha',alp,'EdgeColor',lc);
        else
            h(i)=fill([F(:,i)+x(i);flipud(x(i)-F(:,i))],[U(:,i);flipud(U(:,i))],fc(i,:),'FaceAlpha',alp,'EdgeColor',lc);
        end
    end
    hold on
    if setX == 0
        if plotmean == 1
            p(1)=plot([interp1(U(:,i),F(:,i)+i,MX(:,i)), interp1(flipud(U(:,i)),flipud(i-F(:,i)),MX(:,i)) ],[MX(:,i) MX(:,i)],mc,'LineWidth',2);
        end
        if plotmedian == 1
            p(2)=plot([interp1(U(:,i),F(:,i)+i,MED(:,i)), interp1(flipud(U(:,i)),flipud(i-F(:,i)),MED(:,i)) ],[MED(:,i) MED(:,i)],medc,'LineWidth',2);
        end
    elseif setX == 1
        if plotmean == 1
            p(1)=plot([interp1(U(:,i),F(:,i)+i,MX(:,i))+x(i)-i, interp1(flipud(U(:,i)),flipud(i-F(:,i)),MX(:,i))+x(i)-i],[MX(:,i) MX(:,i)],mc,'LineWidth',2);
        end
        if plotmedian == 1
            p(2)=plot([interp1(U(:,i),F(:,i)+i,MED(:,i))+x(i)-i, interp1(flipud(U(:,i)),flipud(i-F(:,i)),MED(:,i))+x(i)-i],[MED(:,i) MED(:,i)],medc,'LineWidth',2);
        end
    end
end

%% Add legend if requested
if plotlegend==1 & plotmean==1 | plotlegend==1 & plotmedian==1

    if plotmean==1 & plotmedian==1
        L=legend([p(1) p(2)],'Mean','Median');
    elseif plotmean==0 & plotmedian==1
        L=legend([p(2)],'Median');
    elseif plotmean==1 & plotmedian==0
        L=legend([p(1)],'Mean');
    end

    set(L,'box','off','FontSize',14)
else
    L=[];
end

%% Set axis
if setX == 0
    axis([0.5 size(Y,2)+0.5, min(U(:)) max(U(:))]);
elseif setX == 1
    axis([min(x)-0.05*range(x) max(x)+0.05*range(x), min(U(:)) max(U(:))]);
end

%% Set x-labels
xL2={''};
i=1;
for i=1:size(xL,2)
    xL2=[xL2,xL{i},{''}];
end
set(gca,'TickLength',[0 0],'FontSize',12)
box on

if isempty(xL)==0
    set(gca,'XtickLabel',xL2)
end
%-------------------------------------------------------------------------
end %of function

author: dynamicslab

repo: SINDy-PI


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

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

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
clear all; close all; clc

% A=randn(10,10,10);
% model=parafac(A,3);
% [A1,A2,A3]=fac2let(model)

x=-5:0.1:5; y=-6:0.1:6; t=0:0.1:10*pi;
[X,Y,T]=meshgrid(x,y,t);
A=exp(-(X.^2+0.5*Y.^2)).*(cos(2*T))+ ...
    (sech(X).*tanh(X).*exp(-0.2*Y.^2)).*sin(T);

for j=1:length(t)
  pcolor(x,y,A(:,:,j)), shading interp, caxis([-1 1]), drawnow
end

figure(1)
for j=1:8
  subplot(2,4,j)
  pcolor(x,y,A(:,:,8*j-3)), colormap(hot), shading interp, caxis([-1 1]), axis off
end

figure(2)
model=parafac(A,2);
[A1,A2,A3]=fac2let(model);
subplot(3,1,1), plot(y,A1,'Linewidth',[2])
subplot(3,1,2), plot(x,A2,'Linewidth',[2])
subplot(3,1,3), plot(t,A3,'Linewidth',[2])

subplot(3,1,1), set(gca,'Xtick',[-6 0 6],'Fontsize',[15])
subplot(3,1,2), set(gca,'Xtick',[-5 0 5],'Fontsize',[15])
subplot(3,1,3), set(gca,'Xlim',[0 10*pi],'Xtick',[0 5*pi 10*pi],'Xticklabels',{'0','5\pi','10\pi'},'Fontsize',[15])


text:
function [d,i0,i1,j0,j1,t0,t1,anglefrom1to2] = dnose2center_pair(trx,fly1,fly2)

% initialize
d = nan(1,trx(fly1).nframes);
anglefrom1to2 = [];

% get start and end frames of overlap
t0 = max(trx(fly1).firstframe,trx(fly2).firstframe);
t1 = min(trx(fly1).endframe,trx(fly2).endframe);

% indices for these frames
i0 = t0 + trx(fly1).off;
i1 = t1 + trx(fly1).off;
j0 = t0 + trx(fly2).off;
j1 = t1 + trx(fly2).off;

% no overlap
if t1 < t0,
  return;
end

% position of nose1
xnose1 = trx(fly1).x_mm(i0:i1) + 2*trx(fly1).a_mm(i0:i1).*cos(trx(fly1).theta_mm(i0:i1));
ynose1 = trx(fly1).y_mm(i0:i1) + 2*trx(fly1).a_mm(i0:i1).*sin(trx(fly1).theta_mm(i0:i1));

% position of center2
x_mm2 = trx(fly2).x_mm(j0:j1);
y_mm2 = trx(fly2).y_mm(j0:j1);

dx = x_mm2-xnose1;
dy = y_mm2-ynose1;

z = sqrt(dx.^2 + dy.^2);
d(i0:i1) = z;

% anglefrom1to2
if nargout >= 8,
  theta1 = trx(fly1).theta_mm(i0:i1);
  theta2 = atan2(dy,dx);
  anglefrom1to2 = modrange(theta2-theta1,-pi,pi);
end

author: kristinbranson

repo: JAABA


text:
classdef TestObjectnessBING
    %TestObjectnessBING

    methods (Static)
        function test_1
            saliency = cv.ObjectnessBING();

            resdir = tempname();
            cObj = onCleanup(@() deleteDir(resdir));
            saliency.setTrainingPath(get_training_path());
            saliency.setBBResDir(resdir);

            img = imread(fullfile(mexopencv.root(),'test','balloon.jpg'));
            objectnessBoundingBox = saliency.computeSaliency(img);
            validateattributes(objectnessBoundingBox, {'cell'}, {'vector'});
            cellfun(@(bb) validateattributes(bb, {'numeric'}, ...
                {'vector', 'integer', 'numel',4}), objectnessBoundingBox);

            objectnessValues = saliency.getObjectnessValues();
            validateattributes(objectnessValues, {'numeric'}, ...
                {'vector', 'numel',numel(objectnessBoundingBox)});
        end
    end

end

function deleteDir(fname)
    if isdir(fname)
        if mexopencv.isOctave()
            %HACK: dont ask for confirmation
            confirm_recursive_rmdir(false, 'local');
        end
        rmdir(fname, 's');
    end
end

function training_path = get_training_path()
    training_path = fullfile(mexopencv.root(),'test','ObjectnessTrainedModel');
    if ~isdir(training_path)
        % download trained models from GitHub
        files = {
            'ObjNessB2W8HSV.idx.yml.gz'
            'ObjNessB2W8HSV.wS1.yml.gz'
            'ObjNessB2W8HSV.wS2.yml.gz'
            'ObjNessB2W8I.idx.yml.gz'
            'ObjNessB2W8I.wS1.yml.gz'
            'ObjNessB2W8I.wS2.yml.gz'
            'ObjNessB2W8MAXBGR.idx.yml.gz'
            'ObjNessB2W8MAXBGR.wS1.yml.gz'
            'ObjNessB2W8MAXBGR.wS2.yml.gz'
        };
        mkdir(training_path);
        for i=1:numel(files)
            url = 'https://cdn.rawgit.com/opencv/opencv_contrib/3.2.0/modules/saliency/samples/ObjectnessTrainedModel/';
            urlwrite([url files{i}], fullfile(training_path,files{i}));
        end
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef PTKAirwayGenerator < handle
    % PTKAirwayGenerator. Creates an artifical airway tree using a volume-filling
    % algorithm.
    %
    %     This class generates an artificial model of the airway tree. You
    %     specify the lung volume to fill, and provide a starting tree (such as
    %     that produced by a CT region-growing algorithm). The airway growing
    %     starts at the endpoints of the starting tree and grows into the
    %     provided volume. The resulting tree includes the starting tree.
    %
    %     You can also choose to grow specific parts of the tree into specific
    %     volumes, so for example a lobar region can be exclusively allocated to
    %     a lobar bronchus and its descendents.
    %
    %     This code has in part been adapted from C++ code by Rafel Bordas which
    %     forms part of the Chaste project at the University of Oxford.
    %     The algorithm is derived from Tawhai et al. (2004), although some
    %     changes have been made to the algorithm.
    %
    %     Syntax:
    %         airway_generator = PTKAirwayGenerator(
    %             lung_mask,              % A binary mask of the whole lung volume
    %             centreline_tree,          % A PTKModelTree produced from PTKAirwayCentreline
    %             point_limit_voxels,     % Branches will terminate if the size of the region they grow into in voxels is less than this limit
    %             approx_number_points,   %
    %             reporting               % A CoreReportingInterface object for error, warning and progress reporting
    %         )
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (Constant)
        LengthLimitMm = 1.2 % A branch is terminated if its length is less than this value. Tawhai et al 2004 use value 2.
        PointLimitVoxels = 1 % A branch is terminated if the point cloud of the apex has fewer points than this limit
        NumberOfGenerationsToReallocate = 20  % If ReallocatePointsAtEachGeneration is set to true, this is the number of generations for which points in the volume will be reallocated at each generation
        PointNumberMultiple = 1 % The desired number of grid points is obtained by multiplying approx_number_points by this value
        MaximumGenerationNumber = 25     % All branches of the output tree will terminate if they extend beyond this generation number
        InitialTerminatingBranchLengthLimit = 6
        BranchingFraction = 0.4 % The fraction a branch extends towards the centre of the point cloud. Value 0.4 from Tawhai et al., 2004
        ReallocatePointsAtEachGeneration = true % If true, points will be reassigned to apices at each generation up to the generation number set in NumberOfGenerationsToReallocate
        AngleLimitRadians = 60*(pi/180.0); % Value 60 degrees from Tawhai et al, 2004
        PointDistanceLimit = 5 % Points greather than a voxel distance of this number multiplied by the number of cloud points are removed from the apex cloud
        BranchLengthToParentRatioLimit = 1 % Child branches cannot be longer than this factor multiplied by ther parent's length
        BranchLengthToParentGenerationLimit = 20 % The BranchLengthToParentRatioLimit parameter starts taking effect from this generation number
    end

    properties
        AirwayTree
        InitialApexImage
        GridSpacingMm
    end

    methods
        function obj = PTKAirwayGenerator(lung_mask, centreline_tree, approx_number_points, reporting)

            % Compute the grid spacing. We choose to have more grid points than
            % the minimum required as a finer grid will give better branching
            approx_number_grid_points = PTKAirwayGenerator.PointNumberMultiple*approx_number_points;
            obj.GridSpacingMm = lung_mask.ComputeResamplingGridSpacing(approx_number_grid_points);

            % Initialise airway tree
            initial_airway_tree = obj.CreateInitialTreeFromSegmentation(centreline_tree, PTKAirwayGenerator.MaximumGenerationNumber, reporting);
            obj.RemoveSmallTerminatingAirways(initial_airway_tree, PTKAirwayGenerator.InitialTerminatingBranchLengthLimit, reporting);
            obj.AirwayTree = initial_airway_tree;
        end

        function delete(~)
        end

        % Starting from the initial airway tree generated from AddTree(),
        function GrowTree(obj, growth_volume, starting_segment, reporting)
            obj.InitialApexImage = PTKAirwayGenerator.GrowTreeUsingThisGridSpacing(obj.AirwayTree, growth_volume, starting_segment, obj.GridSpacingMm, reporting);
        end
    end

    methods (Static, Access = private)
        function initial_apex_image = GrowTreeUsingThisGridSpacing(airway_tree, growth_volume, starting_segment, grid_spacing_mm, reporting)

            reporting.PushProgress;

            resampled_volume = PTKAirwayGenerator.CreatePointCloud(growth_volume, grid_spacing_mm);
            disp(['Number of seed points:' int2str(sum(resampled_volume.RawImage(:)))]);

            initial_apex_image = PTKAirwayGenerator.Grow(resampled_volume, airway_tree, starting_segment, reporting);

            reporting.PopProgress;
        end

        % Use CreateInitialTreeFromSegmentation to create an initial airway tree
        % from the airway centreline results
        function airway_tree = CreateInitialTreeFromSegmentation(segmented_centreline_tree, maximum_generation_number, reporting)
            airway_tree = PTKAirwayGrowingTree;
            airway_tree.CentrelineTreeSegment = segmented_centreline_tree;
            segments_to_do = airway_tree;
            while ~isempty(segments_to_do)
                segment = segments_to_do(end);
                segments_to_do(end) = [];

                % Get the first and last voxel coordinates from the segment of
                % the centreline airway tree
                centreline_segment = segment.CentrelineTreeSegment;
                if isempty(segment.Parent)
                    first_point = centreline_segment.Centreline(1);
                else
                    first_point = centreline_segment.Parent.Centreline(end);
                end
                end_point = centreline_segment.Centreline(end);

                segment.StartPoint = first_point;
                segment.EndPoint = end_point;

                if isempty(centreline_segment.GenerationNumber)
                    error('program error');
                end

                if centreline_segment.GenerationNumber >= maximum_generation_number
                    reporting.ShowWarning('PTKAirwayGenerator:SegmentedBranchesExcluded', 'Initial branches have been excluded due to the maximum generation parameter', []);
                else

                    if ~isempty(centreline_segment.Children)
                        % Add a new branch to the tree for each child
                        for child = centreline_segment.Children
                            % Create a new segment
                            new_segment = PTKAirwayGrowingTree(segment);
                            new_segment.CentrelineTreeSegment = child;
                            new_segment.IsGenerated = false;
                            segments_to_do = [segments_to_do, new_segment];
                        end
                    end
                end

            end
        end

        function RemoveSmallTerminatingAirways(initial_airway_tree, length_limit, reporting)
            branches_to_do = initial_airway_tree;
            has_changed = true;
            iteration_number = 0;
            while (has_changed)
                has_changed = false;
                iteration_number = iteration_number + 1;
                while ~isempty(branches_to_do)
                    branch = branches_to_do(end);
                    branches_to_do(end) = [];
                    if isempty(branch.Children)
                        branch_length = MimImageCoordinateUtilities.DistanceBetweenPoints(branch.StartPoint, branch.EndPoint);
                        if branch_length < length_limit
                            if ~isempty(branch.Parent)
                                reporting.ShowWarning('PTKAirwayGenerator:SegmentedBranchesBelowLimit', 'Initial branches have been excluded due to their length being below the limit', []);
                                branch.Parent.RemoveChildren;
                                has_changed = true;
                            end
                        end

                    else
                        branches_to_do = [branches_to_do, branch.Children];
                    end
                end
            end
        end

        function in_volume = IsInsideVolume(point_mm, lung_volume)
            global_coordinates = MimImageCoordinateUtilities.GetGlobalCoordinatesForPoints(point_mm, lung_volume);
            if ~lung_volume.IsPointInImage(global_coordinates)
                in_volume = false;
            else
                in_volume = lung_volume.GetVoxel(global_coordinates);
            end
        end

        function apices = GetApicesBelowThisBranchForTerminalSegments(centreline_branches, airway_tree, generation_number, reporting)
            apices = [];

            % Find the starting segments
            segments_to_do = [];
            for centreline_branch = centreline_branches
                branch = airway_tree.FindCentrelineBranch(centreline_branch, reporting);
                segments_to_do = [segments_to_do, branch];
            end

            % Now search the tree below the starting segments and create an apex
            % for each one at the correct generation number
            while ~isempty(segments_to_do)
                segment = segments_to_do(end);
                segments_to_do(end) = [];
                children = segment.Children;
                segments_to_do = [segments_to_do, children];

                % Now search the tree below the starting segments and create an apex
                % for each one that is a terminal segment
                if isempty(children) && isempty(segment.IsTerminal)
                    % Add apices for airway growing
                    new_apex = PTKAirwayGeneratorApex(segment, PTKCoords, true);
                    apices = [apices new_apex];
                end
            end
        end

        function apices = GetApicesBelowThisBranchForThisGeneration(start_centreline_branches, airway_tree, generation_number, reporting)
            apices = [];

            % Find the starting segments
            segments_to_do = [];
            for centreline_branch = start_centreline_branches
                branch = airway_tree.FindCentrelineBranch(centreline_branch, reporting);
                segments_to_do = [segments_to_do, branch];
            end

            % Now search the tree below the starting segments and create an apex
            % for each one at the correct generation number
            while ~isempty(segments_to_do)
                segment = segments_to_do(end);
                segments_to_do(end) = [];
                children = segment.Children;
                segments_to_do = [segments_to_do, children];

                % Find all branches of this generation which have not been
                % terminaetd by the algorithm. This includes branches which
                % already have child branches - these will not be grown, but
                % points will be allocated to their apices so that they may
                % later be available for their child branches to grow
                if segment.GenerationNumber == generation_number
                    if isempty(segment.IsTerminal)
                        is_growing_apex = isempty(children);
                        new_apex = PTKAirwayGeneratorApex(segment, PTKCoords, is_growing_apex);
                        apices = [apices new_apex];
                    end
                end
            end
        end

        function generation_number = GetMinimumTerminalGeneration(apices)
            generation_number = [];
            for apex = apices
                apex_generation = apex.AirwayGrowingTreeSegment.GenerationNumber;
                if isempty(generation_number)
                    generation_number = apex_generation;
                else
                    generation_number = min(generation_number, apex_generation);
                end
            end
        end

        % Find the smallest generation number of a branch in the airway tree
        % which is still growing. This will return an empty variable if no more
        % airways are growing
        function generation_number = GetMinimumActiveTerminalGeneration(airway_tree, start_centreline_branches, reporting)
            % Find the starting segments
            segments_to_do = [];
            for centreline_branch = start_centreline_branches
                branch = airway_tree.FindCentrelineBranch(centreline_branch, reporting);
                segments_to_do = [segments_to_do, branch];
            end

            generation_number = [];
            while ~isempty(segments_to_do)
                next_branch = segments_to_do(end);
                child_branches = next_branch.Children;

                % A branch is still growing if it has no child branches and it
                % has not been terminated by the growing algorithm
                if isempty(child_branches) && isempty(next_branch.IsTerminal)

                    % We want the smallest generation number for the tree
                    if (isempty(generation_number) || (generation_number > next_branch.GenerationNumber))
                        generation_number = next_branch.GenerationNumber;
                    end
                end

                segments_to_do(end) = [];
                segments_to_do = [segments_to_do child_branches];
            end
        end

        function centre = GetCentreOfMass(point_cloud)
            cloud_points = point_cloud.Coords;
            centre = mean(cloud_points, 1);
        end

        function new_image = GetImageFromApex(template, apices)
            new_image = zeros(template.ImageSize, 'uint16');
            colour = 1;
            for apex = apices
                cloud_points = apex.PointCloud.Coords;
                if ~isempty(cloud_points)
                    cloud_points = MimImageCoordinateUtilities.PTKCoordinatesToCoordinatesMm(cloud_points);
                    point_coord = template.CoordinatesMmToGlobalCoordinates(cloud_points);
                    point_coord = template.GlobalToLocalCoordinates(point_coord);
                    indices = sub2ind(template.ImageSize, point_coord(:,1), point_coord(:,2), point_coord(:,3));
                    new_image(indices) = colour;
                end
                colour = colour + 1;
            end
        end

        function centreline_indices_local = CentrelinePointsToLocalIndices(centreline_points, template_image)
            centreline_indices_global = MimImageCoordinateUtilities.GetGlobalIndicesForPoints(centreline_points, template_image);
            centreline_indices_local = template_image.GlobalToLocalIndices(centreline_indices_global);
        end

        function closest_point = FindClosestPointInCloud(point, cloud)
            point_coords = [point.CoordX, point.CoordY, point.CoordZ];
            num_points_in_cloud = size(cloud.Coords, 1);
            distance = cloud.Coords - repmat(point_coords, [num_points_in_cloud, 1]);
            distance = sqrt(distance(:,1).^2 + distance(:,2).^2 + distance(:,3).^2);
            [~, min_index] = min(distance, [], 1);
            closest_point = cloud.Coords(min_index, :);
        end

        % Compute a normalised direction vector for a plane parallel to the
        % given vector and the direction of the given branch. If the two vectors
        % are parallel then the direction of the parent branch is used. If this
        % does not exist or is still parallel, then we choose a guaranteed
        % non-parallel vector using the null space.
        function plane_normal = GetValidPlaneNormal(vector_to_com, branch)

            if ~isempty(branch.Parent)
                parent_branch_direction = branch.Parent.Direction;
                plane_normal = cross(vector_to_com, parent_branch_direction);
            else
                branch_direction = branch.Direction;
                plane_normal = cross(vector_to_com, branch_direction);
            end

            % This is an implementation of the algorithm described in Tawhai et
            % al., 2004 for finding the plane. However, this does not work since
            % the parent branch will often be collinear with the centre of mass
            % point (since branches are created in the direction of the centre
            % of mass).
%             branch_direction = branch.Direction;
%
%             % Find the normal to the plane which passes through the centre
%             % of mass and the branch start and end points
%             plane_normal = cross(vector_to_com, branch_direction);
%
%             % If the branch direction is parallel to the vector then we try the
%             % parent direction
%             if isequal(plane_normal, [0 0 0])
%                 disp('parallel - correcting');
%                 if ~isempty(branch.Parent)
%                     parent_branch_direction = branch.Parent.Direction;
%                     plane_normal = cross(vector_to_com, parent_branch_direction);
%                 end
%             end

            % If the vectors are still parallel then we use the null space to
            % find another vector which is perpendicular
            if isequal(plane_normal, [0 0 0])
                disp('still parallel - correcting');
                plane_null_space = null(vector_to_com/norm(vector_to_com));
                perpendicular_vector = plane_null_space(:, 1);
                plane_normal = cross(vector_to_com, perpendicular_vector);
            end

            plane_normal = plane_normal / norm(plane_normal);
        end

        % Divides the point cloud in two, based on a plane perpendicular to the
        % direction of the last two branches
        function [this_plane, other_plane] = SplitPointCloud(point_cloud, normal, origin)
            points_coords = point_cloud.Coords;
            p = - dot(normal, origin);

            % Find which axis is most parallel to the plane
            [~, ordered_directions] = sort(normal, 'descend');
            main_direction = ordered_directions(1);
            other_directions = ordered_directions(2:3);
            plane_points = - (p + normal(other_directions(1))*points_coords(:, other_directions(1)) + normal(other_directions(2))*points_coords(:, other_directions(2)))/normal(main_direction);
            in_plane = points_coords(:, main_direction) < plane_points;
            this_plane = PTKCoords(points_coords(in_plane, :));
            other_plane = PTKCoords(points_coords(~in_plane, :));
        end

        function resampled_volume = CreatePointCloud(growth_volume, grid_spacing_mm)
            resampled_volume = growth_volume.Copy;
            grid_spacing = [grid_spacing_mm, grid_spacing_mm, grid_spacing_mm];
            resampled_volume.Resample(grid_spacing, '*nearest')
        end

        % Checks the branch angle of a proposed growth centre and adjusts it within tolerance, if necessary
        function end_point = CheckBranchAngleLengthAndAdjust(start_point, end_coords, parent_direction, parent_length_mm, generation_number)
            % calculate vector from apex start to the centre
            start_coords = [start_point.CoordX, start_point.CoordY, start_point.CoordZ];
            new_direction = end_coords - start_coords;

            % Record the branching length
            branch_length = norm(new_direction);

            new_direction = new_direction / norm(new_direction);

            old_direction = parent_direction;
            old_direction = old_direction / norm(old_direction);

            % determine branch angle
            dot_product = dot(new_direction, parent_direction);
            branch_angle = acos(dot_product/(norm(new_direction)*norm(parent_direction)));

            % If the branch angle is above the limit then set it to the limit
            if (branch_angle > PTKAirwayGenerator.AngleLimitRadians)
                perpendicular_one = cross(new_direction, old_direction);
                perpendicular_one = perpendicular_one/norm(perpendicular_one);
                perpendicular_two = cross(old_direction, perpendicular_one);
                perpendicular_two = perpendicular_two/norm(perpendicular_two);
                old_direction = old_direction*cos(PTKAirwayGenerator.AngleLimitRadians);
                perpendicular_two = perpendicular_two*sin(PTKAirwayGenerator.AngleLimitRadians);
                new_direction = old_direction + perpendicular_two;
            end

            % Reduce the branch length as required
            branch_length = branch_length*PTKAirwayGenerator.BranchingFraction;
            if generation_number >= PTKAirwayGenerator.BranchLengthToParentGenerationLimit
                branch_length = min(branch_length, parent_length_mm*PTKAirwayGenerator.BranchLengthToParentRatioLimit);
            end
            new_direction = new_direction * branch_length;

            % Update the centre
            end_coords = start_coords + new_direction;
            end_point = PTKPoint(end_coords(1), end_coords(2), end_coords(3));
        end

        function new_apex = AddBranchAndApex(parent_branch, cloud, lung_volume, lung_mask, generation_number, image_size, reporting)
            new_apex = [];

            % Determine the centre of the point cloud
            centre_of_mass = PTKAirwayGenerator.GetCentreOfMass(cloud);

            % The new branch starts at the end of the parent branch
            new_branch_start_point = parent_branch.EndPoint;

            % Fetch the director vector of the parent branch
            parent_direction = parent_branch.Direction;

            % Calculate the end point of the new branch. This is subject both to
            % a length limit and an angle limit.
            new_branch_end_point = PTKAirwayGenerator.CheckBranchAngleLengthAndAdjust(new_branch_start_point, centre_of_mass, parent_direction, parent_branch.Length, generation_number);

            if ~PTKAirwayGenerator.IsInsideVolume(new_branch_end_point, lung_mask)
                reporting.ShowWarning('PTKAirwayGenerator:OutsideVolume', 'Branches terminated because they grew outside the lung volume', []);
                parent_branch.IsTerminal = true;
                return;
            end

            % Create the new branch in AirwayTree
            new_branch = PTKAirwayGrowingTree(parent_branch);
            new_branch.StartPoint = new_branch_start_point;
            new_branch.EndPoint = new_branch_end_point;
            new_branch.IsGenerated = true;

            new_branch_length = MimImageCoordinateUtilities.DistanceBetweenPoints(new_branch_end_point, new_branch_start_point);

            % Find closest point in cloud and delete
            closest_point = PTKAirwayGenerator.FindClosestPointInCloud(new_branch_end_point, cloud);
            closest_point_global_coords = MimImageCoordinateUtilities.PTKCoordinatesToCoordinatesMm(closest_point);
            closest_point_global_coords = lung_volume.CoordinatesMmToGlobalCoordinates(closest_point_global_coords);
            lung_volume.SetVoxelToThis(closest_point_global_coords, 0);

            % Check if the number of points in the cloud is less than the
            % required minimum
            if (size(cloud.Coords, 1) <= PTKAirwayGenerator.PointLimitVoxels)
                new_branch.IsTerminal = true;
                reporting.ShowWarning('PTKAirwayGenerator:BelowLengthThreshold', 'Branches terminated because the number of points was below the threshold', []);
            else
                if (new_branch_length < PTKAirwayGenerator.LengthLimitMm)
                    new_branch.IsTerminal = true;
                    reporting.ShowWarning('PTKAirwayGenerator:BelowLengthThreshold', 'Branches terminated because their length was below the threshold', []);
                else
                    % Create new growth branch
                    new_apex = PTKAirwayGeneratorApex(new_branch, cloud, true);
                end
            end
        end

        function apices = AssignPointCloudToApices(apices, resampled_volume, reporting)
            if numel(apices) == 0
               reporting.Error('PTKAirwayGenerator:AssignPointCloudToApices', 'No starting points');
            end
            sample_image = zeros(resampled_volume.ImageSize, 'uint16');
            apex_start_points = zeros(numel(apices), 3);
            for apex_number = 1 : numel(apices)
                apex = apices(apex_number);
                end_point = apex.AirwayGrowingTreeSegment.EndPoint;
                end_point_global_coordinates = MimImageCoordinateUtilities.GetGlobalCoordinatesForPoints(end_point, resampled_volume);
                end_point = resampled_volume.GlobalToLocalCoordinates(end_point_global_coordinates);
                apex_start_points(apex_number, :) = [end_point(1), end_point(2), end_point(3)];
                current_value = sample_image(end_point(1), end_point(2), end_point(3));
                if current_value ~= 0
                    reporting.ShowWarning('PTKAirwayGenerator:ApexAlreadyUsed', 'The start point for allocating voxels to apices could not be assigned as the point has already been used by another apex', []);
                end
                sample_image(end_point(1), end_point(2), end_point(3)) = apex_number;
            end

            bw_image = sample_image > 0;

            [~, IDX] = bwdist(bw_image);
            mapped_image = sample_image(IDX);
            mapped_image_masked = zeros(size(bw_image), 'uint16');
            point_indices = find(resampled_volume.RawImage);
            mapped_image_masked(point_indices) = mapped_image(point_indices);

            for apex_number = 1 : numel(apices)
                local_indices = find(mapped_image_masked == apex_number);
                if ~isempty(local_indices)
                    [di, dj, dk] = ind2sub(size(mapped_image_masked), local_indices);

                    % Calculate the distance from these new points to the start
                    % point
                    root_point = apex_start_points(apex_number, :);
                    di = di - root_point(1);
                    dj = dj - root_point(2);
                    dk = dk - root_point(3);
                    dist_voxels = sqrt(di.^2 + dj.^2 + dk.^2);
                    points_too_far_away = dist_voxels > PTKAirwayGenerator.PointDistanceLimit*numel(local_indices);

                    % This is a check for distance from root point compared to
                    % parent length
%                     di = (di - root_point(1))*resampled_volume.VoxelSize(1);
%                     dj = (dj - root_point(2))*resampled_volume.VoxelSize(2);
%                     dk = (dk - root_point(3))*resampled_volume.VoxelSize(3);
%                     dist_mm = sqrt(di.^2 + dj.^2 + dk.^2);
%                     last_branch_length = apices(apex_number).AirwayGrowingTreeSegment.Length;
%                     points_too_far_away = dist_mm > PTKAirwayGenerator.PointDistanceLimit*last_branch_length;

                    number_of_points_to_remove = sum(uint8(points_too_far_away));
                    if number_of_points_to_remove > 0
                        disp([int2str(number_of_points_to_remove) ' points removed as they were too far away (generation' int2str(apices(apex_number).AirwayGrowingTreeSegment.GenerationNumber) ')']);
                    end
                    local_indices = local_indices(~points_too_far_away, :);

                    if ~isempty(local_indices)
                        global_indices =  resampled_volume.LocalToGlobalIndices(local_indices);
                        [ic, jc, kc] = resampled_volume.GlobalIndicesToCoordinatesMm(global_indices);
                        [ptk_x, ptk_y, ptk_z] = MimImageCoordinateUtilities.CoordinatesMmToPTKCoordinates(ic, jc, kc);

                        apices(apex_number).PointCloud.Coords = [ptk_x, ptk_y, ptk_z];
                    else
                        apices(apex_number).PointCloud.Coords = [];
                    end
                else
                    apices(apex_number).PointCloud.Coords = [];
                end
            end
        end

        function [apex_1, apex_2] = GrowApex(current_apex, lung_volume, lung_mask, generation_number, image_size, reporting)
            branch = current_apex.AirwayGrowingTreeSegment;

            if isempty(current_apex.PointCloud.Coords)
                branch.IsTerminal = true;
                reporting.ShowWarning('PTKAirwayGenerator:EmptyPointCloud', 'Branches terminated because they had an empty point cloud', []);
                apex_1 = [];
                apex_2 = [];
            elseif size(current_apex.PointCloud.Coords, 1) < 2
                branch.IsTerminal = true;
                reporting.ShowWarning('PTKAirwayGenerator:TooFewPointsInCloud', 'Branches terminated because there was only one point in the point cloud', []);
                apex_1 = [];
                apex_2 = [];
            else

                % Find the end coordinates for this branch
                branch_end_point = branch.EndPoint;
                branch_end_coords = [branch_end_point.CoordX, branch_end_point.CoordY, branch_end_point.CoordZ];

                % Find the centre of mass for the point cloud assigned to this
                % apex
                centre_of_mass = PTKAirwayGenerator.GetCentreOfMass(current_apex.PointCloud);

                % Create a vector from the end branch point to the centre of
                % mass
                vector_to_com = centre_of_mass - branch_end_coords;

                % Get a normalised vector perpendicular to the plane passing
                % through the centre of pass and the start and end points of the
                % branch
                plane_normal = PTKAirwayGenerator.GetValidPlaneNormal(vector_to_com, branch);

                % Split the point cloud into two, creating a new apex for each if there are enough points
                [cloud1, cloud2] = PTKAirwayGenerator.SplitPointCloud(current_apex.PointCloud, plane_normal, branch_end_coords);

                if isempty(cloud1.Coords) || isempty(cloud2.Coords)
                    reporting.ShowWarning('PTKAirwayGenerator:EmptyPointCloud', 'Branches terminated because the point clouds could not be divided into two nonempty clouds', []);
                    branch.IsTerminal = true;
                    apex_1 = [];
                    apex_2 = [];
                    return;
                end

                % Create new branches
                apex_1 = PTKAirwayGenerator.AddBranchAndApex(branch, cloud1, lung_volume, lung_mask, generation_number, image_size, reporting);
                apex_2 = PTKAirwayGenerator.AddBranchAndApex(branch, cloud2, lung_volume, lung_mask, generation_number, image_size, reporting);
            end
        end

        function [apices, initial_apex_image] = Grow(lung_volume, airway_tree, starting_segment, reporting)

            % We need to keep a copy of the lung mask for checking that airways
            % are inside. The original image will have points removed as airways are grown.
            lung_mask = lung_volume.Copy;

            reporting.ShowProgress('Growing branches');
            reporting.UpdateProgressValue(0);
            first_run = true;
            image_size = lung_volume.ImageSize;

            % Set this to something non-empty for the first run
            apices = 1;
            generation_number = 1;

            % Step through generations one by one and terminate when there are
            % no active growing branches left
            while ~isempty(generation_number)

                % Find the lowest generation to grow
                generation_number = PTKAirwayGenerator.GetMinimumActiveTerminalGeneration(airway_tree, starting_segment, reporting);
                if (first_run)
                    min_generation_number = generation_number;
                end

                % Terminate when there are no active growing branches left
                if isempty(generation_number)
                    return;
                end

                % Terminate when we have exceeded a generation threshold, and
                % report how many branches were incomplete
                if generation_number >= PTKAirwayGenerator.MaximumGenerationNumber
                    reporting.ShowMessage('PTKAirwayGenerator:IncompleteApices', [num2str(length(apices)) ' branches were incomplete when the terminal generation was reached']);
                    return
                end

                % Progress reporting
                reporting.ShowMessage('PTKAirwayGenerator:GenerationNumber', ['Generation:' int2str(generation_number)]);
                reporting.UpdateProgressValue(round(100*(generation_number - min_generation_number)/(PTKAirwayGenerator.MaximumGenerationNumber - min_generation_number)));

                % Allocate/reallocate points to the apices
                if (first_run) || ((generation_number < PTKAirwayGenerator.NumberOfGenerationsToReallocate) && (PTKAirwayGenerator.ReallocatePointsAtEachGeneration))

                    is_this_the_last_reallocation = (generation_number + 1 == PTKAirwayGenerator.NumberOfGenerationsToReallocate) || (~PTKAirwayGenerator.ReallocatePointsAtEachGeneration);

                    if is_this_the_last_reallocation
                        apices = PTKAirwayGenerator.GetApicesBelowThisBranchForTerminalSegments(starting_segment, airway_tree, generation_number, reporting);
                    else
                        apices = PTKAirwayGenerator.GetApicesBelowThisBranchForThisGeneration(starting_segment, airway_tree, generation_number, reporting);
                    end

                    num_apices = numel(apices);
                    apices = PTKAirwayGenerator.AssignPointCloudToApices(apices, lung_volume, reporting);
%                     disp(['Generation:' int2str(generation_number) ' Number of apices:' int2str(num_apices) ' Apices post-allocation:' int2str(numel(apices))]);
                    if (first_run)
                        initial_apex_image = PTKAirwayGenerator.GetImageFromApex(lung_volume, apices);
                        first_run = false;
                    end
                end

                apices_next_generation = PTKAirwayGeneratorApex.empty(50000,0);
                for apex = apices

                    % Ignore non-growing apices - they are just there to help
                    % with the point reallocation
                    if (apex.IsGrowingApex)
                        % After the final point reallocation, we may have apices
                        % with higher generation numbers than the current one - this
                        % will happen if the segmented tree has a higher number of
                        % generations than the number of generations to reallocate.
                        % We ignore any apices with a higher generation number than
                        % the current one, but save these for processing later
                        if (generation_number == apex.AirwayGrowingTreeSegment.GenerationNumber)
                            [new_apex_1, new_apex_2] = PTKAirwayGenerator.GrowApex(apex, lung_volume, lung_mask, generation_number, image_size, reporting);
                            if ~isempty(new_apex_1)
                                apices_next_generation(end+1) = new_apex_1;
                            end
                            if ~isempty(new_apex_2)
                                apices_next_generation(end+1) = new_apex_2;
                            end
                        else
                            apices_next_generation(end + 1) = apex;
                        end
                    end
                end

                % The next generation of apices is used when point reallocation
                % no longer happens (either the point reallocation generation
                % has been exceeded, or the reallocation flag is switched off).
                % The list of apices includes newly generated apices, plus any
                % terminal branhces from the segmented airways which are waiting
                % to be grown
                apices = apices_next_generation;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function plotByColorMap(time,velocity,colorVar)
% TODO
%
if ndims(time) ~= ndims(velocity)
    error('Time and Velocity are not same size')
end

% Is there a current figure?
h = get(0,'CurrentFigure');
if isempty(h)
  h = figure;
end

% set colormap environment
colormap('autumn');
cmap = colormap;
[cmapVar, cax] = mcColorMap(colorVar,cmap);
figure(h), set(h,'Colormap',cmap), caxis(cax), colorbar, hold on

% add lines to plot
for i = 1:length(time)
   line(time{i},velocity{i},'color',cmapVar(i,:))
end

% add labels
xlabel('Time (s)')
ylabel('Angular Velocity (rpm)')

author: Sable

repo: mcbench-benchmarks


text:
function dat = calc(e,dat)
   if isa(dat,'cell')
     for i=1:length(dat)
       dat{i}=feval(e.type,e,dat{i});
     end
   else
     dat=feval(e.type,e,dat);
   end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [M,dsym,VoxelSize]=track_clust(TrackFile,Lthresh,J)

% This script computes the distance between each pair of streamlines
% Input
% TrackFile: Streamlines file computed in Diffusion Toolkit
% Lthresh: length threshold of streamlines
% J:Number of points resampled for each streamline
% Output
% M: x y z coordinates for each point on each streamline
% dsym: symmetric distance matrix between pairs of streamlines
% VoxelSize: voxel size in mm

[header,tracks] = trk_read(TrackFile);
VoxelSize=header.voxel_size(1);
fprintf('Total number of streamlines %d\n',length(tracks));

lengths = trk_length(tracks);
ind_ori=find(lengths>=Lthresh);
tracks=tracks(ind_ori);
fprintf('Keep %d streamlines longer than %d mm\n',length(tracks),Lthresh);

% Number of points on each streamline
for i=1:length(tracks)
    nP(i)=tracks(i).nPoints;
end

% Resample so that each streamline comprises exactly J points
fprintf('Resampling to %d points...\n',J);
frst=0;
for i=1:length(tracks)
    x=tracks(i).matrix(:,1); y=tracks(i).matrix(:,2); z=tracks(i).matrix(:,3);
    t=[0;cumsum(sqrt(diff(x).^2+diff(y).^2+diff(z).^2))];
    t=t/t(end);
    ti=linspace(0,1,J);
    xx=spline(t,x,ti);
    yy=spline(t,y,ti);
    zz=spline(t,z,ti);
    tracks_ds(i).matrix=[xx',yy',zz'];
    tracks_ds(i).nPoints=J;
    show_progress(i,length(tracks),frst); frst=1;
end
tracks=tracks_ds; clear track_ds;

% Coordinates(x,y,z)for each point on each streamline
M=zeros(J,3,length(tracks));
for i=1:length(tracks)
    m=tracks(i).matrix;
    M(:,:,i)=m;
end

% Compute the distance across each pair of streamlines
fprintf('Computing distance matrix...\n');  % Lauren et al. 2007
d=zeros(size(M,3),size(M,3)); % distance matrix, asymmetric
frst=0;
for i=1:size(M,3)
    tmp=zeros(J,3,J);
    for k=1:J
        tmp(:,:,k)=repmat(M(k,:,i),J,1);
    end
    for j=1:size(M,3)
        dval=squeeze(sqrt(sum((repmat(M(:,:,j),1,1,J)-tmp).^2,2)));
        d(i,j)=mean(min(dval));
    end
    show_progress(i,size(M,3),frst); frst=1;
end

% Make the distance matrix symmetric
ind_upper=find(triu(ones(size(M,3),size(M,3)),1));
ind_lower=find(tril(ones(size(M,3),size(M,3)),-1));

% Take the mean of the two possible distances (i->j and j->i)
d_avg=mean([d(ind_upper),d(ind_lower)],2);
dsym=zeros(size(M,3),size(M,3));
dsym(ind_upper)=d_avg;
dsym=dsym+dsym';

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function [ y, n ] = r8_upak ( x )

%*****************************************************************************80
%
%% R8_UPAK unpacks an R8 into a mantissa and exponent.
%
%  Discussion:
%
%    This function unpacks a floating point number x so that
%
%      x = y * 2.0^n
%
%    where
%
%      0.5 <= abs ( y ) < 1.0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real X, the number to be unpacked.
%
%    Output, real Y, the mantissa.
%
%    Output, integer N, the exponent.
%
  absx = abs ( x );
  n = 0;
  y = 0.0;

  if ( x == 0.0 )
    return
  end

  while ( absx < 0.5 )
    n = n - 1;
    absx = absx * 2.0;
  end

  while ( 1.0 <= absx )
    n = n + 1;
    absx = absx * 0.5;
  end

  if ( x < 0.0 )
    y = - absx;
  else
    y = + absx;
  end

  return
end


text:
function [P,S] = computeSMFromWProjective( W, method, isCalibrated )
% Compute SfM from measurements for projective camera
%
% This function should not be called directly: it does not contain
% pre-processing steps (removal of calibration matrix if known) or
% post-procssing steps (affine/metric upgrades) and bundle adjustment
% computeSMFromW should always be used instead.
%
% Essential matrix decomposition
% Reference: HZ2, p259, Result 9.19, and p. 294
% nFrame==2 && isCalibrated
% fast
%
% Normalized 8-point algorithm
% Reference: HZ2, p279 and alg 11.1 p. 282
% nFrame==2 && ~isCalibrated
% fast
%
% Projective Factorization
% Reference: HZ2, p445, Algorithm 18.2
% Sturm and Triggs ECCV 96
% A Factorization Based Algorithm for Multi-Image Projective Structure and
% Motion
% nFrame>2 && method==0
%
% Projective Factorization
% Reference: Iterative Extensions of the Sturm/Triggs Algorithm:
% Convergence and Nonconvergence, from Oliensis, Hartley, PAMI 07
% nFrame>2 && method==Inf
% slower
%
% If there are any missing entries:
% Low-Rank Matrix Fitting Based on Subspace Perturbation Analysis
% with Applications to Structure from Motion
% Hongjun Jia, Aleix M. Martinez, PAMI 08
%
% Returns the structure and camera parameters in an Animation object
%
% USAGE
%   anim = computeSMFromWProjective( W, method )
%
% INPUTS
%  W            - [ 2 x nPoint x nFrame ] 2D projected features
%  method       - [inf] method for performing SFM (see above for details)
%  isCalibrated - flag indicating if the camera is calibrated (only used
%                 for 2 frames)
%
% OUTPUTS
%  P          - [ 3 x 4 x nFrame ] projection matrices
%  S          - [ 3 x nPoint ] structure matrix
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

nFrame = size(W,3); nPoint = size(W,2);

P=zeros(3,4,nFrame); P(:,:,1)=eye(3,4);

WIsnan=isnan(W);
hasAnyNan=any(WIsnan(:));

% Normalized 8-point algorithm
% Reference: HZ2, p279 and alg 11.1 p. 282
if nFrame==2
  % Normalize input data
  [x T]=normalizePoint(W(:,:,1),Inf);
  [xp Tp]=normalizePoint(W(:,:,2),Inf);

  A=[xp([1 1],:).*x; xp(1,:); xp([2 2],:).*x;xp(2,:);x;ones(1,nPoint)]';

  [U,S,V]=svd(A,0); F=reshape(V(:,end),[3,3])';
  [U,S,V]=svd(F,0);
  F=U*diag([S(1,1) S(2,2) 0])*V';

  F=[ Tp; 0 0 1 ]'*F*[ T; 0 0 1 ];

  if ~isCalibrated
    P(:,:,2) = convertPF([],F,true);
    S = computeSFromWM( true, W, P, 'method', 0 );
  else
    % Essential matrix decomposition
    % Reference: HZ2, p259, Result 9.19, and p. 294
    [ U disc V ] = svd(F);
    if det(U)<0; U=-U; end; if det(V)<0; V=-V; end

    % Check which of the 4 possibilities gives a point in front of both
    % cameras
    WW = [ 0 -1 0; 1 0 0; 0 0 1 ];
    maxCount=0;
    for i = 1 : 4
      switch i
        case 1,
          P(:,:,2) = [ U*WW*V' U(:,3) ];
        case 2,
          P(:,:,2) = [ U*WW*V' -U(:,3) ];
        case 3,
          P(:,:,2) = [ U*WW'*V' U(:,3) ];
        case 4,
          P(:,:,2) = [ U*WW'*V' -U(:,3) ];
      end

      C = zeros(3,2); v = C;
      for j = 1 : 2
        % camera center (Reference: HZ2, p158-161)
        M = P(:,1:3,j);
        C(:,j)=-M\P(:,4,j);
        % principal axis
        v(:,j) = det(M)*M(3,:)';
      end
      % compute the image of a point
      STmp = computeSFromWM( true, W, P, 'method', 0);
      maxCountTmp = sum((v(:,1)'*bsxfun(@minus,STmp,C(:,1)))>0) + ...
        sum((v(:,2)'*bsxfun(@minus,STmp,C(:,2)))>0);
      if maxCountTmp>maxCount
        P2 = P(:,:,2); maxCount = maxCountTmp; S=STmp;
      end
    end
    P(:,:,2) = P2;
  end
end

% Projective Factorization
% Reference: Iterative Extensions of the Sturm/Triggs Algorithm:
% Convergence and Nonconvergence, from Oliensis, Hartley, PAMI 07
if nFrame>2 && ismember(method,[ 0 Inf ]) && ~hasAnyNan
  % Normalize coordinates
  [ W T ] = normalizePoint(W,-Inf);

  % Initialize
  lam = ones( 1, nPoint, nFrame );
  WSquaredSum = sum(W.^2,1);
  C0Const = sum(WSquaredSum(:));
  for n=1:50
    % Stage 1
    if method==0
      for k = 1 : 10
        % normalize each column and then each row
        lam=bsxfun(@rdivide, lam, sqrt(sum(lam.^2,2)));
        lam=bsxfun(@rdivide, lam, sqrt(sum(lam.^2,3)));
      end
    end

    % get the best rank 4 approximation
    % Wkm1 is for Wk minus 1
    [Wkm1,Wkm1Hat]=projSturmTriggsWkm1Wkm1Hat(lam,W);
    WWkm1Hat = sum(W.*Wkm1Hat,1);

    if method==Inf
      if n==1
        mu = norm(Wkm1(:)-Wkm1Hat(:))^2/norm(Wkm1(:))^4*1.1;
      else
        C0 = mu*C0Const;
        C1 = mu*sum(WWkm1Hat(:));
        C2 = WWkm1Hat.^2./WSquaredSum; C2 = mu*sum(C2(:));
        C3 = mu*sum(Wkm1Hat(:).^2);
      end
    end

    % Stage 2
    % get the optimal lambdas
    lam=WWkm1Hat./WSquaredSum;

    % Stage 3
    if method==Inf && n>1
      a = roots( [ C0, -(C0^2-2*C1), -(2*C0*C3-C2), ...
        -(4*C1*C3-2*C2*C0), C0*C3^2-2*C2*C3, ...
        2*C1*C3^2-C2^2, C2*C3^2 ] );
      a = real(a(abs(imag(a))<1e-10)); a = a(a>=0);

      if length(a)>1
        zkm1=C3*C0/C2;
        if abs(zkm1-1)>eps
          a = a( sign(zkm1-1)==sign(a/sqrt(C3)-1) );
        else a = sqrt(C2/C0)+1;
        end
      end
      lam = ( a + lam )*sqrt(a/(a^2*C0+2*a*C1+C2));
      %        [a,sqrt(a/(a^2*C0+2*a*C1+C2)), mean(lam(:)), std(lam(:))]
    end
  end

  % De-normalize coordinates
  Wk = bsxfun(@times,lam,W);
  [ U S V ] = svd( reshape(permute(Wk,[1,3,2]),3*nFrame,nPoint),...
    'econ' );
  P = U(:,1:4)*S(1:4,1:4); S = V(:,1:4)';

  PTmp = P;
  S = normalizePoint(S,4);
  P = zeros(3,4,nFrame);
  for i=1:nFrame; P(:,:,i) = T(:,:,i)\PTmp(3*i-2:3*i,:); end
end

% Low-Rank Matrix Fitting Based on Subspace Perturbation Analysis
% with Applications to Structure from Motion
% Hongjun Jia, Aleix M. Martinez, PAMI 08
if nFrame>2 && ismember(method,[ 0 Inf ]) && hasAnyNan
  % Normalize coordinates
  [ q TW ] = normalizePoint(W, -Inf);
  HHat=reshape(permute(q,[1,3,2]),3*nFrame,nPoint);
  TInv=1./permute(sum(q.*q,1),[3,2,1]);
  Sk=HHat;
  for k = 1 : 30
    % Sk is HHat (homogeneous W) but scaled by lambda, initialized to 1
    [PStack,S]=lowRankDecomposition(Sk,4);

    % make sure the columns of PStack are orthonormal
    [PStack,disc,disc]=svd(PStack,'econ');

    P = permute(reshape(PStack,3,nFrame,4),[1,3,2]);
    %sumTot=0;
    for j=1:nPoint
      mask=~isnan(TInv(:,j));
      if all(~mask); continue; end
      Cj=permute(sum(bsxfun(@times,q(:,j,mask),P(:,:,mask)),1),[3,2,1]);
      % Solve the maximum of the Rayleigh quotient
      % it is the highest eigenvector of A*v=lam*(1./TjInv)*v
      % or heighest eigen value of TjInv*A*v=lam*v
      A=Cj*Cj';
      [lamjTmp,val]=eigs(bsxfun(@times,TInv(mask,j),A),1,'lm');
      lamj=zeros(1,1,nFrame); lamj(mask)=lamjTmp;
      %sumTot=sumTot+lamjTmp'*Cj*Cj'*lamjTmp/(lamjTmp'*diag(1./ ...
      % TInv(mask,j))*lamjTmp);
      Sk(:,j)=reshape(bsxfun(@times,q(:,j,:),lamj),3*nFrame,1);
    end
    %sumTot;
  end
  [PStack,S]=lowRankDecomposition(Sk,4);

  PTmp = PStack;
  S = normalizePoint(S,4);
  P = zeros(3,4,nFrame);
  for i=1:nFrame; P(:,:,i) = TW(:,:,i)\PTmp(3*i-2:3*i,:); end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ Wkm1, Wkm1Hat ]=projSturmTriggsWkm1Wkm1Hat(lam,W)
% perform the computation of Wk and its rank 4 approximation
% in Storm Triggs
nFrame=size(W,3); nPoint=size(W,2);
Wkm1 = bsxfun(@times,lam,W);
[ U S V ] = svd( reshape(permute(Wkm1,[1,3,2]),3*nFrame,nPoint), 'econ' );
Wkm1Hat = permute(reshape(U(:,1:4)*S(1:4,1:4)*V(:,1:4)',3,nFrame,...
  nPoint),[1,3,2]);
end

author: vrabaud


text:
function g = gold()
  % GOLD return a gold color
  %
  % g = gold()
  %
  % Outputs:
  %   g  1 by 3 gold rgb values
  %

  g = [255 228  58]/255;

end

author: alecjacobson

repo: gptoolbox


text:
function b = nlfilter_cerr(varargin)
%NLFILTER General sliding-neighborhood operations.
%   B = NLFILTER(A,[M N],FUN) applies the function FUN to each M-by-N
%   sliding block of the grayscale image A. FUN is a function that accepts
%   an M-by-N matrix as input and returns a scalar:
%
%       C = FUN(X)
%
%   FUN must be a FUNCTION_HANDLE.
%
%   C is the output value for the center pixel in the M-by-N block X.
%   NLFILTER calls FUN for each pixel in A. NLFILTER zero pads the M-by-N
%   block at the edges, if necessary.
%
%   B = NLFILTER(A,'indexed',...) processes A as an indexed image, padding
%   with ones if A is of class single or double and zeros if A is of class
%   logical, uint8, or uint16.
%
%   Class Support
%   -------------
%   The input image A can be of any class supported by FUN. The class of B
%   depends on the class of the output from FUN. When A is grayscale, it
%   can be any numeric type or logical. When A is indexed, it can be
%   logical, uint8, uint16, single or double.
%
%   Remarks
%   -------
%   NLFILTER can take a long time to process large images. In some cases,
%   the COLFILT function can perform the same operation much faster.
%
%   Example
%   -------
%   This example produces the same result as calling MEDFILT2 with a 3-by-3
%   neighborhood:
%
%       A = imread('cameraman.tif');
%       fun = @(x) median(x(:));
%       B = nlfilter(A,[3 3],fun);
%       imshow(A), figure, imshow(B)
%
%   See also BLOCKPROC, COLFILT, FUNCTION_HANDLE.

%   Copyright 1993-2009 The MathWorks, Inc.
%   $Revision: 5.20.4.10 $  $Date: 2009/12/28 04:16:40 $

% Obsolete syntax:
%   B = NLFILTER(A,[M N],FUN,P1,P2,...) passes the additional parameters
%   P1,P2,..., to FUN.
%

[a, nhood, fun, params, padval] = parse_inputs(varargin{:});

% Expand A
[ma,na] = size(a);
aa = mkconstarray(class(a), padval, size(a)+nhood-1);
aa(floor((nhood(1)-1)/2)+(1:ma),floor((nhood(2)-1)/2)+(1:na)) = a;

% Find out what output type to make.
rows = 0:(nhood(1)-1);
cols = 0:(nhood(2)-1);
b = mkconstarray(class(feval(fun,aa(1+rows,1+cols),params{:})), 0, size(a));

% Apply fun to each neighborhood of a
f = waitbar(0,'Applying neighborhood operation...');
for i=1:ma,
    for j=1:na,
        x = aa(i+rows,j+cols);
        b(i,j) = feval(fun,x,params{:});
    end
    waitbar(i/ma)
end
close(f)

%%%
%%% Function parse_inputs
%%%
function [a, nhood, fun, params, padval] = parse_inputs(varargin)

blockSizeParamNum = 2;

switch nargin
    case {0,1,2}
        eid = sprintf('Images:%s:tooFewInputs',mfilename);
        msg = 'Too few inputs to NLFILTER';
        error(eid,'%s',msg);
    case 3
        if (strcmp(varargin{2},'indexed'))
            eid = sprintf('Images:%s:tooFewInputsIfIndexedImage',mfilename);
            msg = 'Too few inputs to NLFILTER';
            error(eid,'%s',msg);
        else
            % NLFILTER(A, [M N], 'fun')
            a = varargin{1};
            nhood = varargin{2};
            fun = varargin{3};
            params = cell(0,0);
            padval = 0;
        end

    otherwise
        if (strcmp(varargin{2},'indexed'))
            % NLFILTER(A, 'indexed', [M N], 'fun', P1, ...)
            a = varargin{1};
            nhood = varargin{3};
            fun = varargin{4};
            params = varargin(5:end);
            padval = 1;
            blockSizeParamNum = 3;

        else
            % NLFILTER(A, [M N], 'fun', P1, ...)
            a = varargin{1};
            nhood = varargin{2};
            fun = varargin{3};
            params = varargin(4:end);
            padval = 0;
        end
end

if (isa(a,'logical') || isa(a,'uint8') || isa(a,'uint16'))
    padval = 0;
end

% Validate 2D input image
iptcheckinput(a,{'logical','numeric'},{'2d'},mfilename,'A',1);

% Validate neighborhood
iptcheckinput(nhood,{'numeric'},{'integer','row','positive','nonnegative','nonzero'},mfilename,'[M N]',blockSizeParamNum);
if (numel(nhood) ~= 2)
    eid = sprintf('Images:%s:invalidBlockSize',mfilename);
    msg = 'Invalid block size, expected [M N].';
    error(eid,'%s',msg);
end

fun = fcnchk(fun,length(params));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% the following script calculates and plots
% the skin friction factors for smooth pipes for
% laminar, turbulent and transitonal flows
% for Reynnolds number up to 100000
%  CALLED FUNCTION: ffsmooth
% ---------------------------------------------------------------
% The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

clear
format long g;                  % set the format of the calculations
n=2000; %number of points to be plotted

RE=linspace(600,100000,n);

for i=1:n
   ffsm(i)=ffsmooth(RE(i)); %calculate the friction factors
end

figure(1)
h=get(0,'CurrentFigure');
set(h,'Name','Friction Factor for Smooth Pipes','NumberTitle','off');
orient landscape;
loglog(RE,ffsm)
%plot(RE,ffsm)
set(findobj('Type','line'),'Color','k','LineWidth',1.2);  %set all the lines black, and the line width to 1.2 from 0.5 (default)
grid on;
axis ([600 10^5 0.015 0.1]);
%axis tight;
xlabel('{\bf Reynolds number} {\it ( Re={UD_h}/\nu )}');
ylabel('{\bfFriction factor} {\it f}');
title ('{\bf FRICTION FACTOR FOR SMOOTH PIPES}','FontSize',13);
set(gca,'ytick',[0.015 0.02 0.025 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1]);
text(1.05*10^7,0.032,'\it0.006','FontSize',7);
text(1000,0.0205,'{\it To be used for laminar, transitional and turbulent flows,}','FontSize',8);
text(1000,0.0195,'{\it  in smooth pipes at Reynolds numbers below 100000}','FontSize',8);
text(0.99*10^3,0.06,'{\it Laminar flow}','FontSize',8,'Rotation',-62);
text(1.05*10^4,0.0295,'{\it Turbulent flow}','FontSize',8,'Rotation',-28);
text(2400,0.031,'{\it Flow Transition}','FontSize',8,'Rotation',50);
text(4*10^4,0.013,'MATLAB script by Tibor Balint, April 1998','FontSize',6);
return
% -------------- end of the script ----------------------


text:
function ThrowFcn2(~,~, game, res1, N, M, FinishBtn, SelectBtn)

%% GET OFFENDER/DEFENDER DATA
offender = get(game.OffenseTxt, 'UserData');
defender = get(game.DefenseTxt, 'UserData');

Noff = get(game.board(offender).Text, 'UserData');
Ndef = get(game.board(defender).Text, 'UserData');

offOcc = get(game.board(offender).Patch, 'UserData');
defOcc = get(game.board(defender).Patch, 'UserData');

%% SET THROWBUTTON
set(game.ThrowBtnDef, 'Enable', 'off', 'Callback', [])

%% GET NUMBER OF TROOPS TO DEFEND WITH
while true
    N(2) = str2double(inputdlg('Number of defending troops', 'How many troops?', 1, {'1'}));

    if isempty(N(2)) || isnan(N(2))
        continue
    end
    if N(2) > Ndef
        uiwait(msgbox(sprintf('Maximum number of defending troops is %d', Ndef)))
        continue
    end
    if N(2) <= 0
        uiwait(msgbox('Defend with at least 1 unit'))
        continue
    end

    break
end

%% NUMBER OF DICE TO BE USED
M(2) = N(2) * (N(2) <= 2) + 2 * (N(2) > 2);

%% RESULT
res2 = round(rand(1,M(2)) * 5) + 1;
for i = 4 : M(2) + 3
    j = i - 3;
    imshow(sprintf('dice/%d.jpg', res2(j)), 'Parent', game.dice(i));
end
for i = M(2) + 4 : 5
    imshow('dice/7.jpg', 'Parent', game.dice(i));
end

%% INCREMENT COUNTER
set(FinishBtn, 'UserData', get(FinishBtn, 'UserData') + 1);

%% DETERMINE WINNER
X = sort(res1, 'descend');
Y = sort(res2, 'descend');
res = [X(1:min(M)); Y(1:min(M))];

for i = 1:min(M)
    if res(1,i) > res(2,i)
        set(game.board(defender).Text, 'String', num2str(Ndef - 1), ...
                                       'UserData', Ndef - 1);
        Ndef = Ndef - 1;
        playerdata = get(game.Players(defOcc(1)), 'UserData');
        playerdata(1) = playerdata(1) - 1;
        set(game.Players(defOcc(1)), 'UserData', playerdata)
        set(game.Players(defOcc(1)), 'String', info(defOcc(1), game))
    else
        set(game.board(offender).Text, 'String', num2str(Noff - 1), ...
                                       'UserData', Noff - 1);
        Noff = Noff - 1;
        playerdata = get(game.Players(offOcc(1)), 'UserData');
        playerdata(1) = playerdata(1) - 1;
        set(game.Players(offOcc(1)), 'UserData', playerdata)
        set(game.Players(offOcc(1)), 'String', info(offOcc(1), game))
    end
end

if Ndef == 0 % Invasion!
    uiwait(msgbox(sprintf('Victory! You invaded %s!', game.board(defender).Name)))

    % Update Board
    set(game.board(offender).Patch, 'UserData', [offOcc(1), 0], ...
                                    'FaceColor', game.cmap(offOcc(1), :))
    set(game.board(defender).Patch, 'UserData', [offOcc(1), 0], ...
                                    'FaceColor', game.cmap(offOcc(1), :))
    set(game.board(offender).Text, 'String', num2str(Noff - N(1)), ...
                                   'Color', [0 0 0], ...
                                   'UserData', Noff - N(1))
    set(game.board(defender).Text, 'String', num2str(N(1)), ...
                                   'UserData', N(1))

    % Update Scoreboard
    playerdata = get(game.Players(offOcc(1)), 'UserData');
    playerdata(2) = playerdata(2) + 1;
    set(game.Players(offOcc(1)), 'UserData', playerdata)
    set(game.Players(offOcc(1)), 'String', info(offOcc(1), game))
    playerdata = get(game.Players(defOcc(1)), 'UserData');
    playerdata(2) = playerdata(2) - 1;
    set(game.Players(defOcc(1)), 'UserData', playerdata)
    set(game.Players(defOcc(1)), 'String', info(defOcc(1), game))
end

%% RE-ACTIVATE FINISH/SELECT/OFFENSE BUTTON AND RESET TEXT
set(game.ThrowBtnOff, 'Enable', 'on')
set(game.txtH, 'String', 'Offender may throw the dice')
set(FinishBtn, 'Enable', 'on')
set(SelectBtn, 'Enable', 'on')


text:
function [ n_data, x, y, fxy ] = beta_log_values ( n_data )

%*****************************************************************************80
%
%% BETA_LOG_VALUES returns some values of the logarithm of the Beta function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Log[Beta[x]]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  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, Y, the arguments of the function.
%
%    Output, real FXY, the value of the function.
%
  n_max = 17;

  fxy_vec = [ ...
      0.1609437912434100E+01, ...
      0.9162907318741551E+00, ...
      0.5108256237659907E+00, ...
      0.2231435513142098E+00, ...
      0.1609437912434100E+01, ...
      0.9162907318741551E+00, ...
      0.0000000000000000E+00, ...
     -0.1791759469228055E+01, ...
     -0.3401197381662155E+01, ...
     -0.4941642422609304E+01, ...
     -0.6445719819385578E+01, ...
     -0.3737669618283368E+01, ...
     -0.5123963979403259E+01, ...
     -0.6222576268071369E+01, ...
     -0.7138866999945524E+01, ...
     -0.7927324360309794E+01, ...
     -0.9393661429103221E+01 ];

  x_vec = [ ...
     0.2E+00, ...
     0.4E+00, ...
     0.6E+00, ...
     0.8E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     2.0E+00, ...
     3.0E+00, ...
     4.0E+00, ...
     5.0E+00, ...
     6.0E+00, ...
     6.0E+00, ...
     6.0E+00, ...
     6.0E+00, ...
     6.0E+00, ...
     7.0E+00 ];

  y_vec = [ ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     0.2E+00, ...
     0.4E+00, ...
     1.0E+00, ...
     2.0E+00, ...
     3.0E+00, ...
     4.0E+00, ...
     5.0E+00, ...
     2.0E+00, ...
     3.0E+00, ...
     4.0E+00, ...
     5.0E+00, ...
     6.0E+00, ...
     7.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;
    y = 0.0;
    fxy = 0.0;
  else
    x = x_vec(n_data);
    y = y_vec(n_data);
    fxy = fxy_vec(n_data);
  end

  return
end


text:
function [DataMat, ChannelMat] = in_data_besa(DataFile)
% IN_DATA_BESA: Read BESA EEG files.
%
% USAGE:  OutputData = in_data_besa( DataFile )
%
% INPUT:
%     - DataFile : Full path to a recordings file.
% OUTPUT:
%     - DataMat : Brainstorm data (recordings) structure

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

% Get format
[fPath, fBase, fExt] = bst_fileparts(DataFile);
% Initialize returned structure
DataMat = db_template('DataMat');
DataMat.Comment  = fBase;
DataMat.Device   = 'BESA';
DataMat.DataType = 'recordings';
DataMat.nAvg     = 1;
ChanNames = {};

% Open file
fid = fopen(DataFile, 'r');
if (fid == -1)
    error('Cannot open file.');
end

% Switch according to file format
switch lower(fExt)
    case {'.avr', '.mul'}
        % Read header (first line)
        hdr = fgetl(fid);
        % Split to get all the parameters
        hdr = str_split(hdr, ' =');
        % Multiplexed/averaged files
        if strcmpi(fExt, '.mul')
            nTime     = str2num(hdr{2});
            nChannels = str2num(hdr{4});
            timeStart = str2num(hdr{6}) / 1000;  % Convert to seconds
            timeStep  = str2num(hdr{8}) / 1000;  % Convert to seconds
        else
            nTime     = str2num(hdr{2});
            timeStart = str2num(hdr{4}) / 1000;  % Convert to seconds
            timeStep  = str2num(hdr{6}) / 1000;  % Convert to seconds
            nChannels = 0;
        end
        % Read second line: Either the channel names or the first sensor
        hdr = fgetl(fid);
        % If there are alphabetical characters in the line: sensor names
        if any(ismember(double(hdr), double('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz()[]/\_@')))
            % Replace spaces in channel names with _ (eg. "TP 1" to "TP_1")
            for iDig = 0:9
                hdr = strrep(hdr, sprintf(' %d', iDig), sprintf('_%d', iDig));
            end
            % Try to get the channel names
            ChanNames = str_split(hdr, ' ');
        % Else: It was the first sensor
        else
            % Restart reading the file at the beginning
            fseek(fid, 0, 'bof');
            hdr = fgetl(fid);
        end
        % Read the recordings in multiplexed/averaged mode
        if strcmpi(fExt, '.mul')
            DataMat.F = fscanf(fid, '%f', [nChannels, nTime]);
        else
            DataMat.F = fscanf(fid, '%f', [nTime, Inf])';
        end

    case {'.mux'}
        % Skip three lines
        hdr = fgetl(fid);
        hdr = fgetl(fid);
        hdr = fgetl(fid);
        % Read the recordings, line by line
        allLines = {};
        while 1
            newLine = fgetl(fid);
            if ~ischar(newLine)
                break;
            end
            allLines{end+1} = str2num(newLine);
        end
        % Concatenate everything
        DataMat.F = cat(1, allLines{:})';

        % Ask for time window
        res = java_dialog('input', {'Start time (in miliseconds):', 'Sampling frequency'}, ...
                                    'Time definition (in Hz)', [], {'0','1000'});
        if isempty(res) || (length(str2num(res{1})) ~= 1) || (length(str2num(res{2})) ~= 1)
            DataMat = [];
        else
            timeStart = str2num(res{1}) / 1000;
            timeStep  = 1 / str2num(res{2});
        end
    otherwise
        error(['Unsupported file extension: ' fExt]);
end
% Close file
fclose(fid);

% Rebuild time vector
DataMat.Time = timeStart + (0:size(DataMat.F,2)-1) .* timeStep;
% No bad channels defined in those files: all good
nChannels = size(DataMat.F,1);
DataMat.ChannelFlag = ones(nChannels, 1);

% Try to build a channel file
if ~isempty(ChanNames) && (length(ChanNames) == nChannels)
    % Default channel structure
    ChannelMat = db_template('channelmat');
    ChannelMat.Comment = 'BESA channels';
    ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);
    % For each channel
    for i = 1:nChannels
        if ~isempty(ChanNames{i})
            ChannelMat.Channel(i).Name = ChanNames{i};
        elseif (length(ChannelMat.Channel) > 99)
            ChannelMat.Channel(i).Name = sprintf('E%03d', i);
        else
            ChannelMat.Channel(i).Name = sprintf('E%02d', i);
        end
        ChannelMat.Channel(i).Type    = 'EEG';
        ChannelMat.Channel(i).Loc     = [0; 0; 0];
        ChannelMat.Channel(i).Orient  = [];
        ChannelMat.Channel(i).Weight  = 1;
        ChannelMat.Channel(i).Comment = [];
    end
else
    ChannelMat = [];
end

author: brainstorm-tools


text:
function [ Temp,WVapour,Geopot,Pressure,longrid,latgrid,xx,yy,lon0360_flag] = aps_load_merra(file,merra_model,coeff)
% loading merra data and output the variables for TRAIN
%
%     Copyright (C) 2016  Bekaert David
%     davidbekaert.com
%
%     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 2 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, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% modifications
% DB    10/04/2016      Base this code on aps_load_era.m modular approach
% DB    01/05/2016      Include MERRA2 support
% DB    12/11/2017      Change MERRA to be nc4

% debug figure to test and validate dataloading.
debug_fig = 0;

% missing values in the hdf files of MERRA
missing_value = 999999986991104;

% specify the coeficients in case they are not given
if nargin<3
    coeff.Rv= 461.495;           % [j/kg/K] water vapour
    coeff.Rd= 287.05;            % [j/kg/K] specific gas constants dry air
end
if ~isfield(coeff,'Rd') || ~isfield(coeff,'Rv')
   error('Coefficients Rd and Rv are not spcified')
end

%%
if strcmpi(merra_model,'merra')

    [a,b,file_ext] = fileparts(file);

    if strcmpi(file_ext,'.nc4')
        % open the netcdf
        ncid = netcdf.open(file,'NC_NOWRITE');

        % read netcdf variables and get number of variables
        [numdims,numvars,numglobalatts,unlimdimid] = netcdf.inq(ncid);

         % loading all the variables from the netcdf
         for i = 0:numvars-1
            [varname, xtype, dimids, numatts] = netcdf.inqVar(ncid,i);
            flag = 0;
            for j = 0:numatts - 1
                attname1 = netcdf.inqAttName(ncid,i,j);
                attname2 = netcdf.getAtt(ncid,i,attname1);

                if strcmp('add_offset',attname1)
                    offset = attname2;
                end

                if strcmp('scale_factor',attname1)
                    scale = attname2;
                    flag = 1;
                end
            end
            if flag
                eval([varname '= double(netcdf.getVar(ncid,i))*scale + offset;']);
            else
                eval([varname '= double(netcdf.getVar(ncid,i));']);
            end
            clear varname xtype dimids numatts scale offset
         end

         % close nedtcdf
        netcdf.close(ncid)

        % definition of QV
        qv = QV;

        % definition of geopotential
        % Convert Geopotential Height to geopotential
        % important this is needed to be consitent with the other weather modules
        g0 = 9.80665;
        Geopot = H.*g0;

        % define geocoordinates
        lons = XDim;
        lats = YDim;
        clear XDim YDim
        % define temperature
        Temp = T;
        clear T
        % define pressure
        Plevs = Height;
        clear Height

        % this is the time-stamps
        TIME./60;

    elseif  strcmpi(file_ext,'.hdf')
        %% OLD HDF5 files, seems no longer to work for new MERRA website
        % keep it for being backward compatible with earlier TRAIN version
        hdf_file = hdfinfo(file);
        Temp = hdfread(file,'/t');                      % temperature in [K]
        qv  = hdfread(file,'/qv');                      % specific humidity
        H =  hdfread(file,'/h');
        % Convert Geopotential Height to geopotential
        % important this is needed to be consitent with the other weather modules
        g0 = 9.80665;
        Geopot = H.*g0;                                 % geopotential
        Psurface = hdfread(file,'/ps');                 % surface pressure in [hPa]
        Plevs = hdfread(file,'/levels')';               % Pressure in [hPa]
        lons = hdfread(file,'/longitude')';
        lats = hdfread(file,'/latitude')';
    end
elseif strcmpi(merra_model,'merra2')

    % open the netcdf
    ncid = netcdf.open(file,'NC_NOWRITE');

    % read netcdf variables and get number of variables
    [numdims,numvars,numglobalatts,unlimdimid] = netcdf.inq(ncid);

     % loading all the variables from the netcdf
     for i = 0:numvars-1
        [varname, xtype, dimids, numatts] = netcdf.inqVar(ncid,i);
        flag = 0;
        for j = 0:numatts - 1
            attname1 = netcdf.inqAttName(ncid,i,j);
            attname2 = netcdf.getAtt(ncid,i,attname1);

            if strcmp('add_offset',attname1)
                offset = attname2;
            end

            if strcmp('scale_factor',attname1)
                scale = attname2;
                flag = 1;
            end
        end
        if flag
            eval([varname '= double(netcdf.getVar(ncid,i))*scale + offset;']);
        else
            eval([varname '= double(netcdf.getVar(ncid,i));']);
        end
        clear varname xtype dimids numatts scale offset
     end

     % close nedtcdf
    netcdf.close(ncid)

    % definition of QV
    qv = QV;

    % definition of geopotential
    % Convert Geopotential Height to geopotential
    % important this is needed to be consitent with the other weather modules
    g0 = 9.80665;
    Geopot = H.*g0;

    % define geocoordinates
    lons = lon;
    lats = lat;
    clear lon lat
    % define temperature
    Temp = T;
    clear T
    % define pressure
    Plevs = lev;
    clear lev

    % this is the time-stamps
    time./60;
end
%%
% update the no data values with NaN's
Temp(Temp==missing_value)=NaN;
qv(qv==missing_value)=NaN;
Geopot(Geopot==missing_value)=NaN;

% number of lon lat grid nodes
n_latitude_points = size(lats,1);
n_longitude_points = size(lons,1);

% generate the pressure levels
Pressure = repmat(Plevs,[1,n_latitude_points,n_longitude_points]);
Pressure = permute(Pressure,[2,3,1]);

% change the order of variables such its consistent with the pressure levels
if strcmpi(merra_model,'merra2')
    Temp = permute(Temp,[2,1,3]);
    qv = permute(qv,[2,1,3]);
    Geopot = permute(Geopot,[2,1,3]);

elseif  strcmpi(merra_model,'merra')
    if strcmpi(file_ext,'.nc4')
        Temp = permute(Temp,[2,1,3]);
        qv = permute(qv,[2,1,3]);
        Geopot = permute(Geopot,[2,1,3]);
    elseif strcmpi(file_ext,'.hdf')
        Temp = permute(Temp,[2,3,1]);
        qv =  permute(qv,[2,3,1]);
        Geopot = permute(Geopot,[2,3,1]);
    end
end

% computation of the saturated water vapour e
% use mixing ratio r = E*e/(p-e), where E = Rd/Rv
% and specific humidity qv = r/(1+r)
E = coeff.Rd./coeff.Rv;
WVapour= qv.*Pressure./(E.*(1-qv)+qv);
clear qv

% Get list of points to look at in analysis
[xx,yy] = meshgrid(1:n_longitude_points,1:n_latitude_points);

% generate the lon lat grid
latgrid = repmat(lats,[1,length(Plevs),n_longitude_points]);
latgrid = permute(latgrid,[1,3,2]);
longrid = repmat(lons,[1,length(Plevs),n_latitude_points]);
longrid = permute(longrid,[3,1,2]);

% inform about the organisation of the longitudes
if sum(lons>180)>1
    lon0360_flag = 'y';
else
    lon0360_flag = 'n';
end

% validation plots
if debug_fig ==1
    figure('position',[ 3         628        1402         586]);
    subplot(2,5,1)
    imagesc(Temp(:,:,end))
    colorbar
    axis xy
    axis equal
    axis tight
    title('temp upper atmo')
    subplot(2,5,2)
    imagesc(Temp(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('temp lower atmo')

    subplot(2,5,3)
    imagesc(Pressure(:,:,end))
    colorbar
    axis xy
    axis equal
    axis tight
    title('pressure upper atmo')
    subplot(2,5,4)
    imagesc(Pressure(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('pressure lower atmo')

    subplot(2,5,6)
    imagesc(WVapour(:,:,end))
    colorbar
    axis xy
    axis equal
    axis tight
    title('water vapour upper atmo')
    subplot(2,5,7)
    imagesc(WVapour(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('water vapour lower atmo')

    subplot(2,5,8)
    imagesc(Geopot(:,:,end))
    colorbar
    axis xy
    axis equal
    axis tight
    title('geopotential upper atmo')
    subplot(2,5,9)
    imagesc(Geopot(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('geopotential lower atmo')

    subplot(2,5,5)
    imagesc(latgrid(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('lat')
    subplot(2,5,10)
    imagesc(longrid(:,:,1))
    colorbar
    axis xy
    axis equal
    axis tight
    title('lon')
end

%
%
%
%  ncid = netcdf.open(file,'NC_NOWRITE');
% [numdims,numvars,numglobalatts,unlimdimid] = netcdf.inq(ncid);
% [dimname, dimlen] = netcdf.inqDim(ncid,0);
%
%         for k=1:numvars
%             [dimname, dimlen] = netcdf.inqVar(ncid,k-1);
%             fprintf([num2str(k-1) ' - ' dimname '\n'])
%         end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN