text:
classdef PTKTreeUtilities < handle
    % PTKTreeUtilities. Utility functions related to tree structures
    %
    %
    %     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.
    %

    methods (Static)
        function permutations = GetAllBranchPermutations(start_branches, max_generations)
            permutations = [];
            for branch = start_branches
                if ~isempty(branch.Children) && max_generations > 0
                    new_permutations = PTKTreeUtilities.GetAllBranchPermutations(branch.Children, max_generations - 1);
                else
                    new_permutations = [];
                end
                new_permutations{end + 1} = branch;
                if isempty(permutations)
                    permutations = new_permutations;
                else
                    new_list = [];
                    for current_permutation = permutations
                        for new_permutation = new_permutations;
                            combined = [current_permutation{1}, new_permutation{1}];
                            new_list{end + 1} = combined;
                        end
                    end
                    permutations = new_list;
                end
            end
        end

        % Returns a branch which is the ancestor of both specifed branches
        function ancestor = FindCommonAncestor(branch_1, branch_2)
            ancestor = branch_1;
            while ~ancestor.ContainsBranch(branch_2);
                ancestor = ancestor.Parent;
                if isempty(ancestor)
                    return
                end
            end
        end

        % Adds additional branches to the tree to ensure a minimum number of generations
        function ExpandTree(trees, number_of_generations)
            if (number_of_generations > 0)
                for tree = trees
                    if isempty(tree.Children)
                        new_child_1 = PTKTreeModel(tree);
                        new_child_2 = PTKTreeModel(tree);
                        new_child_1.Radius = 0.7*tree.Radius;
                        new_child_2.Radius = 0.7*tree.Radius;
                    end
                    for child = tree.Children
                        PTKTreeUtilities.ExpandTree(child, number_of_generations - 1);
                    end
                end
            end
        end

        % Computes all possible divisions of the tree defined by start_branches
        % into the number of subtrees specified by number_of_branches, where
        % each subtree is specified by a single root branch.
        % Will return an empty matrix if the tree could not be exactly divided
        % into the requested number of branches
        function new_permutations = GetBranchPermutationsForBranchNumber(start_branches, number_of_generations_to_search, number_of_branches_to_find, reporting)

            % Make a copy of the tree - so we can extend it with artificial
            % branches where necessary
            start_branches_copy = [];
            for start_branch = start_branches
                next_start_branches_copy = PTKTreeModel.SimpleTreeCopy(start_branch);
                start_branches_copy = [start_branches_copy, next_start_branches_copy];
            end

            PTKTreeUtilities.ExpandTree(start_branches_copy, number_of_generations_to_search);

            % Find all combinations of branches
            branch_permutations = PTKTreeUtilities.GetAllBranchPermutations(start_branches_copy, number_of_generations_to_search);

            % Remove all combinations which don't have the right number of branches
            new_permutations = [];
            for index = 1 : length(branch_permutations)
                this_permutation = branch_permutations{index};
                if numel(this_permutation) == number_of_branches_to_find
                    new_permutations{end + 1} = this_permutation;
                end
            end
        end

        function largest_branches = GetLargestBranchesFromPermutations(permutations)

            % Find the minimum branch radius for each branch permutation
            min_radius = [];
            for index = 1 : length(permutations)
                this_permutation = permutations{index};
                [sorted_permutation, sorted_radii] = PTKTreeUtilities.SortBranchesByRadiusValues(this_permutation);

                % Choosing min radius
                min_radius(index) = sorted_radii(1);

                permutations{index} = sorted_permutation;
            end

            % Choose the permutation with the largest value of the minimum radius
            [~, sort_indices] = sort(min_radius);

            largest_branch_index = sort_indices(end);
            largest_branches = permutations{largest_branch_index};
        end

        function pruned_permutations = RemovePermutationsWithFalseBranches(permutations)
            pruned_permutations = [];
            for index = 1 : length(permutations)
                next_permutation = permutations{index};
                source_branches = PTKTreeUtilities.BranchesToSourceBranches(next_permutation);
                source_branches = unique(source_branches);

                if numel(source_branches) == numel(next_permutation)
                    pruned_permutations{end + 1} = next_permutation;
                end
            end

        end

        function source_branches = BranchesToSourceBranches(branches)
            source_branches = PTKTreeModel.empty();
            for index = 1 : length(branches)
                next_branch = branches(index);
                while isempty(next_branch.BranchProperties)
                    next_branch = next_branch.Parent;
                end
                source_branches(end+1) = next_branch.BranchProperties.SourceBranch;
            end
        end

        function [sorted_branches, sorted_radii] = SortBranchesByRadiusValues(branches)
            radius_values = [];
            for segment_index = 1 : length(branches)
                radius_values(segment_index) = branches(segment_index).Radius;
            end
            [~, sort_indices] = sort(radius_values);
            sorted_branches = branches(sort_indices);
            sorted_radii = radius_values(sort_indices);
        end

        function sorted_segments = OrderSegmentsByCentroidI(segments_to_order, template)
            centroids_i = [];
            for i = 1 : length(segments_to_order)
                centroids_i(end + 1) = PTKTreeUtilities.GetICentroid(segments_to_order(i), template);
            end

            [~, sorted_indices] = sort(centroids_i);
            sorted_segments = segments_to_order(sorted_indices);
        end

        function sorted_child_indices = OrderByCentroidI(start, template)
            child_indices = start.Children;
            centroids_i = [];
            for i = 1 : length(child_indices)
                centroids_i(end + 1) = PTKTreeUtilities.GetICentroid(child_indices(i), template);
            end

            [~, sorted_indices] = sort(centroids_i);
            sorted_child_indices = child_indices(sorted_indices);
        end

        function sorted_segments = OrderSegmentsByCentroidDistanceFromDiagonalPlane(segments_to_order, template)
            centroids_dp = [];
            for i = 1 : length(segments_to_order)
                centroids_dp(end + 1) = PTKTreeUtilities.GetDPCentroid(segments_to_order(i), template);
            end

            [~, sorted_indices] = sort(centroids_dp);
            sorted_segments = segments_to_order(sorted_indices);
        end

        function centroid_i = GetICentroid(start, template)
            tree = PTKTreeUtilities.CentrelinePointsToLocalIndices(start.GetCentrelineTree, template);
            centroid = PTKTreeUtilities.GetCentroid(tree, template);
            centroid_i = centroid(1);
        end

        function centroid_dp = GetDPCentroid(start_segments, template)
            tree_points = [];
            for segment = start_segments
                tree_points = [tree_points segment.GetCentrelineTree];
            end
            tree = PTKTreeUtilities.CentrelinePointsToLocalIndices(tree_points, template);
            centroid = PTKTreeUtilities.GetCentroid(tree, template);
            centroid_dp = - centroid(3) - centroid(1);
        end

        function centroid = GetCentroid(indices, template)
            [i, j, k] = ind2sub(template.ImageSize, indices);
            centroid = zeros(1, 3);
            centroid(1) = mean(i);
            centroid(2) = mean(j);
            centroid(3) = mean(k);
        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 k_distance = GetKDistance(branch)
            start_centreline_point = branch.GetCentrelineTree;
            start_centreline_point = start_centreline_point(1);
            start_k = start_centreline_point.CoordZ;
            tree_points = branch.GetCentrelineTree;
            k_coords = [tree_points.CoordZ];
            max_k = max(k_coords);
            min_k = min(k_coords);
            if abs(max_k - start_k) > abs(min_k - start_k)
                k_distance = max_k - start_k;
            else
                k_distance = min_k - start_k;
            end
        end

        function voxels = GetCentrelineVoxelsForTheseBranches(start_branches, template)
            voxels = [];
            for index = 1 : numel(start_branches)
                voxels = cat(2, voxels, PTKTreeUtilities.CentrelinePointsToLocalIndices(start_branches(index).GetCentrelineTree, template)');
            end
        end

        function end_voxel = GetEndmostPoint(start_branches)
            while ~isempty(start_branches)
                final_generation = start_branches;
                start_branches = [];
                for next_branch = start_branches
                    start_branches = [start_branches, next_branch{1}.Children];
                end
            end

            centreline_length = [];
            for index = 1 : numel(final_generation)
                centreline_length(index) = numel(final_generation(index).Centreline);
            end
            [~, max_index] = max(centreline_length);
            end_voxel = final_generation(max_index).Centreline(end);
        end

        function voxels = GetCentrelineVoxelsForTheseBranchesExtended(start_branches, template)
            voxels = [];
            if isempty(start_branches)
                return;
            end

            for index = 1 : numel(start_branches)
                voxels = cat(2, voxels, PTKTreeUtilities.CentrelinePointsToLocalIndices(start_branches(index).GetCentrelineTree, template)');
                parent = start_branches(index).Parent;
                while ~isempty(parent)
                    centreline_indices = PTKTreeUtilities.CentrelinePointsToLocalIndices(parent.Centreline, template)';
                    voxels = cat(2, voxels, centreline_indices);
                    parent = parent.Parent;
                end
            end

            % Add nearest neighbours to the list of voxels, otherwise it is possible for
            % a diagnoally-connected centreline segment to pass through a
            % diagnonally-connected airway segment
            voxels = MimImageCoordinateUtilities.AddNearestNeighbours(voxels, template);
        end

        function PruneTree(tree_root)
            % remove each end branch from the tree
            bronchi_to_do = CoreStack(tree_root);
            while ~bronchi_to_do.IsEmpty
                next_bronchus = bronchi_to_do.Pop;
                children = next_bronchus.Children;

                if isempty(children)
                    next_bronchus.CutFromTree;
                else
                    bronchi_to_do.Push(children);
                end
            end

        end

        function PruneTreeByLength(tree_root, min_length)
            % Recursvely prune each end branch of the tree which is less than the minimum length

            bronchi_to_do = CoreStack(tree_root);
            while ~bronchi_to_do.IsEmpty
                next_bronchus = bronchi_to_do.Pop;
                children = next_bronchus.Children;

                if isempty(children)

                    branch_length = next_bronchus.LengthMm;
                    if branch_length < min_length
                      next_bronchus.CutFromTree;
                      bronchi_to_do.Push(next_bronchus.Parent);
                    end
                else
                    bronchi_to_do.Push(children);
                end
            end

        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function pass = test_multipleOutputs(pref)
% This test ensures that the STRINGPARSER class is doing the correct thing when
% called with multiple outputs.

%% BVPS

% Strings we want to test.
str = { 'u(-1) = 1'
    'u(0)=1, u(2) + x = 3'
    'u = 1,v = 0,w = 3'
    'x + diff(sin(v))'
    'sum(u,0,.5) = 0'
    'fred(sin(x-y),u) = fred(cos(x-z),u)'};

% Store the number of strings
n = numel(str);
anFun = cell(n, 1);
indVarNames = anFun; varNames = anFun; pdeVarNames = anFun;
eigVarNames = anFun;
commaSeparated = zeros(n, 1);

for k = 1:n
    [anFun{k}, indVarNames{k}, varNames{k}, pdeVarNames{k}, eigVarNames{k}, ...
        commaSeparated(k)] = stringParser.str2anon(str{k}, 'bvp');
end

correctAnFun = {
    'feval(u,-1)-1'
    'feval(u,0)-1;feval(u,2)+x-3'
    'u-1;v;w-3'
    'x+diff(sin(v))'
    'sum(u,0,.5)'
    'fred(@(x,y)sin(x-y),u)-fred(@(x,z)cos(x-z),u)'};

correctIndVar = {
    {'', ''}
    {'x', ''}
    {'', ''}
    {'x', ''}
    {'', ''}
    {'x', ''}
    };

correctVarNam = {
    {'u'}
    {'u'}
    {'u'; 'v'; 'w'}
    {'v'}
    {'u'}
    {'u'}
    };

correctCommaSep = [0; 1; 1; 0; 0; 0];

passBVP = zeros(n, 1);
for k = 1:n
    passBVP(k) = ( strcmp(anFun{k}, correctAnFun{k}) && ...
        all( strcmp(indVarNames{k}, correctIndVar{k}) ) && ...
        all( strcmp(varNames{k}, correctVarNam{k}) ) && ...
        isempty(pdeVarNames{k}) && isempty(eigVarNames{k}) && ...
        ( commaSeparated(k) == correctCommaSep (k) ) );
end

%% EIGS

str = {'v''''-lambda*v = 0'
    'u''''+u'' = lam*(u + u'') + x*u'};
n2 = numel(str);
anFun = cell(n2, 1);
indVarNames = anFun; varNames = anFun; pdeVarNames = anFun;
eigVarNames = anFun;
commaSeparated = zeros(n2, 1);

for k = 1:n2
    [anFun{k}, indVarNames{k}, varNames{k}, pdeVarNames{k}, eigVarNames{k}, ...
        commaSeparated(k)] = stringParser.str2anon(str{k}, 'eig');
end

correctAnFun = {{'diff(v,2)'
    'v'}
    {'diff(u,2)+diff(u)-x.*u'
    'u+diff(u)'}};

correctIndVar = {
    {'', ''}
    {'x', ''}
    };

correctVarNam = {
    {'v'}
    {'u'}
    };

correctEigNam = {
    {'lambda'}
    {'lam'}
    };

correctCommaSep = [0; 0];
passEIG = zeros(n2, 1);
for k = 1:n2
    passEIG(k) = ( strcmp(anFun{k}{1}, correctAnFun{k}{1}) && ...
        strcmp(anFun{k}{2}, correctAnFun{k}{2}) && ...
        all( strcmp(indVarNames{k}, correctIndVar{k}) ) && ...
        all( strcmp(varNames{k}, correctVarNam{k}) ) && ...
        isempty(pdeVarNames{k}) && ...
        strcmp(eigVarNames{k}, correctEigNam{k}) && ...
        ( commaSeparated(k) == correctCommaSep (k) ) );
end

%% PDES

str = {'u_t+x*u'''' = u'''};
n3 = numel(str);

anFun = cell(n3, 1);
indVarNames = anFun; varNames = anFun; pdeVarNames = anFun;
eigVarNames = anFun;
commaSeparated = zeros(n2, 1);

for k = 1:n3
    [anFun{k}, indVarNames{k}, varNames{k}, pdeVarNames{k}, eigVarNames{k}, ...
        commaSeparated(k)] = stringParser.str2anon(str{k}, 'pde');
end
%%
correctAnFun = {'-x.*diff(u,2)+diff(u)'};
correctIndVar = {{'x', 't'}};
correctVarNam = {{'u'}};
correctPdeNam = {{'u_t'}};
correctCommaSep = 0;

passPDE = zeros(n3, 1);
for k = 1:n3
    passPDE(k) = ( strcmp(anFun{k}, correctAnFun{k}) && ...
        all( strcmp(indVarNames{k}, correctIndVar{k}) ) && ...
        all( strcmp(varNames{k}, correctVarNam{k}) ) && ...
        strcmp(pdeVarNames{k}, correctPdeNam{k}) && ...
        isempty(eigVarNames{k}) && ...
        ( commaSeparated(k) == correctCommaSep (k) ) );
end

%% Concatenate all types
pass = [passBVP; passEIG; passPDE];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function code = getfourcc(fourcc)

code = [char(bitshift(fourcc,0,8)) char(bitshift(fourcc,-8,8)) char(bitshift(fourcc,-16,8)) char(bitshift(fourcc,-24,8))];

%   $Revision: 1.1.6.2 $  $Date: 2004/03/30 13:06:56 $
%   Copyright 1984-2003 The MathWorks, Inc.

author: kristinbranson

repo: JAABA


text:
classdef GroupClass < TreeNodeClass

    properties % (Access = private)
        version;
        subjs;
    end

    properties % (Access = private)
        outputFilename
        oldDerivedPaths
        derivedPathBidsCompliant
        initsaveflag
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function obj = GroupClass(varargin)
            obj@TreeNodeClass(varargin);

            obj.InitVersion();
            obj.oldDerivedPaths = {obj.path, [obj.path, 'homerOutput']};
            obj.derivedPathBidsCompliant = 'derivatives/homer';
            obj.initsaveflag = false;

            obj.type    = 'group';
            obj.subjs   = SubjClass().empty;

            obj.outputFilename = obj.cfg.GetValue('Output File Name');
            if isempty(obj.outputFilename)
                obj.outputFilename = 'groupResults.mat';
            end

            if nargin==0
                return;
            end
            if nargin>0
                if ischar(varargin{1}) && strcmp(varargin{1},'copy')
                    return;
                elseif isa(varargin{1}, 'GroupClass')
                    obj.Copy(varargin{1});
                    return;
                end
                if isa(varargin{1}, 'FileClass')
                    obj.name = varargin{1}.ExtractNames();
                else
                    obj.name = varargin{1};
                end
            end
            if nargin>1
                obj.iGroup = varargin{2};
            end

            if isempty(obj.name)
                % Derive obj name from the name of the root directory
                curr_dir = pwd;
                k = sort([findstr(curr_dir,'/') findstr(curr_dir,'\')]); %#ok<*FSTR>
                obj.name = curr_dir(k(end)+1:end);
            end
        end

        % ----------------------------------------------------------------------------------
        function InitVersion(obj)
            obj.SetVersion(getVernum('DataTree'));
        end

        % ----------------------------------------------------------------------------------
        function SetVersion(obj, vernum)
            obj.version = vernum;
        end

        % ----------------------------------------------------------------------------------
        function vernum = GetVersion(obj)
            vernum = obj.version;
        end

        % ----------------------------------------------------------------------------------
        function filename = GetFilename(obj)
            filename = obj.outputFilename;
        end

        % ----------------------------------------------------------------------------------
        function res = CompareVersions(obj, obj2)
            v1 = versionstr2num(obj.version);
            v2 = versionstr2num(obj2.version);
            if v1 == v2
                res = 0;
            elseif v1 < v2
                res = 1;
            elseif v1 > v2
                res = -1;
            end

            % TBD: not sure how to handle this. For now just return 0
            res = 0;
        end

        % ----------------------------------------------------------------------------------
        % Groups obj1 and obj2 are considered equivalent if their names
        % are equivalent and their subject sets are equivalent.
        % ----------------------------------------------------------------------------------
        function B = equivalent(obj1, obj2)
            B=1;
            if ~strcmp(obj1.name, obj2.name)
                B=0;
                return;
            end
            for i=1:length(obj1.subjs)
                j = existSubj(obj1, i, obj2);
                if j==0 || (obj1.subjs(i) ~= obj2.subjs(j))
                    B=0;
                    return;
                end
            end
            for i=1:length(obj2.subjs)
                j = existSubj(obj2, i, obj1);
                if j==0 || (obj2.subjs(i) ~= obj1.subjs(j))
                    B=0;
                    return;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2, conditional)
            % Copy GroupClass object obj2 to GroupClass object obj. Conditional option applies
            % only to all the runs under this group. If == 'conditional' ONLY derived data,
            % that is, only from procStream but NOT from acquired data is copied for all the runs.
            %
            % Conversly unconditional copy copies all properties in the runs under this group
            if nargin==3 && strcmp(conditional, 'conditional')
                if obj.Mismatch(obj2)
                    return
                end
                for i = 1:length(obj.subjs)
                    j = obj.existSubj(i,obj2);
                    if j>0
                        obj.subjs(i).Copy(obj2.subjs(j), 'conditional');
                    elseif i<=length(obj2.subjs)
                        obj.subjs(i).Copy(obj2.subjs(i), 'conditional');
                    else
                        obj.subjs(i).Mismatch();
                    end
                end
                obj.Copy@TreeNodeClass(obj2, 'conditional');
            else
                for i=1:length(obj2.subjs)
                    obj.subjs(i) = SubjClass(obj2.subjs(i));
                end
                obj.Copy@TreeNodeClass(obj2);
            end
        end

        % --------------------------------------------------------------
        function CopyStims(obj, obj2)
            obj.CondNames = obj2.CondNames;
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).CopyStims(obj2.subjs(ii));
            end
        end

        % ----------------------------------------------------------------------------------
        function CopyFcalls(obj, varargin)
            if isa(varargin{1},'TreeNodeClass')
                procStream  = varargin{1}.procStream;
                type        = varargin{1}.type;
            elseif isa(varargin{1},'ProcStreamClass')
                procStream  = varargin{1};
                type        = varargin{2};
            end

            % Copy default procStream function call chain to all uninitialized nodes
            % in the group
            switch(type)
                case 'group'
                    obj.procStream.CopyFcalls(procStream);
                case 'subj'
                    for jj = 1:length(obj.subjs)
                        obj.subjs(jj).procStream.CopyFcalls(procStream);
                    end
                case 'sess'
                    for ii = 1:length(obj.subjs)
                        for jj = 1:length(obj.subjs(ii).sess)
                            obj.subjs(ii).sess(jj).procStream.CopyFcalls(procStream);
                        end
                    end
                case 'run'
                    for ii = 1:length(obj.subjs)
                        for jj = 1:length(obj.subjs(ii).sess)
                            for kk = 1:length(obj.subjs(ii).sess(jj).runs)
                                obj.subjs(ii).sess(jj).runs(kk).procStream.CopyFcalls(procStream);
                            end
                        end
                    end
            end
        end

        % ----------------------------------------------------------------------------------
        function Add(obj, subj, sess, run)
            [~,f,e] = fileparts(subj.GetName());
            if strcmp(f, obj.name)
                msg{1} = sprintf('WARNING: The subject being added (%s) has the same name as the group (%s) containing it. ', [f,e], obj.name);
                msg{2} = sprintf('The subject names should not have the same name as the group folder, otherwise ');
                msg{3} = sprintf('it may cause incorrect results in processing.');
                obj.logger.Write('%s\n', [msg{:}]);
            end

            % Add subject to this group
            jj = 0;
            for ii = 1:length(obj.subjs)
                if strcmp(obj.subjs(ii).GetName(), subj.GetName())
                    jj=ii;
                    break;
                end
            end
            if jj==0
                jj = length(obj.subjs)+1;
                subj.SetIndexID(obj.iGroup, jj);
                subj.SetPath(obj.path);                      % Inherit root path from group
                obj.subjs(jj) = subj;
                obj.logger.Write('   Added subject  "%s"  to group  "%s" .\n', obj.subjs(jj).GetFileName, obj.GetFileName);
            end

            % Add sess to subj
            obj.subjs(jj).Add(sess, run);
            obj.children = obj.subjs;
        end

        % ----------------------------------------------------------------------------------
        function list = DepthFirstTraversalList(obj)
            list{1} = obj;
            for ii=1:length(obj.subjs)
                list = [list; obj.subjs(ii).DepthFirstTraversalList()];
            end
        end

        % #########################################################################

        % -----------------------------------------------------------------------------
        function [g, s, e, r] = GetInitialFuncCallChain(obj)
            g = obj;
            s = obj.subjs(1);
            e = obj.subjs(1).sess(1);
            r = obj.subjs(1).sess(1).runs(1);
            for ii = 1:length(obj.subjs)
                if ~obj.subjs(ii).procStream.IsEmpty()
                    s = obj.subjs(ii);
                end
                for jj = 1:length(obj.subjs(ii).sess)
                    if ~obj.subjs(ii).sess(jj).procStream.IsEmpty()
                        e = obj.subjs(ii).sess(jj);
                    end
                    for kk = 1:length(obj.subjs(ii).sess(jj).runs)
                        if ~obj.subjs(ii).sess(jj).runs(kk).procStream.IsEmpty()
                            r = obj.subjs(ii).sess(jj).runs(kk);
                        end
                    end
                end
            end

            % Generate procStream defaults at each level with which to initialize
            % any uninitialized procStream.input
            g.CreateProcStreamDefault();
        end

        % ----------------------------------------------------------------------------------
        function [fnameFull, status] = InitProcStreamLevel(obj, g, o, fnameFull)
            status = 0;

            procStreamDefault = o.GetProcStreamDefault();

            % If any of the tree nodes still have unintialized procStream input, ask
            % user for a config file to load it from
            if o.procStream.IsEmpty()

                if ~ispathvalid(fnameFull)
                    [fnameFull, autoGenDefaultFile] = g.procStream.GetConfigFileName(fnameFull, obj.path);
                else
                    autoGenDefaultFile = false;
                end

                % If user did not provide procStream config filename and file does not exist
                % then create a config file with the default contents
                if ~exist(fnameFull, 'file')
                    procStreamDefault.SaveConfigFile(fnameFull, o.type);
                end

                obj.logger.Write('Attempting to load %s-level proc stream from %s\n', o.type, fnameFull);

                % Load file to the first empty procStream in the dataTree at each processing level
                err = o.LoadProcStreamConfigFile(fnameFull);

                % If proc stream input is still empty it means the loaded config
                % did not have valid proc stream input. If that's the case we
                % Load a default proc stream input
                if err ~= 0
                    [~, fname, ext] = fileparts(fnameFull);
                    msg{1} = sprintf('Some functions at the %s-level failed to load from selected proc stream config file,  "%s".  ', o.type, [fname, ext]);
                    msg{2} = sprintf('The functions may be obsolete or contain errors.  Loading default %s proc stream ...\n', o.type);
                    obj.logger.Write([msg{:}]);
                    g.CopyFcalls(procStreamDefault, o.type);
                    status = 1;

                % Otherwise the non-default processing stream loaded from file to this group and to first subject
                % disseminate it to all subjects and all runs in this group
                else

                    obj.logger.Write('Loading proc stream from %s\n', fnameFull);
                    g.CopyFcalls(o.procStream, o.type);

                end
                obj.initsaveflag = true;

            end
        end

        % ----------------------------------------------------------------------------------
        function ErrorCheckInitErr(obj, procStreamCfgFile, status);
            if ~all(status==0)
                [~, fname, ext] = fileparts(procStreamCfgFile);
                levels = '';
                if status(1)~=0
                    levels = 'group';
                end
                if status(2)~=0
                    if isempty(levels)
                        levels = 'subject';
                    else
                        levels = [levels, ', subject'];
                    end
                end
                if status(3)~=0
                    if isempty(levels)
                        levels = 'session';
                    else
                        levels = [levels, ', session'];
                    end
                end
                if status(4)~=0
                    if isempty(levels)
                        levels = 'run';
                    else
                        levels = [levels, ', run'];
                    end
                end
                k = find(levels == ',');
                if ~isempty(k)
                    levels = sprintf('%s and %s levels', levels(1:k-1), levels(k+2:end));
                    procStreamLabelPlural = 'streams';
                else
                    levels = sprintf('%s level', levels);
                    procStreamLabelPlural = 'stream';
                end
                msg{1} = sprintf('\nWARNING: There were errors loading user functions from "%s" at the %s.\n', [fname, ext], levels);
                msg{2} = sprintf('These functions may have changed since this file was created or have errors in the help\n');
                msg{3} = sprintf('sections of the processing stream description. Replacing %s processing\n', levels);
                msg{4} = sprintf('%s with default processing %s ... \n', procStreamLabelPlural, procStreamLabelPlural);
                % if strcmpi(obj.cfg.GetValue('Include Archived User Functions'), 'No')
                %       MenuBox(msg, {'OK'});
                % end
                obj.logger.Write([msg{:}]);
            end
        end

        % ----------------------------------------------------------------------------------
        function InitProcStream(obj, procStreamCfgFile, options)
            if isempty(obj)
                return;
            end
            if ~exist('procStreamCfgFile','var')
                procStreamCfgFile = '';
            end
            if ~exist('options','var')
                options = '';
            end

            if optionExists(options, 'noloadconfig')
            	return
            end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Find out if we need to ask user for processing options
            % config file to initialize procStream.fcalls at the
            % run, subject or group level. First try to find the proc
            % input at each level from the saved derived data.
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            [g, s, e, r] = obj.GetInitialFuncCallChain();

            %
            [procStreamCfgFile, status(1)] = obj.InitProcStreamLevel(g, g, procStreamCfgFile);
            [procStreamCfgFile, status(2)] = obj.InitProcStreamLevel(g, s, procStreamCfgFile);
            [procStreamCfgFile, status(3)] = obj.InitProcStreamLevel(g, e, procStreamCfgFile);
            [~,                 status(4)] = obj.InitProcStreamLevel(g, r, procStreamCfgFile);

            obj.ErrorCheckInitErr(procStreamCfgFile, status);
        end

        % ---------------------------------------------------------------
        function PrintProcStream(obj)
            fcalls = obj.procStream.GetFuncCallChain();
            obj.logger.Write('Group processing stream:\n');
            for ii = 1:length(fcalls)
                obj.logger.Write('%s\n', fcalls{ii});
            end
            obj.logger.Write('\n');
            obj.subjs(1).PrintProcStream();
        end

       % ----------------------------------------------------------------------------------
        function FreeMemoryRecursive(obj)
            if isempty(obj)
                return
            end
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).FreeMemoryRecursive();
            end
            obj.FreeMemory();
        end

        % ----------------------------------------------------------------------------------
        function LoadRecursive(obj)
            if isempty(obj)
                return
            end
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).LoadRecursive();
            end
            obj.Load();
        end

        % ----------------------------------------------------------------------------------
        function LoadInputVars(obj, tHRF_common)
            obj.inputVars = [];
            for iSubj = 1:length(obj.subjs)
            % Set common tHRF: make sure size of tHRF, dcAvg and dcAvg is same for
            % all subjects. Use smallest tHRF as the common one.
                obj.subjs(iSubj).procStream.output.SettHRFCommon(tHRF_common, obj.subjs(iSubj).name, obj.subjs(iSubj).type);

                obj.inputVars.dodAvgSubjs{obj.subjs(iSubj).iSubj}    = obj.subjs(iSubj).procStream.output.GetVar('dodAvg');
                obj.inputVars.dodAvgStdSubjs{obj.subjs(iSubj).iSubj} = obj.subjs(iSubj).procStream.output.GetVar('dodAvgStd');
                obj.inputVars.dcAvgSubjs{obj.subjs(iSubj).iSubj}     = obj.subjs(iSubj).procStream.output.GetVar('dcAvg');
                obj.inputVars.dcAvgStdSubjs{obj.subjs(iSubj).iSubj}  = obj.subjs(iSubj).procStream.output.GetVar('dcAvgStd');
                obj.inputVars.tHRFSubjs{obj.subjs(iSubj).iSubj}      = obj.subjs(iSubj).procStream.output.GetTHRF();
                obj.inputVars.nTrialsSubjs{obj.subjs(iSubj).iSubj}   = obj.subjs(iSubj).procStream.output.GetVar('nTrials');

                obj.subjs(iSubj).FreeMemory();
        end
        end

        % ----------------------------------------------------------------------------------
        function Calc(obj, options)
            if ~exist('options','var') || isempty(options)
                options = 'overwrite';
            end

            if strcmpi(options, 'overwrite')
                % Recalculating result means deleting old results, if
                % option == 'overwrite'
                obj.procStream.output.Flush();
            end
            if obj.DEBUG
                obj.logger.Write('Calculating processing stream for group %d\n', obj.iGroup)
            end

            % Calculate all subjs in this session
            tHRF_common = {};
            for iSubj = 1:length(obj.subjs)
                obj.subjs(iSubj).Calc();

                % Find smallest tHRF among the subjs and make this the common one.
                tHRF_common = obj.subjs(iSubj).procStream.output.GeneratetHRFCommon(tHRF_common);
            end

            % Update call application GUI using it's generic Update function
            if ~isempty(obj.updateParentGui)
                obj.updateParentGui('DataTreeClass', [obj.iGroup, obj.iSubj, obj.iSess, obj.iRun]);
            end

            % Load all the output valiraibles that might be needed by procStream.Calc() to calculate proc stream for this group
            obj.LoadInputVars(tHRF_common);

            Calc@TreeNodeClass(obj);

            if obj.DEBUG
                obj.logger.Write('Completed processing stream for group %d\n', obj.iGroup);
                obj.logger.Write('\n');
            end

        end

        % ----------------------------------------------------------------------------------
        function Print(obj, indent)
            obj.logger.Write('\n');
            if ~exist('indent', 'var')
                indent = 0;
            end
            Print@TreeNodeClass(obj, indent);
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).Print(indent);
            end
        end

        % ----------------------------------------------------------------------------------
        % Deletes derived data in procStream.output
        % ----------------------------------------------------------------------------------
        function Reset(obj, option)
            if ~exist('option','var')
                option = 'down';
            end
            if exist([obj.path, obj.outputDirname, obj.outputFilename],'file')
                delete([obj.path, obj.outputDirname, obj.outputFilename]);
            end
            if strcmp(option, 'down')
                for jj = 1:length(obj.subjs)
                    obj.subjs(jj).Reset();
                end
            end
            Reset@TreeNodeClass(obj);
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            for ii = 1:length(obj.subjs)
                if ~obj.subjs(ii).IsEmpty()
                    b = false;
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmptyOutput(obj)
            b = true;
            if isempty(obj)
                return;
            end
            for ii = 1:length(obj.subjs)
                if ~obj.subjs(ii).IsEmptyOutput()
                    b = false;
                    break;
                end
            end
        end

    end   % Public methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public Save/Load methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            for ii = 1:length(obj.subjs)
                nbytes = nbytes + obj.subjs(ii).MemoryRequired();
            end
            nbytes = nbytes + obj.procStream.MemoryRequired();
        end

        % ----------------------------------------------------------------------------------
        function err = LoadSubBranch(obj)
            err = -1;
            if isempty(obj)
                return;
            end
            err1 = obj.procStream.Load([obj.path, obj.GetOutputFilename()]);
            err2 = obj.subjs(1).LoadSubBranch();
            if err1==0 && err2==0
                err = 0;
            end
        end

        % ----------------------------------------------------------------------------------
        function FreeMemorySubBranch(obj)
            if isempty(obj)
                return;
            end
            obj.subjs(1).FreeMemorySubBranch()
        end

        % ----------------------------------------------------------------------------------
        function err = Load(obj, options)
            err = -1;
            if isempty(obj)
                return;
            end
            if ~exist('options','var')
                options = '';
            end

            % If this group has been loaded, then no need to go through the whole Load function. Instead
            % default to the generic TreeNodeClass.Load method.
            if ~optionExists(options, {'init','reload'})
                err = obj.Load@TreeNodeClass();
                return;
            end

            obj.BackwardCompatability();

            group = [];
            if ispathvalid([obj.path, obj.outputDirname, obj.outputFilename],'file')
                g = load([obj.path, obj.outputDirname, obj.outputFilename]);

                % Do some basic error checks on saved derived data contents
                if isproperty(g, 'group') && isa(g.group, 'GroupClass')
                    if isproperty(g.group, 'version')
                        if ismethod(g.group, 'GetVersion')
                            % obj.logger.Write('Saved group data, version %s exists\n', g.group.GetVersionStr());
                            group = g.group;
                        end
                    end
                end
            end

            % Copy saved group to current group if versions are compatible. obj.CompareVersions==0
            % means the versions of the saved group and current one are equal.
            if ~isempty(group) && obj.CompareVersions(group) >= 0
                % Do a conditional copy of group from saved processing output file. Conditional copy copies ONLY
                % derived data, that is, only from procStream but NOT acqruired. We do not want to
                % overwrite the real acquired data loaded from acquisition files
                hwait = waitbar(0,'Loading group');
                obj.Copy(group, 'conditional');
                close(hwait);
            else
                obj.initsaveflag = true;
            end
            err = 0;
        end

        % ----------------------------------------------------------------------------------
        function Save(obj, hwait)
            if ~exist('hwait','var')
                hwait = [];
            end
            if obj.initsaveflag == false
                obj.initsaveflag = true;
                return
            end

            obj.logger.Write('Saving processed data in %s\n', [obj.path, obj.outputDirname, obj.outputFilename]);

            if ishandle(hwait)
                obj.logger.Write('Auto-saving processing results ...\n', obj.logger.ProgressBar(), hwait);
            end

            group = GroupClass(obj);
            try
                obj.CreateOutputDir();
                save([obj.pathOutputAlt, obj.outputDirname, obj.outputFilename], 'group');
            catch ME
                MessageBox(ME.message);
                obj.logger.Write(ME.message);
            end
            obj.initsaveflag = true;
        end

        % ----------------------------------------------------------------------------------
        function CreateOutputDir(obj)
            if ispathvalid([obj.pathOutputAlt, obj.outputDirname])
                return;
            end
            mkdir([obj.pathOutputAlt, obj.outputDirname]);
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).CreateOutputDir();
            end
        end

        % ----------------------------------------------------------------------------------
        function SaveAcquiredData(obj)
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).SaveAcquiredData();
            end
        end

        % ----------------------------------------------------------------------------------
        function b = AcquiredDataModified(obj)
            b = false;
            for ii = 1:length(obj.subjs)
                if obj.subjs(ii).AcquiredDataModified()
                    b = true;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function varval = GetVar(obj, varname)
            % First call the common code for all levels
            varval = obj.GetVar@TreeNodeClass(varname);

            % Now call the group specific part
            if isempty(varval)
                varval = obj.subjs(1).GetVar(varname);
            end
        end

    end  % Public Save/Load methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public Set/Get methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function probe = GetProbe(obj)
            probe = obj.subjs(1).GetProbe();
        end

        % ----------------------------------------------------------------------------------
        function SD = GetSDG(obj,option)
            if exist('option','var')
                SD = obj.subjs(1).GetSDG(option);
            else
                SD = obj.subjs(1).GetSDG();
            end
        end

        % ----------------------------------------------------------------------------------
        function bbox = GetSdgBbox(obj)
            bbox = obj.subjs(1).GetSdgBbox();
        end

        % ----------------------------------------------------------------------------------
        function wls = GetWls(obj)
            wls = obj.subjs(1).GetWls();
        end

        % ----------------------------------------------------------------------------------
        function [iDataBlks, ich] = GetDataBlocksIdxs(obj, ich)
            if nargin<2
                ich = [];
            end
            [iDataBlks, ich] = obj.subjs(1).GetDataBlocksIdxs(ich);
        end

        % ----------------------------------------------------------------------------------
        function n = GetDataBlocksNum(obj)
            n = obj.subjs(1).GetDataBlocksNum();
        end

        % ----------------------------------------------------------------------------------
        function aux = GetAuxiliary(obj) %#ok<MANU>
            aux = [];
        end

    end      % Public Set/Get methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Conditions related methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function SetConditions(obj)
            if isempty(obj)
                return;
            end

            % First get global et of conditions across all runs and
            % subjects
            CondNames = {};
            for ii=1:length(obj.subjs)
                obj.subjs(ii).SetConditions();
                CondNames = [CondNames, obj.subjs(ii).GetConditions()];
            end
            obj.CondNames    = unique(CondNames);

            % Now that we have all conditions, set the conditions across
            % the whole group to these
            for ii=1:length(obj.subjs)
                obj.subjs(ii).SetConditions(obj.CondNames);
            end
        end

        % ----------------------------------------------------------------------------------
        function [fn_error, missing_args, prereqs] = CheckProcStreamOrder(obj)
            missing_args = {};
            fn_error = 0;
            prereqs = '';
            for i = 1:length(obj.subjs)
                [fn_error, missing_args, prereqs] = obj.subjs(i).CheckProcStreamOrder;
                if ~isempty(missing_args)
                    return
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function CondNames = GetConditionsActive(obj)
            CondNames = obj.CondNames;
            for ii=1:length(obj.subjs)
                CondNamesSubj = obj.subjs(ii).GetConditionsActive();
                for jj=1:length(CondNames)
                    k = find(strcmp(['-- ', CondNames{jj}], CondNamesSubj));
                    if ~isempty(k)
                        CondNames{jj} = ['-- ', CondNames{jj}];
                    end
                end
            end
        end

    end

    methods

        % ----------------------------------------------------------------------------------
        function r = ListOutputFilenames(obj, options)
            if ~exist('options','var')
                options = '';
            end
            r = obj.GetOutputFilename(options);
            fprintf('%s %s\n', obj.path, r);
            for ii = 1:length(obj.subjs)
                obj.subjs(ii).ListOutputFilenames(options);
            end
        end

        % ---------------------------------------------------------------
        function CleanUpOutput(obj, filesObsolete)
            for jj = 1:length(filesObsolete)
                renameFlag = false;
                for ii = 1:length(filesObsolete(jj).files)
                    if isempty(filesObsolete(jj).files(ii).namePrev)
                        continue;
                    end
                    renameFlag = true;
                end

                % If something changed in the folder structure
                if renameFlag
                    msg{1} = sprintf('Previous Homer3 processing output exists but is now inconsistent with the current ');
                    msg{2} = sprintf('data files. This output should be regenerated in the new Homer3 session to reflect the new file/folder names. ');
                    msg{3} = sprintf('The existing Homer processing output will be moved to %s. Is this okay?', obj.GetArchivedOutputDirname());
                    q = MenuBox(msg, {'YES','NO'});
                    if q==1
                        if isempty(obj.outputDirname)
                            movefile('*.mat', obj.GetArchivedOutputDirname())
                            movefile('*.txt', obj.GetArchivedOutputDirname())
                        else
                            movefile(obj.outputDirname, obj.GetArchivedOutputDirname())
                        end
                        obj.Save();
                    end
                end
            end
        end

        % -----------------------------------------------------------------
        function name = GetArchivedOutputDirname(obj)
            n = 1;
            addon = '_old';
            if isempty(obj.outputDirname)
                base = 'homerOutput';
            else
                base = filesepStandard(obj.outputDirname, 'file');
            end
            name = sprintf('%s%s%d', base, addon, n);
            while ispathvalid(name)
                n = n+1;
                name = sprintf('%s%s%d', base, addon, n);
            end
        end

        % --------------------------------------------------------------------------
        function ApplyParamEditsToAllSubjects(obj, iFcall, iParam, val)
            for jj = 1:length(obj.subjs)
                obj.subjs(jj).procStream.EditParam(iFcall, iParam, val);
            end
        end

        % --------------------------------------------------------------------------
        function ApplyParamEditsToAllSessions(obj, iFcall, iParam, val)
            for jj = 1:length(obj.subjs)
                obj.subjs(jj).ApplyParamEditsToAllSessions(iFcall, iParam, val);
            end
        end

        % --------------------------------------------------------------------------
        function ApplyParamEditsToAllRuns(obj, iFcall, iParam, val)
            for jj = 1:length(obj.subjs)
                obj.subjs(jj).ApplyParamEditsToAllRuns(iFcall, iParam, val);
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Private methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public)

        % ----------------------------------------------------------------------------------
        % Check whether subject k'th subject from this group exists in group G and return
        % its index in G if it does exist. Else return 0.
        % ----------------------------------------------------------------------------------
        function j = existSubj(obj, k, G)
            j=0;
            for i=1:length(G.subjs)
                if strcmp(obj.subjs(k).name, G.subjs(i).name)
                    j=i;
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function BackwardCompatability(obj)
            for jj = 1:length(obj.oldDerivedPaths)
                oldDerivedPath = filesepStandard(obj.oldDerivedPaths{jj});
                if ispathvalid([oldDerivedPath, 'groupResults.mat'])
                    try
                        g = load([oldDerivedPath, 'groupResults.mat']);
                    catch
                        g = [];
                    end

                    % Do not try to restore old data older than Homer3
                    if isempty(g)
                        return;
                    end
                    if ~isproperty(g,'group')
                        return;
                    end
                    if ~isa(g.group, 'GroupClass')
                        return;
                    end

                    % Do not try to restore old data if there is already data
                    % in the new format
                    if obj.HaveOutput()
                        return;
                    end

                    oldDerivedPathRel = pathsubtract_startup(oldDerivedPath, obj.path);
                    if isempty(oldDerivedPathRel)
                        return
                    end
                    msg{1} = sprintf('Detected derived data in older Homer3 folder "%s" ', oldDerivedPathRel);
                    if pathscompare(obj.derivedPathBidsCompliant, obj.outputDirname, 'nameonly')
                        msg{2} = sprintf('The current derived output folder, "%s", is BIDS compliant. ', ...
                            filesepStandard(obj.derivedPathBidsCompliant, 'filesepwide:nameonly'));
                    else
                        msg{2} = '.';
                    end
                    msg = [msg{:}];
                    obj.logger.Write('Backward Compatability:   %s\n', msg);

                    % If we're here it means that old format homer3 data exists
                    % AND NO new homer3 format data exists
                    q = MenuBox(sprintf('%s Do you want to move %s to the new folder?', msg, oldDerivedPathRel), {'Yes','No'});
                    if q==1
                        if ispathvalid([obj.path, obj.outputDirname])
                            try
                                rmdir([obj.path, obj.outputDirname], 's')
                            catch
                                MenuBox(sprintf('ERROR:  Could not remove new derived folder'), 'OK');
                                return
                            end
                        end
                        obj.logger.Write('Moving %s to %s\n', oldDerivedPath, [obj.path, obj.outputDirname]);
                        if ~ispathvalid([obj.path, obj.outputDirname])
                            mkdir([obj.path, obj.outputDirname]);
                        end

                        try
                            movefile(oldDerivedPath, [obj.path, obj.outputDirname])
                        catch
                        end

                        if ispathvalid([obj.path, obj.outputDirname, obj.outputFilename])
                            if ~strcmp(obj.outputFilename, 'groupResults.mat')
                                obj.logger.Write('Moving %s to %s\n', [obj.path, obj.outputDirname, 'groupResults.mat'], ...
                                    [obj.path, obj.outputDirname, obj.outputFilename])
                                movefile([obj.path, obj.outputDirname, 'groupResults.mat'], [obj.path, obj.outputDirname, obj.outputFilename])
                            end
                        end
                    end
                    break
                end

            end
        end

    end  % Private methods

end % classdef GroupClass < TreeNodeClass

author: BUNPC


text:
function PTKSaveTreeAsNodes(tree_root, file_path, filename_prefix, coordinate_system, template_image, reporting)
    % PTKSaveTreeAsNodes. Exports a tree structure into node and element files
    %
    % Two text files are produced; a node file and an element file.
    %
    % The node file describes points which comprise the start point of the
    % tree (at the top of the trachea), each bifurcation point, and each
    % end point. The first line contains the number of node points.
    % Subequent lines are in the following text format:
    % node_number, x_coord, y_coord, z_coord, radius_approx, is_bifurcation_point
    %
    % Note the radius approximation is for the element ending at this node.
    %
    % The element file describes which nodes are linked together.
    % The first line contains the number of element points.
    % Subsequent lines are in the following text format:
    % node_1, node_2
    % where node_1 and node_2 are connected by the tree
    %
    %     Syntax
    %     ------
    %
    %         PTKSaveTreeAsNodes(tree_root, file_path, base_filename, reporting)
    %
    %             tree_root       is the root branch in a PTKTreeModel structure
    %             file_path       is the path where the node and element files
    %                             are to be stored
    %             filename_prefix is the filename prefix. The node and element
    %                             files will have '_node.txt' and '_element.txt'
    %                             appended to this prefix before saving.
    %             coordinate_system  a MimCoordinateSystem enumeration
    %                             specifying the coordinate system to use
    %             template_image  A PTKImage providing voxel size and image size
    %                             parameters
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     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.
    %

    if nargin < 4
        reporting.Error('PTKSaveTreeAsNodes:BadArguments', 'No coordinate_system parameter specified');
    end

    if ~isa(coordinate_system, 'MimCoordinateSystem')
        reporting.Error('PTKSaveTreeAsNodes:BadArguments', 'coordinate_system parameter is not of type MimCoordinateSystem');
    end

    reporting.ShowProgress('Saving tree');
    reporting.UpdateProgressValue(0);

    node_file_name = fullfile(file_path, [filename_prefix '_node.txt']);
    element_file_name = fullfile(file_path, [filename_prefix '_element.txt']);

    linear_branch_list = tree_root.GetBranchesAsList;

    num_branches = length(linear_branch_list);

    % Create the node file (contains the coordinates and radius)
    node_file_handle = fopen(node_file_name, 'w');
    if node_file_handle == -1
        reporting.Error('PTKSaveTreeAsNodes:CreateFileFailed', ['I could not create the file ' node_file_name '. Please check the disk permissions to ensure I have write access.']);
    end

    element_file_handle = fopen(element_file_name, 'w');
    if element_file_handle == -1
        fclose(node_file_handle);
        reporting.Error('PTKSaveTreeAsNodes:CreateFileFailed', ['I could not create the file ' element_file_name '.']);
    end

    % The first file entry is the number of nodes
    fprintf(node_file_handle, '%u\r\n', num_branches);
    fprintf(element_file_handle, '%u\r\n', num_branches);

    % Create the first node from the trachea start point. The first element
    % connects this point to itself
    first_point = tree_root.StartPoint;

    dicom_coordinates = MimImageCoordinateUtilities.ConvertFromPTKCoordinates([first_point.CoordX, first_point.CoordY, first_point.CoordZ], coordinate_system, template_image);
    start_x_mm = dicom_coordinates(1);
    start_y_mm = dicom_coordinates(2);
    start_z_mm = dicom_coordinates(3);

    start_radius_mm = tree_root.Radius;
    start_density_mgml = tree_root.Density;

    PrintNodeToFile(node_file_handle, 0, start_x_mm, start_y_mm, start_z_mm, start_radius_mm, start_density_mgml, false);
    PrintElementToFile(element_file_handle, 0, 0);

    % The first branch in the tree will be converted to a node at its endpoint.
    % We need this to be connected to the first node we created above, so set
    % its parent index to zero, the index of the first node
    linear_branch_list(1).TemporaryIndex = 0;
    current_node_index = 1;

    for branch = linear_branch_list
        reporting.UpdateProgressStage(current_node_index, num_branches);

        parent_node_index = branch.TemporaryIndex;

        % Write branch to node file (the end point will be the node coordinate).
        PrintBranchToFileAsNode(node_file_handle, current_node_index, branch, coordinate_system, template_image);

        % Write (parent, current_node_index) to element file
        PrintElementToFile(element_file_handle, parent_node_index, current_node_index)

        % Set the parent node index for each child branch to be the current branch
        % node index
        if ~isempty(branch.Children)
            for child = branch.Children
                child.TemporaryIndex = current_node_index;
            end
        end

        current_node_index = current_node_index + 1;
    end

    % Close the files
    fclose(node_file_handle);
    fclose(element_file_handle);

    reporting.CompleteProgress;
end

function PrintElementToFile(fid, parent_index, child_index)
    fprintf(fid, '%u,%u\r\n', parent_index, child_index);
end

function PrintBranchToFileAsNode(fid, node_index, branch, coordinate_system, template_image)

    % We define nodes using the end coordinate of each branch, since this is the
    % bifurcation point
    last_point = branch.EndPoint;

    dicom_coordinates = MimImageCoordinateUtilities.ConvertFromPTKCoordinates([last_point.CoordX, last_point.CoordY, last_point.CoordZ], coordinate_system, template_image);
    x_mm = dicom_coordinates(1);
    y_mm = dicom_coordinates(2);
    z_mm = dicom_coordinates(3);

    radius_mm = branch.Radius;
    density_mgml = branch.Density;

    is_final_node = isempty(branch.Children);
    PrintNodeToFile(fid, node_index, x_mm, y_mm, z_mm, radius_mm, density_mgml, is_final_node);
end

function PrintNodeToFile(fid, node_index, x_mm, y_mm, z_mm, radius_mm, density, is_final_node)
    if is_final_node
        is_final_node = 'y';
    else
        is_final_node = 'n';
    end
    output_string = sprintf('%u,%6.6g,%6.6g,%6.6g,%2.2g,%6.6g,%c\r\n', node_index, x_mm, y_mm, z_mm, radius_mm, density, is_final_node);
    fprintf(fid, regexprep(output_string, ' ', ''));
end


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

choice = [0.01 0.03 0.1 0.3 1 3 10 30]';
minError = Inf;
curC = Inf;
cur_sigma = Inf;

for i = 1:8
	for j = 1:8
		model = svmTrain(X, y, choice(i), @(x1, x2) gaussianKernel(x1, x2, choice(j)));
		predictions = svmPredict(model,Xval);
		error = mean(double(predictions ~= yval));
		if error < minError
			minError = error;
			curC = choice(i);
			cur_sigma = choice(j);
		end
	end
end

C = curC;
sigma = cur_sigma;

steps = [ 0.01 0.03 0.1 0.3 1 3 10 30 ];
minError = Inf;
minC = Inf;
minSigma = Inf;

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

end

author: AvaisP


text:
function data=getMultipleAmps_PerCondition(view,selectedROIs,projPhase)
%function data=getMultipleAmps_PerCondition(view,selectedROIs)
% Returns the same data as you would get from plotMultileAmps_PerCondition.
%
%Reference scan is the current scan
mrGlobals;

refScan = getCurScan(view);
if (ieNotDefined('selectedROIs'))
    % Select ROIs
    nROIs=size(view.ROIs,2);
    roiList=cell(1,nROIs);
    for r=1:nROIs
        roiList{r}=view.ROIs(r).name;
    end

    selectedROIs = find(buttondlg('ROIs to Plot',roiList));
end

nROIs=length(selectedROIs);
if (nROIs==0)
    error('No ROIs selected');
end

if (ieNotDefined('projPhase'))
    projPhase=NaN;
end

nscans = numScans(view);
ROIamps=zeros(nscans,nROIs);
ROIseZ=zeros(nscans,nROIs);
ROImeanPhs=zeros(nscans,nROIs);

for r=1:nROIs

    n=selectedROIs(r);
    view = selectROI(view,n); % is there another way?
    [meanAmps,meanPhs,seZ] = vectorMeans(view);
    ROIamps(:,r)=meanAmps(:);
    ROIseZ(:,r)=seZ(:);
    ROImeanPhs(:,r)=meanPhs(:);
    xstr{r}=int2str(r);
    roiName{r}=view.ROIs(selectedROIs(r)).name;
    fprintf(['\n#%d :',roiName{r}],r);

end

if (~isnan(projPhase))
    ROIprojAmp=ROIamps.*(cos(projPhase-ROImeanPhs));
    data.projAmp=ROIprojAmp;
end

data.amps =ROIamps;
data.seZ= ROIseZ;
data.meanPhs= ROImeanPhs;
data.nROIs=nROIs;
data.ROIname=roiName;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pass = test_fevalm( pref )
% Test diskfun/fevalm

if ( nargin == 0)
    pref = chebfunpref;
end

tol = 100*pref.cheb2Prefs.chebfun2eps;
rng(2016);

% Check empty diskfun:
f = diskfun;
t = pi*(2*rand(5,1) - 1);
r = 2*rand(5,1)-1;
B = fevalm(f, t, r);
pass(1) = isempty( B );

% Check rank 1 diskfun:
f = diskfun(@(t,r) r.*sin(t), 'polar');
t = pi*(2*rand(5,1) - 1);
r = rand(5,1);
[tt, rr] = meshgrid( t, r);
A = feval(f, tt, rr, 'polar');
B = fevalm(f,t,r);
pass(2) = norm( A - B ) < tol;

% Check essentially one dimensional function:
f = diskfun(@(t,r) exp(-r.^2), 'polar');
t = pi*(2*rand(5,1) - 1);
r = rand(5,1);
[tt, rr] = meshgrid( t, r);
A = feval(f, tt, rr, 'polar');
B = fevalm(f, t, r);
pass(3) = norm( A - B ) < tol;

end


text:
function id = mapWindowIndices(map,mapWindow)
%
% id = mapWindowIndices(map,mapWindow)
%
% Returns indices of vector map within mapWindow
%
% if mapWindow(1)<mapWindow(2) returns mapWindow(1) <= map <= mapWindow(2)
% if mapWindow(1)>mapWindow(2) returns map >= mapWindow(1) or map <= mapWindow(2)
%
% djh, 7/98

if diff(mapWindow)>0
  id = find(map>=mapWindow(1) & map<=mapWindow(2));
else
  id = find(map>=mapWindow(1) | map<=mapWindow(2));
end


text:
function x = elastix_transf_imcoord2(x, t)
% elastix_transf_imcoord  Convert 2D pixel coordinates using
% elastix/transformix transform.
%
% Let x be the (x, y)-coordinates of a pixel in an image (instead of
% row/column). This function computes the coordinates that the same pixel
% will have in the new image after transformix or elastix transforms.
%
% XT = elastix_transf_imcoord2(X, T)
%
%   X is a two-column matrix with the Cartesian coordinates of one or more
%   points in the input image.
%
%   T is a struct produced by elastix with the details of the image
%   transformation. See help elastix for details.
%
%   XT has the same size as X, and contains the coordinates of the same
%   points in the output image.
%
% See also: elastix, transformix, elastix_read_file2param,
% elastix_read_reg_output, elastix_write_param2file,
% elastix_compose_afftransf.

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

% check arguments
narginchk(2, 2);
nargoutchk(0, 1);

if (size(x, 2) ~= 2)
    error('X must have two columns (2D points)')
end
if (~isstruct(t))
    error('T must be a struct with the transform parameters')
end
if (isfield(t, 'Direction') && any(t.Direction ~= [1 0 0 1]))
    error('Not implemented for Direction ~= [1 0 0 1]')
end

% defaults
if (~isfield(t, 'CenterOfRotationPoint'))
    t.CenterOfRotationPoint = [0 0];
end

% select type of transform
switch (t.Transform)

    case 'SimilarityTransform' % similarity transform

        if (length(t.TransformParameters) ~= 4)
            error('SimilarityTransform must have 4 parameters')
        end

        % transformation nomenclature
        s =     t.TransformParameters(1);
        theta = t.TransformParameters(2);
        tx =    t.TransformParameters(3);
        ty =    t.TransformParameters(4);
        cx =    t.CenterOfRotationPoint(1);
        cy =    t.CenterOfRotationPoint(2);

        % transform points
        %
        % note that this is the inverse of the similarity transform. The
        % reason is the way that ITK, and by extension elastix/transformix,
        % implement image transformations
        x = [cos(theta) sin(theta); -sin(theta) cos(theta)] ...
            * [x(:, 1)' - tx - cx; x(:, 2)' - ty - cy] / s;
        x = x';
        x(:, 1) = x(:, 1) + cx;
        x(:, 2) = x(:, 2) + cy;

    otherwise

        error('Transform not implemented')

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function s = s_to_rot13 ( s )

%*****************************************************************************80
%
%% S_TO_ROT13 "rotates" the alphabetical characters in a string by 13 positions.
%
%  Discussion:
%
%    Two applications of the routine will return the original string.
%
%  Examples:
%
%    Input:                      Output:
%
%    abcdefghijklmnopqrstuvwxyz  nopqrstuvwxyzabcdefghijklm
%    Cher                        Pure
%    James Thurston Howell       Wnzrf Guhefgba Ubjryy
%    0123456789                  56789012345
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   24 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string to be "rotated".
%
%    Output, string S, the rotated string.
%
  s_length = length ( s );

  for i = 1 : s_length
    s(i) = ch_to_rot13 ( s(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This file is part of TREEQSM.
%
% TREEQSM is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function QSM = simplify_qsm(QSM,MaxOrder,SmallRadii,ReplaceIterations,Plot,Disp)

% ---------------------------------------------------------------------
% SIMPLIFY_QSM.M   Simplifies cylinder QSMs by restricting the maximum
%                       branching order, by removing thin branches, and by
%                       replacing two concecutive cylinders with a longer cylinder
%
% Version 2.0.0
% Latest update     4 May 2022
%
% Copyright (C) 2015-2022 Pasi Raumonen
% ---------------------------------------------------------------------
%
% Inputs:
% QSM           QSM-structure, output of treeqsm.m, must contain only one model
% MaxOrder      Maximum branching order, higher order branches removed
% SmallRadii    Minimum acceptable radius for a branch at its base
% ReplaceIterations Number of iterations for replacing two concecutive
%                     cylinders inside one branch with one longer cylinder
% Plot          If true/1, then plots the cylinder models before and
%                 after the simplification
% Disp          If Disp == 1, then display the simplication results
%                 (the number of cylinders after each step). If
%                 Disp == 2, then display also the treedata results for
%                 the original and simplified QSMs. If Disp == 0, then
%                 nothing is displayed.
%
% Output:
% Modified QSM      NOTICE: cylinder, branch and treedata are modified.

% Changes from version 1.1.0 to 2.0.0, 4 May 2022:
% 1) Added modification of branch and treedata structures based on the
%     modified cylinders
% 2) Added input for plotting and displaying the results
% 3) Corrected some bugs that could cause errors in some special cases

if nargin <= 4
  Plot = 0;
  Disp = 1;
elseif nargin <= 5
  Disp = 1;
end

if Disp == 2
  inputs = QSM.rundata.inputs;
  display_treedata(QSM.treedata,inputs)
end
% Plot the cylinder model before the simplification
if Plot
  plot_cylinder_model(QSM.cylinder,'branch',1,20,1)
end

%% Maximum branching order
c = QSM.cylinder;
nc = size(c.radius,1);
if Disp >= 1
  disp([' ',num2str(nc),' cylinders originally'])
end

% Cylinders with branching order up to MaxBranchOrder
SmallOrder = c.BranchOrder <= MaxOrder;
N = fieldnames(c);
n = max(size(N));
for i = 1:n
  c.(N{i}) = c.(N{i})(SmallOrder,:);
end

% Modify topology information
Ind = (1:1:nc)';
m = nnz(SmallOrder);
Ind(SmallOrder) = (1:1:m)';
I = c.parent > 0;
c.parent(I) = Ind(c.parent(I));
I = c.extension > 0;
c.extension(I) = Ind(c.extension(I));

if Disp == 1
  nc = nnz(SmallOrder);
  disp([' ',num2str(nc),' cylinders after branching order simplification'])
end

%% Small branches
if nargin >= 3 && SmallRadii > 0

  nc = size(c.radius,1);
  % Determine child branches
  BPar = QSM.branch.parent;
  nb = size(BPar,1);
  BChi = cell(nb,1);
  for i = 1:nb
    P = BPar(i);
    if P > 0
      BChi{P} = [BChi{P}; i];
    end
  end

  % Remove branches whose radii is too small compared to its parent
  Large = true(nc,1);
  Pass = true(nb,1);
  for i = 1:nb
    if Pass(i)
      if QSM.branch.diameter(i) < SmallRadii
        B = i;
        BC = BChi{B};
        while ~isempty(BC)
          B = [B; BC];
          BC = vertcat(BChi{BC});
        end
        Pass(B) = false;
        m = length(B);
        for k = 1:m
          Large(c.branch == B(k)) = false;
        end
      end
    end
  end

  % Modify topology information
  Ind = (1:1:nc)';
  m = nnz(Large);
  Ind(Large) = (1:1:m)';
  I = c.parent > 0;
  c.parent(I) = Ind(c.parent(I));
  I = c.extension > 0;
  c.extension(I) = Ind(c.extension(I));

  % Update/reduce cylinders
  for i = 1:n
    c.(N{i}) = c.(N{i})(Large,:);
  end

  if Disp >= 1
    nc = nnz(Large);
    disp([' ',num2str(nc),' cylinders after small branch simplification'])
  end
end

%% Cylinder replacing
if nargin >= 4 && ReplaceIterations > 0

  % Determine child cylinders
  nc = size(c.radius,1);
  CChi = cell(nc,1);
  for i = 1:nc
    P = c.parent(i);
    if P > 0
      PE = c.extension(P);
      if PE ~= i
        CChi{P} = [CChi{P}; i];
      end
    end
  end

  % Replace cylinders
  for j = 1:ReplaceIterations

    nc = size(c.radius,1);
    Ind = (1:1:nc)';
    Keep = false(nc,1);
    i = 1;
    while i <= nc
      t = 1;
      while i+t <= nc && c.branch(i+t) == c.branch(i)
        t = t+1;
      end
      Cyls = (i:1:i+t-1)';
      S = c.start(Cyls,:);
      A = c.axis(Cyls,:);
      L = c.length(Cyls);
      if t == 1 % one cylinder in the branch
        Keep(i) = true;
      elseif ceil(t/2) == floor(t/2) % even number of cylinders in the branch
        I = (1:2:t)'; % select 1., 3., 5., ...
        % Correct radii, axes and lengths
        E = S(end,:)+L(end)*A(end,:);
        S = S(I,:);
        m = length(I);
        if m > 1
          A = [S(2:end,:); E]-S(1:end,:);
        else
          A = E-S(1,:);
        end
        L = sqrt(sum(A.*A,2));
        A = [A(:,1)./L A(:,2)./L A(:,3)./L];
        cyls = Cyls(I);
        Keep(cyls) = true;
        V = pi*c.radius(Cyls).^2.*c.length(Cyls);
        J = (2:2:t)';
        V = V(I)+V(J);
        R = sqrt(V./L/pi);
        c.radius(cyls) = R;

      else % odd number of cylinders
        I = [1 2:2:t]'; % select 1., 2., 4., 6., ...
        % Correct radii, axes and lengths
        E = S(end,:)+L(end)*A(end,:);
        S = S(I,:);
        l = L(1);
        a = A(I,:);
        m = length(I);
        if m > 2
          a(2:end,:) = [S(3:end,:); E]-S(2:end,:);
        else
          a(2,:) = E-S(2,:);
        end
        A = a;
        L = sqrt(sum(A.*A,2));
        L(1) = l;
        A(2:end,:) = [A(2:end,1)./L(2:end) A(2:end,2)./L(2:end) A(2:end,3)./L(2:end)];
        cyls = Cyls(I);
        Keep(cyls) = true;
        V = pi*c.radius(Cyls).^2.*c.length(Cyls);
        J = (3:2:t)';
        V = V(I(2:end))+V(J);
        R = sqrt(V./L(2:end)/pi);
        c.radius(cyls(2:end)) = R;
      end

      if t > 1
        % Modify cylinders
        c.length(cyls) = L;
        c.axis(cyls,:) = A;
        % Correct branching/topology information
        c.PositionInBranch(cyls) = (1:1:m)';
        c.extension(cyls) = [cyls(2:end); 0];
        c.parent(cyls(2:end)) = cyls(1:end-1);
        par = c.parent(cyls(1));
        if par > 0 && ~Keep(par)
          par0 = c.parent(par);
          if Keep(par0) && c.extension(par0) == par
            c.parent(cyls(1)) = par0;
          end
        end

        % Correct child branches
        chi = vertcat(CChi{Cyls});
        if ~isempty(chi)
          par = c.parent(chi);
          J = Keep(par);
          par = par(~J)-1;
          c.parent(chi(~J)) = par;

          par = c.parent(chi);
          rp = c.radius(par);
          sp = c.start(par,:);
          ap = c.axis(par,:);
          lc = c.length(chi);
          sc = c.start(chi,:);
          ac = c.axis(chi,:);
          ec = sc+[lc.*ac(:,1) lc.*ac(:,2) lc.*ac(:,3)];
          m = length(chi);
          for k = 1:m
            [d,V,h,B] = distances_to_line(sc(k,:),ap(k,:),sp(k,:));
            V = V/d;
            sc(k,:) = sp(k,:)+rp(k)*V+B;
          end
          ac = ec-sc;
          [ac,lc] = normalize(ac);
          c.length(chi) = lc;
          c.start(chi,:) = sc;
          c.axis(chi,:) = ac;
        end
      end

      i = i+t;
    end
    % Change topology (parent, extension) indexes
    m = nnz(Keep);
    Ind(Keep) = (1:1:m)';
    I = c.parent > 0;
    c.parent(I) = Ind(c.parent(I));
    I = c.extension > 0;
    c.extension(I) = Ind(c.extension(I));

    % Update/reduce cylinders
    for i = 1:n
      c.(N{i}) = c.(N{i})(Keep,:);
    end

    if j < ReplaceIterations
      % Determine child cylinders
      nc = size(c.radius,1);
      CChi = cell(nc,1);
      for i = 1:nc
        P = c.parent(i);
        if P > 0
          PE = c.extension(P);
          if PE ~= i
            CChi{P} = [CChi{P}; i];
          end
        end
      end
    end
  end

  if Disp >= 1
    nc = size(c.radius,1);
    disp([' ',num2str(nc),' cylinders after cylinder replacements'])
  end
end
if Disp >= 1
  nc = size(c.radius,1);
  disp([' ',num2str(nc),' cylinders after all simplifications'])
end

%% Updata the QSM
% Update the branch
branch = branches(c);

% Update the treedata
inputs = QSM.rundata.inputs;
inputs.plot = 0;
% Display
if Disp == 2
  inputs.disp = 2;
else
  inputs.disp = 0;
end
treedata = update_tree_data(QSM,c,branch,inputs);

% Update the cylinder, branch, and treedata of the QSM
QSM.cylinder = c;
QSM.branch = branch;
QSM.treedata = treedata;

% Plot the cylinder model after the simplification
if Plot
  plot_cylinder_model(QSM.cylinder,'branch',2,20,1)
end

end % End of main function

function display_treedata(treedata,inputs)
%% Generate units for displaying the treedata
Names = fieldnames(treedata);
n = size(Names,1);
Units = zeros(n,3);
m = 23;
for i = 1:n
  if ~inputs.Tria && strcmp(Names{i},'CrownVolumeAlpha')
    m = i;
  elseif inputs.Tria && strcmp(Names{i},'TriaTrunkLength')
    m = i;
  end
  if strcmp(Names{i}(1:3),'DBH')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(end-2:end),'ume')
    Units(i,:) = 'L  ';
  elseif strcmp(Names{i}(end-2:end),'ght')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(end-2:end),'gth')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(1:3),'vol')
    Units(i,:) = 'L  ';
  elseif strcmp(Names{i}(1:3),'len')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(end-2:end),'rea')
    Units(i,:) = 'm^2';
  elseif strcmp(Names{i}(1:3),'loc')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(end-4:end),'aConv')
    Units(i,:) = 'm^2';
  elseif strcmp(Names{i}(end-5:end),'aAlpha')
    Units(i,:) = 'm^2';
  elseif strcmp(Names{i}(end-4:end),'eConv')
    Units(i,:) = 'm^3';
  elseif strcmp(Names{i}(end-5:end),'eAlpha')
    Units(i,:) = 'm^3';
  elseif strcmp(Names{i}(end-2:end),'Ave')
    Units(i,:) = 'm  ';
  elseif strcmp(Names{i}(end-2:end),'Max')
    Units(i,:) = 'm  ';
  end
end
%% Display treedata
disp('------------')
disp('  Tree attributes before simplification:')
for i = 1:m
  v = change_precision(treedata.(Names{i}));
  if strcmp(Names{i},'DBHtri')
    disp('  -----')
    disp('  Tree attributes from triangulation:')
  end
  disp(['  ',Names{i},' = ',num2str(v),' ',Units(i,:)])
end
disp('  -----')
end

author: InverseTampere

repo: TreeQSM


text:
%BZreaction animation
%Belousov-Zhabotinsky Reaction animation
%This MATLAB code is converted from Processing code available in this link
%http://www.aac.bartlett.ucl.ac.uk/processing/samples/bzr.pdf

%version 2. Corrected the drift of pixels as suggested
%           by Jonh.

xres=75; %x resolution
yres=75; %y resolution

a=rand(xres,yres,2);
b=rand(xres,yres,2);
c=rand(xres,yres,2);
c_a = zeros(xres,yres);
c_b = zeros(xres,yres);
c_c = zeros(xres,yres);
p = 1;
q = 2;
img=zeros(xres,yres,3);

for k=1:100
         c_a = 0*c_a;
         c_b = 0*c_b;
         c_c = 0*c_c;
         for m=1:xres
            for n=1:yres
               for mm = m:m+2
                  for nn = n:n+2
                    c_a(m,n) =c_a(m,n)+ a( mod(mm+xres,xres)+1 ,mod(nn+yres,yres)+1 ,p);
                    c_b(m,n) =c_b(m,n)+ b( mod(mm+xres,xres)+1 ,mod(nn+yres,yres)+1 ,p);
                    c_c(m,n) =c_c(m,n)+ c( mod(mm+xres,xres)+1 ,mod(nn+yres,yres)+1 ,p);
                  end
               end
            end
         end

        %correction of pixel drift
        c_a = circshift(c_a,[2 2]);
        c_b = circshift(c_b,[2 2]);
        c_c = circshift(c_c,[2 2]);

        c_a =c_a/ 9.0;
        c_b =c_b/ 9.0;
        c_c =c_c/ 9.0;

        a(:,:,q) = double(uint8(255*(c_a + c_a .* (c_b - c_c))))/255;
        b(:,:,q) = double(uint8(255*(c_b + c_b .* (c_c - c_a))))/255;
        c(:,:,q) = double(uint8(255*(c_c + c_c .* (c_a - c_b))))/255;

        img(:,:,1)=c(:,:,q);
        img(:,:,2)=b(:,:,q);
        img(:,:,3)=a(:,:,q);

        if  p == 1
          p = 2; q = 1;
        else
          p = 1; q = 2;
        end

        image(uint8(255*hsv2rgb(img)))
        axis equal off
        drawnow
end

author: Sable

repo: mcbench-benchmarks


text:
function [kPa] = ftH2O2kPa(ftH2O)
% Convert pressure from feet of water column at 4 degrees to kilopascals
% Chad Greene 2012
kPa = ftH2O*2.98907;


text:
function g = getGroupSuperNonOverlap(I, param)

    N = size(I,1) * size(I, 2);

    g = sparse(zeros(N,1));
    g = diag(g);

    imlab = vl_xyz2lab(vl_rgb2xyz(I)) ;
    imlab = single(imlab);
    slicParam = param.superpixel.slicParam;

    groupCount = 0;

    segments = vl_slic(imlab, slicParam(1,1),  slicParam(1,2)) ;
    segments = segments(:);

    for j = 0:segments(N)
        g(:, groupCount+j+1) = segments==j;
    end
    groupCount  = groupCount+ segments(N)+1;

    % croping G
    g =g(:, 1:groupCount);

end

author: andrewssobral

repo: lrslibrary


text:
function [ n_data, p, q, n, d ] = dedekind_sum_values ( n_data )

%*****************************************************************************80
%
%% DEDEKIND_SUM_VALUES returns some values of the Dedekind sum.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Hans Rademacher, Emil Grosswald,
%    Dedekind Sums,
%    Mathematics Association of America, 1972,
%    LC: QA241.R2.
%
%  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, integer P, Q, the arguments of the function.
%
%    Output, integer N, D, the numerator and denominator of the function value.
%
  n_max = 95;

  p_vec = [ ...
     1,  1,  1,  1,  1,  1,  1,  1,  1,  1, ...
     1,  1,  1,  1,  1,  1,  1,  1,  1,  1, ...
     2,  2,  2,  2,  2,  2,  2,  2,  2,  2, ...
     3,  3,  3,  3,  3,  3,  3,  3,  3,  3, ...
     3,  3,  3,  3,  4,  4,  4,  4,  4,  4, ...
     4,  4,  4,  4,  5,  5,  5,  5,  5,  5, ...
     5,  5,  5,  5,  5,  5,  5,  5,  5,  5, ...
     6,  6,  6,  6,  6,  6,  6,  7,  7,  7, ...
     7,  7,  7,  7,  7,  7,  7,  7,  7,  7, ...
     7,  7,  7,  7,  7 ...
  ];

  q_vec = [ ...
     1,  2,  3,  4,  5,  6,  7,  8,  9, 10, ...
    11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...
     1,  3,  5,  7,  9, 11, 13, 15, 17, 19, ...
     1,  2,  4,  5,  7,  8, 10, 11, 13, 14, ...
    16, 17, 19, 20,  1,  3,  5,  7,  9, 11, ...
    13, 15, 17, 19,  1,  2,  3,  4,  6,  7, ...
     8,  9, 11, 12, 13, 14, 16, 17, 18, 19, ...
     1,  5,  7, 11, 13, 17, 19,  1,  2,  3, ...
     4,  5,  6,  8,  9, 10, 11, 12, 13, 15, ...
    16, 17, 18, 19, 20 ...
  ];

  n_vec = [ ...
     0,  0,  1,  1,  1,  5,  5,  7, 14,  3, ...
    15, 55, 11, 13, 91, 35, 20, 34, 51, 57, ...
     0, -1,  0,  1,  4,  5,  4,  7,  8, 21, ...
     0,  0, -1,  0, -1,  1,  0,  3,  1,  3, ...
     5,  5,  9,  3,  0,  1, -1,  1, -4,  3, ...
    -1, 19,  0, 11,  0,  0, -1,  1, -5, -1, ...
    -1,  4, -5, -1,  0,  3, -5,  1,  2, 11, ...
     0,  1, -5,  5, -4,  5, -9,  0,  0,  1, ...
    -1,  0,  5, -7, -4,  0, -3,  1,  4, -7, ...
    -3,  1, -2,  3,  3 ...
  ];

  d_vec = [ ...
     1,  1, 18,  8,  5, 18, 14, 16, 27,  5, ...
    22, 72, 13, 14, 90, 32, 17, 27, 38, 40, ...
     1, 18,  1, 14, 27, 22, 13, 18, 17, 38, ...
     1,  1,  8,  1, 14, 16,  1, 22, 13, 14, ...
    32, 17, 38,  8,  1, 18,  5, 14, 27, 22, ...
    13, 90,  1, 38,  1,  1, 18,  8, 18, 14, ...
    16, 27, 22, 72,  1, 14, 32, 17, 27, 38, ...
     1,  5, 14, 22, 13, 17, 38,  1,  1, 18, ...
     8,  1, 18, 16, 27,  1, 22, 72, 13, 18, ...
    32, 17, 27, 38,  8 ...
  ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    p = 0;
    q = 0;
    n = 0;
    d = 0;
  else
    p = p_vec(n_data);
    q = q_vec(n_data);
    n = n_vec(n_data);
    d = d_vec(n_data);
  end

  return
end


text:
function planC = CERRImport_newplan(optName, archiveDir, saveFile)
%function planC = CERRImport_newplan(optName,archiveDir,saveDir)
%
%27 Oct 2006    KU   Modified CERRImport.m for import of plans to an existing study.
%
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
global stateS

fname1 = '';

[fname1, pathname] = uigetfile('*.*','Select 0000 file to begin importing.');
if fname1==0;
    disp('RTOG scan aborted.');
    return
else
    archiveDir = pathname;
    if exist('optName')==0
        optName='CERROptions.m';
    end
   optS = opts4Exe(optName);
end

%Start diary and timer for import log.
startTime = now;
tmpFileName = tempname;
diary(tmpFileName);

%%%%Call the import tool: IMPORT occurs here.
    temp_planC = importRTOGDir_newplan(optS, archiveDir, fname1);
%%%%

%Check that there are plans in temp_planC.
if isempty(temp_planC{temp_planC{end}.dose}) ||...
        (isempty(temp_planC{temp_planC{end}.dose}(1).doseArray) && isempty(temp_planC{temp_planC{end}.dose}(1,end).doseArray))
    disp('No plans to import.');
    sentence1 = 'There are no plans to be imported.';
    Zmsgbox=msgbox(sentence1, 'modal');
    waitfor(Zmsgbox);
    clear temp_planC;
    return
end

nplans = 0;
for j = 1:length(temp_planC{temp_planC{end}.dose})
    nplans = nplans + 1;
    if ischar(temp_planC{temp_planC{end}.dose}(1,j).fractionGroupID)
        ID = temp_planC{temp_planC{end}.dose}(1,j).fractionGroupID;
    else
        ID = num2str(temp_planC{temp_planC{end}.dose}(1,j).fractionGroupID);
    end
    newPlanName = ID;

    %Check if plan is already in current study.
    match = 0;
    if length(planC{planC{end}.dose}) >=1
        for k=1:length(planC{planC{end}.dose})
            if isequal(temp_planC{temp_planC{end}.dose}(1,j).doseArray, planC{planC{end}.dose}(1,k).doseArray)
                match = 1;
                nplans = nplans - 1;
                disp(['Plan "', ID, '" already exists in the current study.']);
                disp(['Skipping import of plan "',  ID, '".']);
                sentence1 = (['Plan "', ID, '" is already in the current study. ',...
                    'It will not be imported.']);
                Zmsgbox=msgbox(sentence1, 'modal');
                waitfor(Zmsgbox);
            end
        end
    end

    if match ~= 1   %plan is not yet in current study
        % Check dose for duplicate fractionGroupID names.
        if length(planC{planC{end}.dose}) >=1
            n = 0;
            match = 1;
            while match == 1
                match = 0;
                for k=1:length(planC{planC{end}.dose})
                    if ischar(planC{planC{end}.dose}(1,k).fractionGroupID)
                        idString = planC{planC{end}.dose}(1,k).fractionGroupID;
                    else
                        idString = num2str(planC{planC{end}.dose}(1,k).fractionGroupID);
                    end

                    if strcmpi(idString, newPlanName)
                        match = 1;
                        n=n+1;
                        newPlanName = [ID,' (',num2str(n),')'];
                    end
                end
            end
            temp_planC{temp_planC{end}.dose}(1,j).fractionGroupID = newPlanName;
            if n >= 1
                disp(['A plan with the name "', ID, '" already exists in the current study.']);
                disp(['Changing new plan name to "',  newPlanName, '".']);
                sentence1 = (['A plan with the name "', ID, '" already exists in the current study. ',...
                            'Changing new plan name to "',  newPlanName, '".']);
                Zmsgbox=msgbox(sentence1, 'modal');
                beep
                waitfor(Zmsgbox);
            end

            % Check beam geometry for duplicate fractionGroupID names.
            if ~strcmpi(ID, newPlanName)
                for i=1:length(temp_planC{temp_planC{end}.beamGeometry})
                    if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                        tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                    else
                        tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                    end
                    if strcmpi(tempID, ID)
                        temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID = newPlanName;
                    end
                end
            end

            % Check DVH's for duplicate fractionIDOfOrigin names.
            if ~strcmpi(ID, newPlanName)
                for i=1:length(temp_planC{temp_planC{end}.DVH})
                    if ischar(temp_planC{temp_planC{end}.DVH}(1,i).fractionIDOfOrigin)
                        tempID = temp_planC{temp_planC{end}.DVH}(1,i).fractionIDOfOrigin;
                    else
                        tempID = num2str(temp_planC{temp_planC{end}.DVH}(1,i).fractionIDOfOrigin);
                    end
                    if strcmpi(tempID, ID)
                        temp_planC{temp_planC{end}.DVH}(1,i).fractionIDOfOrigin = newPlanName;
                    end
                end
            end
        end

        % Add associated scan and dose UID's to dose UID fields.
        scanUID = planC{planC{end}.scan}(1).scanUID;
        temp_planC{temp_planC{end}.dose}(1,j).assocScanUID = scanUID;
        temp_planC{temp_planC{end}.dose}(1,j).doseUID = createUID('dose');

        % Add new plan to planC.
        if isempty(planC{planC{end}.dose}) || isempty(planC{planC{end}.dose}(1,end).doseArray)
            planC{planC{end}.dose} = temp_planC{temp_planC{end}.dose}(1,j);
            ok = 0;
            for i=1:length(temp_planC{temp_planC{end}.beamGeometry})     %add beams
                if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                    tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                else
                    tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                end
                if strcmpi(tempID, newPlanName)
                    ok = ok+1;
                    if ok == 1
                        planC{planC{end}.beamGeometry} = temp_planC{temp_planC{end}.beamGeometry}(1,i);
                    else
                        planC{planC{end}.beamGeometry} = [planC{planC{end}.beamGeometry}, temp_planC{temp_planC{end}.beamGeometry}(1,i)];
                    end
                end
            end
        else
            if isfield(planC{planC{end}.dose}(1),'cachedMask')
                temp_planC{temp_planC{end}.dose}(1,j).cachedMask = [];
            end
            planC{planC{end}.dose} = dissimilarInsert(planC{planC{end}.dose}, temp_planC{temp_planC{end}.dose}(1,j));

            for i=1:length(temp_planC{temp_planC{end}.beamGeometry})
                if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                    tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                else
                    tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                end
                if strcmpi(tempID, newPlanName)
                    planC{planC{end}.beamGeometry} = dissimilarInsert(planC{planC{end}.beamGeometry}, temp_planC{temp_planC{end}.beamGeometry}(1,i));
                end
            end
        end

    end
end

if nplans == 0
    disp('No plans to import.');
    sentence1 = 'No plans were found to import!';
    Zmsgbox=msgbox(sentence1, 'modal');
    beep
    waitfor(Zmsgbox);
    return
elseif length(planC{planC{end}.dose})>=1 && ~isempty(planC{planC{end}.dose}(1).doseArray)
    planC{planC{end}.dose}(1,end).doseNumber = length(planC{planC{end}.dose});

    %switch to new dose after a short pause.
    pause(.1);
    sliceCallBack('selectDose', num2str(planC{planC{end}.dose}(1,end).doseNumber));

    %switch to scan #1 if not already selected.
    sliceCallBack('selectScan', num2str(1));

    %Turn on dose and switch to colorwash if not already on.
    if stateS.doseToggle == -1
        stateS.doseToggle = - stateS.doseToggle;
        set(findobj('tag', 'doseToggle'), 'checked', 'on');
    end
    stateS.optS.dosePlotType = 'colorwash';
    stateS.doseChanged = 1;
    stateS.doseSetChanged = 1;
    stateS.doseDisplayChanged = 1;
    stateS.structsChanged = 1;
    sliceCallBack('refresh')

    sentence1 = horzcat('The current study has been updated with the new plan(s). ',...
        'The changes have not been saved.  Please remember to save changes ',...
        'before closing the study.');
    Zmsgbox=msgbox(sentence1, 'modal');
    beep
    waitfor(Zmsgbox);
end

%Stop diary and write it to planC.
diary off;
endTime = now;
logC = file2cell(tmpFileName);
delete(tmpFileName);
indexS = planC{end};

temp_planC{temp_planC{end}.importLog} = [];
temp_planC{temp_planC{end}.importLog}.importLog = logC;
temp_planC{temp_planC{end}.importLog}.startTime = datestr(startTime);
temp_planC{temp_planC{end}.importLog}.endTime = datestr(endTime);

planC{indexS.importLog} = [planC{indexS.importLog}, temp_planC{temp_planC{end}.importLog}];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR