Bojibranch

.gitignore

@@ -7,3 +7,5 @@ docs/_*
 old/
 *.mat
 *.ipynb_checkpoints
+*.nii
+*.nii.gz

AirQuantAddPath.m

@@ -32,7 +32,7 @@
 end
 
 if nargout > 0
-    varargout{1} = AirQuantDir
+    varargout{1} = AirQuantDir;
 end
 % configure expected paths for AirQuant
 results_dir = fullfile(AirQuantDir, 'results');
@@ -42,4 +42,18 @@
     mkdir(results_dir)
 end
 
+function a = AQstruct2array(s)
+% AQSTRUCT2ARRAY Convert structure with doubles to an array. 
+% Incase struct2array is not available.
+% 
+% From
+% https://uk.mathworks.com/matlabcentral/answers/1717910-was-struct2array-removed-from-matlab
+%
+
+% Convert structure to cell
+c = struct2cell(s);
+% Construct an array
+a = [c{:}];
+
+end
 end

library/AirQuant.m

@@ -40,7 +40,7 @@ function toNii(obj, filename, options)
             if isprop(obj,'sourceinfo')
                 info = obj.sourceinfo;
             elseif isprop(obj,'network')
-                info = obj.network.sourceinfo;
+                info = obj.network.header;
             end
 
             if exist('info','var')

library/dataio/CropVol.m

@@ -38,15 +38,15 @@
         % incase lower lims are < index 1, set to min possible
         lims = varargin{2};
         for ii = 1:3
-            if lims(ii) < 1
-                lims(ii) = 1;
+            if lims(ii,1) < 1
+                lims(ii,1) = 1;
             end
         end
         % incase upper lims are > size dim, set to max possible
         j = 1;
-        for ii = 4:6
-            if lims(ii) > s(j)
-                lims(ii) = s(j);
+        for ii = 1:3
+            if lims(ii,2) > s(j)
+                lims(ii,2) = s(j);
             end
             j = j+1;
         end

library/dataio/ReorientVolume.m

@@ -18,7 +18,14 @@
     %
 
     % get affine matrix
-    aff_raw_RAS = meta.Transform.T;
+    aff_raw_RAS_loaded = meta.Transform.T;
+    % set values that are basically 0 to 0
+    zero_tol = aff_raw_RAS_loaded < 1e-6 & aff_raw_RAS_loaded > -1e-6;
+    aff_raw_RAS = aff_raw_RAS_loaded;
+    aff_raw_RAS(zero_tol) = 0;
+    if any(aff_raw_RAS_loaded ~= aff_raw_RAS)
+        warning('Oblique orientation, precision loss within 1e-6')
+    end
     % remove origin information
     aff_raw_RAS(4,1:3) = [0,0,0];
     % remove spacing information

library/graph/skel_2_digraph.m

@@ -1,4 +1,4 @@
-function [digraphout, glinkout, gnode] = skel_2_digraph(skel, method)
+function [digraphout, glinkout, gnode, isloops] = skel_2_digraph(skel, method)
     % Generate digraph from skeleton.
     %
     % Generate digraph from skeleton using skel2graph library.
@@ -36,16 +36,34 @@
     if nargin < 2
         method = 'topnode';
     end
+
+    isloops = 0;
 
     [gadj,gnode,glink] = Skel2Graph3D(skel,1);
     % choose originating node using chosen method
 
 
     % Create digraph with edges in both directions, loop through
     % branches and remove opposing direction to originating node.
-    G = digraph(gadj);
-    bins = conncomp(G);
 
+    edges_og = [[glink.n1]', [glink.n2]'];
+    edges_rev = [[glink.n2]', [glink.n1]'];
+    weights = cellfun(@numel, {glink.point});
+    weights = [weights'; weights'];
+    edges_twoway = [edges_og; edges_rev];
+    Edgetable = table(edges_twoway,weights,'VariableNames',{'EndNodes','Weight'});
+    G = digraph(Edgetable);
+    nodelist = 1:numnodes(G);
+
+    G_loop_check = digraph(gadj);
+    if length(glink) ~= height(G_loop_check.Edges)/2
+        warning('Skeleton appears to contain multiple edges between the same nodes. There may be unexpected behaviour')
+        isloops = 1;
+    end
+
+    [bins,binsize] = conncomp(G);
+    originnode = zeros(max(bins),1);
+
     if isnumeric(method)
         assert(length(unique(bins)) < 2, ['Can only be used with one ' ...
             'connected component.'])
@@ -56,21 +74,77 @@
 
     elseif strcmp(method,'topnode')
             % use most superiour node as origin for each subgraph
-            originnode = zeros(max(bins),1);
             for ii = 1:max(bins)
                 binbool = (bins==ii);
-                nodelist = 1:numnodes(G);
                 binidx = nodelist(binbool);
                 gnodeii = gnode(binidx);
                 [~, binorigin] = max([gnodeii.comz]);
                 originnode(ii) = binidx(binorigin);
             end
+
+    elseif strcmp(method,'carina')
+        % airway specific method - largest cc first
+        [~, I] = max(binsize);
+        binbool = (bins==I);
+        binidx = nodelist(binbool);
+
+        g = subgraph(G, binidx);
+        gnodeI=gnode(binidx);
+        %%% identfy carina by most central position
+        pos = [[gnodeI.comx]', [gnodeI.comy]', [gnodeI.comz]'];
+        center = size(skel)/2;
+        centralist_node=zeros(length(pos),1);
+        for i=1:length(pos)
+            centralist_node(i) = norm(center-pos(i,:));
+        end
+        % rank centralist_node - minimise
+        [~, ~, centralist_node_rank] = unique(centralist_node);
+
+        %%% identify carina by greatest node sum weight
+        weightiest_node = centrality(g,'outdegree',...
+        'importance',g.Edges.Weight);
+
+        % rank weightiest - maximise
+        [~,~,weightiest_node_rank_rev] = unique(weightiest_node);
+        weightiest_node_rank = max(weightiest_node_rank_rev) - ...
+            weightiest_node_rank_rev+1;
+
+        %%% identify carina by graph closeness (stucturally central)
+        closest_node = centrality(g,'outcloseness');
+        [~,~,closest_node_rank_rev] = unique(closest_node);
+        closest_node_rank = max(closest_node_rank_rev) - ...
+            closest_node_rank_rev+1;
+
+        %%% determine carina by rank of all 3 measures
+        all_ranks = centralist_node_rank+weightiest_node_rank+closest_node_rank;
+        [~, carina_node] = min(all_ranks);
+
+        %%% identify origin by most superior neighbor
+        %%% determine trachea by superior pos node to carina
+        % pos of carina neighbors
+        carina_neighbors = predecessors(g, carina_node);
+        neighbor_pos = pos(carina_neighbors,:);
+        [~, neighbor_pos_idx] = max(neighbor_pos(:,3));
+        binorigin=carina_neighbors(neighbor_pos_idx);
+        originnode(I) = binidx(binorigin);
+
+        % remaining uses topnode
+        for ii = 1:max(bins)
+            if ii == I
+                continue
+            end
+            binbool = (bins==ii);
+            binidx = nodelist(binbool);
+            gnodeii = gnode(binidx);
+            [~, binorigin] = max([gnodeii.comz]);
+            originnode(ii) = binidx(binorigin);
+        end
     else
         error('Invalid method')
     end
 
 
-    % BF search from carina node to distal for each disconnected graph.
+    % BF search from origin node to distal for each disconnected graph.
     allnode_discovery = cell(max(bins),1);
     for ii = 1:max(bins)
         allnode_discovery{ii} = bfsearch(G,originnode(ii));
@@ -130,7 +204,7 @@
     end
     labels = [1:length(glink)]';
     Edgetable = table(edges,weights,labels,'VariableNames',{'EndNodes', 'Weight', 'Label'});
-
+    
     digraphout = digraph(Edgetable);
     % add node properties from Gnode
     digraphout.Nodes.comx(:) = [gnode(:).comx];
@@ -139,6 +213,7 @@
     digraphout.Nodes.ep(:) = [gnode(:).ep];
     digraphout.Nodes.label(:) = [1:length(gnode)]';
 
+    try
     % BFS per graph in digraph
     % first convert digraph to graph to get CCs
     bins = conncomp(digraphout,'Type','weak','OutputForm','cell');
@@ -158,12 +233,17 @@
         E = [E; sub_E];
     end
 
-    assert(length(E)==length(glink),'Number of edges in new edge BFS not expected.')
+    assert(length(E)==length(glink))
 
     % convert E indicies to glink indicies
     E_glink = digraphout.Edges.Label(E);
 
     % reorder glink
     glinkout = glink(E_glink);
-
+
+    catch
+        glinkout = glink;
+        warning('Failed to reorder edges by BFS. Likely due to presence loops.')
+    end
+
 end

library/perpinterp/PlaneInterpVol.m

@@ -1,4 +1,4 @@
-function out_plane = PlaneInterpVol(vol, voxdim, point, normal, options)
+function out_plane = PlaneInterpVol(whole_vol, voxdim, whole_point, normal, options)
     % short desc
     %
     % long desc
@@ -28,9 +28,9 @@
 
     % Construct vol grid
     arguments
-        vol 
+        whole_vol 
         voxdim {mustBeNumeric}
-        point (3,1) {mustBeNumeric}
+        whole_point (3,1) {mustBeNumeric}
         normal (3,1) {mustBeNumeric}
         options.plane_sz (1,1) = 40
         options.sample_sz (1,1) = 0.5
@@ -39,6 +39,8 @@
         options.offgrid_val (1,1) {mustBeNumeric} = 0
     end
 
+    [vol,point] = min_tube_vol(whole_vol,whole_point,voxdim,options.plane_sz);
+
     image_sz = single(size(vol));
     [x_domain , y_domain , z_domain] = ...
         meshgrid(1:image_sz(2),1:image_sz(1),1:image_sz(3));
@@ -51,9 +53,12 @@
     plane_grid = Grids_coords_for_plane(basis_vecs(:,3),...
         basis_vecs(:,2), point, options.plane_sz, options.sample_sz);
 
+
     if options.gpu
-        assert(strcmp(options.method,'linear'),'Can only use "linear" method with gpu optimisation.')
-        assert(license('test','distrib_computing_toolbox'),' "Parallel Computing Toolbox required for gpu optimisation.')
+        assert(strcmp(options.method,'linear'),'Can only use "linear" method with gpu optimisation.')        
+        verstruct = struct2cell(ver);
+        assert(contains([verstruct{:}], 'Parallel Computing Toolbox'), "Parallel Computing Toolbox required for gpu optimisation. Otherwise set gpu=0")
+
 
         % use gpu
         x_domain = gpuArray(x_domain);

library/perpinterp/min_tube_vol.m

@@ -0,0 +1,39 @@
+function [min_vol,min_point] = min_tube_vol(vol,point,zooms,plane_size)
+    % Returns the minimum volume required t make perpendicular slices of
+    % tubes.
+
+    arguments
+        vol (:,:,:)
+        point (3,1)
+        zooms (3,1)
+        plane_size (1,1)
+    end
+
+    % identify indices of bounding box
+    min_BB_mm = point - plane_size/2;
+    max_BB_mm = point + plane_size/2;
+
+    % transform to index
+    min_BB = floor(min_BB_mm./zooms);
+    max_BB = ceil(max_BB_mm./zooms);
+
+    % check BB does not exceed vol extremes
+    for ii = 1:length(min_BB)
+        if min_BB(ii) < 1
+            min_BB(ii) = 1;
+        end
+
+        maxdim = size(vol,ii);
+        if max_BB(ii) > maxdim
+            max_BB(ii) = maxdim;
+        end
+    end
+
+    % transform vol
+    min_vol = vol(min_BB(1):max_BB(1),min_BB(2):max_BB(2),...
+        min_BB(3):max_BB(3));
+
+    % transform point, plus zooms due to index 1
+    min_point = point - min_BB.*zooms+zooms;
+end
+

library/util/default_dock_fig.m

@@ -0,0 +1,13 @@
+function default_dock_fig(flag)
+% set default behaviour for new figures to be docked.
+if nargin < 1
+    flag = 1;
+end
+
+if flag == 1
+    set(0,'DefaultFigureWindowStyle','docked')
+else
+    set(0,'DefaultFigureWindowStyle','normal')
+end
+
+end

library/util/fig_save.m

@@ -0,0 +1,16 @@
+function fig_save(fig_handle,path,ext)
+% saves figure as both .mat and .png (determined by ext).
+%
+% path must not have an ext.
+%
+
+if nargin < 3
+    ext = '.png';
+end
+
+% save as MATLAB proprietary .fig
+savefig(fig_handle, path)
+% save as image file
+exportgraphics(fig_handle, strcat(path,ext));
+end
+