diff --git a/.gitignore b/.gitignore index 47d9cf9..f1fc345 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,5 @@ docs/_* old/ *.mat *.ipynb_checkpoints +*.nii +*.nii.gz diff --git a/AirQuantAddPath.m b/AirQuantAddPath.m index 7e56837..61528f2 100644 --- a/AirQuantAddPath.m +++ b/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 diff --git a/data/airquant/.empty b/data/airquant/.empty new file mode 100644 index 0000000..e69de29 diff --git a/data/example/github_demo_raw.nii.gz b/data/example/github_demo_raw.nii.gz deleted file mode 100644 index 3a2e0b6..0000000 Binary files a/data/example/github_demo_raw.nii.gz and /dev/null differ diff --git a/data/example/github_demo_seg.nii.gz b/data/example/github_demo_seg.nii.gz deleted file mode 100644 index 351bb44..0000000 Binary files a/data/example/github_demo_seg.nii.gz and /dev/null differ diff --git a/data/example/github_demo_seg_PTKskel.nii.gz b/data/example/github_demo_seg_PTKskel.nii.gz deleted file mode 100644 index 5ddc2c8..0000000 Binary files a/data/example/github_demo_seg_PTKskel.nii.gz and /dev/null differ diff --git a/library/AirQuant.m b/library/AirQuant.m index 236354b..26667dd 100644 --- a/library/AirQuant.m +++ b/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') diff --git a/library/dataio/CropVol.m b/library/dataio/CropVol.m index 418dcf5..04fd702 100644 --- a/library/dataio/CropVol.m +++ b/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 diff --git a/library/dataio/ReorientVolume.m b/library/dataio/ReorientVolume.m index 5b63591..5812acf 100644 --- a/library/dataio/ReorientVolume.m +++ b/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 diff --git a/library/graph/skel_2_digraph.m b/library/graph/skel_2_digraph.m index d03cdd1..04b167b 100644 --- a/library/graph/skel_2_digraph.m +++ b/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,6 +36,8 @@ if nargin < 2 method = 'topnode'; end + + isloops = 0; [gadj,gnode,glink] = Skel2Graph3D(skel,1); % choose originating node using chosen method @@ -43,9 +45,25 @@ % 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 diff --git a/library/perpinterp/PlaneInterpVol.m b/library/perpinterp/PlaneInterpVol.m index d4a360f..326be65 100644 --- a/library/perpinterp/PlaneInterpVol.m +++ b/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); diff --git a/library/perpinterp/min_tube_vol.m b/library/perpinterp/min_tube_vol.m new file mode 100644 index 0000000..b8e1461 --- /dev/null +++ b/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 + diff --git a/library/util/AQstruct2array.m b/library/util/AQstruct2array.m deleted file mode 100644 index 600a357..0000000 --- a/library/util/AQstruct2array.m +++ /dev/null @@ -1,11 +0,0 @@ -function a = AQstruct2array(s) -% AQSTRUCT2ARRAY Convert structure with doubles to an array. -% Incase struct2array is not available. -% -% -% -% https://uk.mathworks.com/matlabcentral/answers/1717910-was-struct2array-removed-from-matlab -% Convert structure to cell -c = struct2cell(s); -% Construct an array -a = [c{:}]; \ No newline at end of file diff --git a/library/util/default_dock_fig.m b/library/util/default_dock_fig.m new file mode 100644 index 0000000..de64138 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/library/util/fig_save.m b/library/util/fig_save.m new file mode 100644 index 0000000..c80d2fb --- /dev/null +++ b/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 + diff --git a/library/util/grid_preview.m b/library/util/grid_preview.m new file mode 100644 index 0000000..9a5fe01 --- /dev/null +++ b/library/util/grid_preview.m @@ -0,0 +1,50 @@ +function grid = grid_preview(tar_file, n_rows, n_cols, savepath) + % Given a tar file of patches, display a grid of them and save as png. + + if nargin < 4 + savepath = 0; + else + savepath = parse_filename_extension(savepath, '.png'); + end + + % Load the tar file + atemp_dir = tempname; + + % check tar path + untar(tar_file,atemp_dir) + + % Get the number of patches + patches = dir(fullfile(atemp_dir, '*.tif')); + npatches = length(patches); + randi = randperm(npatches); + + % check enough patches for grid + if npatches < n_rows*n_cols + warning('Not enough patches in tar file for grid. Resizing grid.') + n_cols = floor(npatches/n_rows); + end + + % get the size of the patches + patch_size = size(imread(fullfile(atemp_dir,patches(1).name))); + + % synthesize the grid + grid = zeros(patch_size(1)*n_rows, patch_size(2)*n_cols); + k=1; + for i=1:n_rows + for j=1:n_cols + grid((i-1)*patch_size(1)+1:i*patch_size(1), (j-1)*patch_size(2)+1:j*patch_size(2)) = imread(fullfile(atemp_dir,patches(randi(k)).name)); + k=k+1; + end + end + + % delete the temp dir + rmdir(atemp_dir, 's'); + + % balance the grid + grid = imadjust(rescale(grid)); + + % save the grid as png + if savepath + imwrite(single(grid), savepath); + end +end \ No newline at end of file diff --git a/mandocs/README.md b/mandocs/README.md deleted file mode 100644 index 8df383d..0000000 --- a/mandocs/README.md +++ /dev/null @@ -1,26 +0,0 @@ -# AirQuant Docs -By Ashkan Pakzad https://ashkanpakzad.github.io - -See [/readme.md](/readme.md) before reading this document regarding install, about and purpose. - - -This software primarily revolves around the [library/AirQuant.m](library/AirQuant.m) class. Every CT case will exist as a new AirQuant class object. When a case is loaded in, it is first initialised, performing a number of checks and short processes. -Currently there is only one method for taking airway measurements, the [FWHMesl method](https://doi.org/10.1117/12.595283). This first relies on the CT scan to be interpolated at right angles to the principle airway axis. For scientific reading on the concept, please see [K Quan et al.](https://doi.org/10.1117/12.2292306) on which this software is based. - -As most of the code behind this software is in [library/AirQuant.m](library/AirQuant.m) it is recommended that you explore it in MATLAB by collapsing all code and expanding as necessary. Set code collapsing for each section MATLAB preferences. - -Example use can also be found in [scripts/](../scripts/) and example data is in [data/](../data/). - - - - -# Contents: - * [Data, Initialisation & Batching](/docs/basic.md) - prepping data and initialisation an AirQuant object. - * [PTKskel: Skeletonisation](/docs/skel.md) - Skeletonisation and the packaged algorithm. - * [CT Airway Interpolation](/docs/interp.md) - Interpolation of the CT along the airway principal axis. - * [FWHMesl method](/docs/fwhm.md) - Measuring the airway. - * [Tapering Metrics](/docs/taper.md) - Tapering metrics that can be generated. - * [Visualisation](/docs/vis.md) - Methods for visualising global data. - * [Segmental Visualisation](/docs/segvis.md) - Methods for visualising data for individual segments. - - * [Future Ambitions](/docs/future.md) - watch this space. diff --git a/mandocs/basic.md b/mandocs/basic.md deleted file mode 100644 index ee1571c..0000000 --- a/mandocs/basic.md +++ /dev/null @@ -1,111 +0,0 @@ -# Basic use - -## Preparing your data -Your data will need to be in Nifti format (.nii or .nii.gz) and read in using `niftiread` and `niftiinfo`. The airway segmentation should be one completely [26-connected](https://en.wikipedia.org/wiki/Pixel_connectivity) object, AQ takes the largest connected object and discards everything else in the segmentation. - -It is recommended that the packaged skeletonisation algorithm is used `PTKskel`. The resultant skeleton of this algorithm is suited for use with AirQuant. - -## Initialisation -Setting up the AQ object for a case is straightforward. It requires loading up the CT, metadata, airway segmentation and skeleton from [NIFTI format](https://brainder.org/2012/09/23/the-nifti-file-format/). The metadata variable is the output from `niftiinfo`. - -Inputs: -* Original CT -* niftiinfo() header information -* Airway Segmentation -* Airway Centreline/Skeleton* -* filename to save results - -Output: -* AirQuant object class - -*Notes* - -Initialising the AQ class calls for the skeleton to be parsed into a directional-graph such that edges face from proximal to the distal direction. [see etc.] This forms the backbone of the AQ class that allows each airway segment to be individually analysed and processed. In addition nodes are bifurcation (or even trifurcation points). - -Lobes are classified by examining the graph nodes and their relative position, e.g. the node at the end of the right major bronchi will be more right than the node at the end of the left major bronchi by definition. This is largely based on the method described by [Gu et al.](doi.org/10.1155/2012/382806) with the addition of classifying the Left lingular separately. This can be generated graphically, see [Visualisation](/docs/vis.md). - -Each segment's generation is identified by counting the number of nodes from the carina node to origin node of an airway segment. The carina node is first identified by a centrality metric of the first pass to convert the skeleton into a digraph. - -All data is converted into LPS orientation (standard for DICOM) once AQ is initialised. So data inside the AQ object may not align with your raw data outside of AQ. It is necessary to declare an orientation within AQ, most CT data is stored in this orientation already. See [Understanding 3D medical image orientation for programmers](https://medium.com/@ashkanpakzad/understanding-3d-medical-image-orientation-for-programmers-fcf79c7beed0). - -Functions that take longer to run often call the `save` method, saving the current object's state to a ".mat" file as declared by the `savename` property. - -For more information see [library/AirQuant.m](library/AirQuant.m) > methods > %%INITIALISATION > AirQuant. - -*Example* -``` -% names input files, give fullpaths if not in matlab path -% add AirQuant library to path -AirQuantDir = AirQuantAddPath(); -casename = 'github_demo'; -dataset = 'example' -results_dir = fullfile(AirQuantDir,'results',dataset, casename); - -% Get filenames -CT_name = [casename, '_raw.nii.gz']; -seg_name = [casename, '_seg.nii.gz']; -skel_name = [casename, '_seg_PTKskel.nii.gz']; - -% Load CT data as double -meta = niftiinfo(CT_name); -CT = double(niftiread(meta)); - -% Load Airway segmentation and its skeleton as logicals -S = logical(niftiread(seg_name)); -skel = logical(niftiread(skel_name)); - -% Initialise AirQuant -% Parses CT, segmentation and skeleton to compute airway tree graph -% savename is given to automatically save/load results. -savename = fullfile(results_dir, [casename, '_AQ.mat']); -AQ = AirQuant(CT, meta, S, skel, savename); -``` -## Reloading -Previously saved AirQuant objects can be reloaded by calling the AirQuant class with the savename/path only. - -*Example* -``` -% add AirQuant library to path -AirQuantDir = AirQuantAddPath(); -dataset = 'example' -casename = 'github_demo'; -results_dir = fullfile(AirQuantDir,'results', dataset, casename); - -% Load AirQuant object -savename = fullfile(results_dir, [char(casename, '_AQ.mat']); -AQ = AirQuant(savename); -``` - -## Batch Processing -AirQuant comes packed with an extended batch function/script under [library/runAQ.m](library/runAQ.m). -This function is designed to run AirQuant on several cases by providing a configuration structure. -It will also run certain visualisation methods and save them to .fig and .png formats. -For alternative batch scripting, it is recommended that you use this function as a template. -An example configuration script can be found under [scripts/example_config.m](scripts/example_config.m). - -*Notes* - -* This requires all data to be stored exactly within `AirQuant/data/datasetname`. -* The command line output of each case is saved in its own log file under `AirQuant/data/datasetname/casename/casename_log.txt`. -* If any of the CT, segmentation or skeletonisation datafiles do not exist, it will inform and skip but not throw an error. This is to prevent time wasted in extended jobs. -* If the AQ object already exists in results, it will not skip but instead reload and process it. The AQ traversing method regularly saves the AQ object. If the job is interrupted it just picks up from where it left off. - -*Example* -``` -AirQuantDir = AirQuantAddPath(); - -%%% set up the config structure -config = []; -% must be string array, using double dash quotes. -config.casenames = ["github_demo"]; % required -config.dataset = 'example'; - -% suffix of each file (these are defaults if not provided) -config.CTsuf = '_raw.nii.gz'; -config.segsuf= '_seg.nii.gz'; -config.skelsuf = '_seg_PTKskel.nii.gz'; - -%%% pass to AirQuant runner -runAQ(config); - -``` diff --git a/mandocs/future.md b/mandocs/future.md deleted file mode 100644 index 5d4a87e..0000000 --- a/mandocs/future.md +++ /dev/null @@ -1 +0,0 @@ -wait for it.... diff --git a/mandocs/fwhm.md b/mandocs/fwhm.md deleted file mode 100644 index 06ce7cf..0000000 --- a/mandocs/fwhm.md +++ /dev/null @@ -1,61 +0,0 @@ -# FWHMesl method - -This section will cover the methods of AirQuant that measures the airways in the interpolated airway slices by the [Full Width at Half Maximum Edge-cued Segmentation Limited method](https://doi.org/10.1117/12.595283). This is executed for the inner lumen wall edge, wall peak attenuation point and outer wall edge. It is strongly encouraged to read [K Quan et al.](https://doi.org/10.1117/12.2292306) on which this software is based. The two higher level functions are documented first, the lower level methods are mentioned to their purpose in a discussion of the theory below. - -### FindFWHMall(obj) -Ultimately this whole process can be called to be executed on all airways by calling this function. - -This step should only take a few minutes. You may receive some warning signs and errors in its process, these can be largely ignored unless you see it occurring for pretty much every single airway. - -All results are stored in the AQ object itself. See the lower level functions that are called by this method below for more information on each step. - -*Notes* - -This step effectively calls `FindAirwayBoundariesFWHM` to run on every airway branch except the trachea. If the user insists on running this function of the trachea, they can call this function explicitly with the trachea index. - -The interpolation parameters are automatically derived from the CT Voxel Size and stored in the AQ object properties. The user can override these parameters by setting the properties themselves after initialisation and before running this method. `num_rays` and `ray_interval`. See the Theory section for more information. - -*Example* -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -FindFWHMall(AQ) -``` - -### FindAirwayBoundariesFWHM(obj, link_index) -It is rare to request the interpolation for just one airway, but this function is available if so, it is intended as a lower level method that can be called by `AirwayImageAll`. It may be useful for testing an individual case and calling some segmental visualisations. - -Unlike the interpolation methods, the AQ object is not saved at the end of processing of each one as this could be significantly slow. It is recommended that the user call the `save` function themselves if they are only processing measuremnts and not computing metrics in the same session. - -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -FindAirwayBoundariesFWHM(AQ, 41) -``` - - -# Theory -There are 3 listed steps to this process. The first two involve numerous loops of code but are very fast. - -* Raycasting from centre outwards. -* Identifying stop points. -* Ellipse fitting. - -The FWHMesl method in essence treats the 1D intesity profile across the wall as a bell shaped peak from which the inner edge is idenitified by the left [FWHM](https://en.wikipedia.org/wiki/Full_width_at_half_maximum) and the outer edge by the right FWHM. However, this profile can be noisy and is therefore made more robust by using information from the interpolated segmentation where the wall is expected to be near the edge of the segmentation. Note that the interpolated segmentation is no longer binary but is now continuous between [0, 1]. This figure illustrates how the method works to supplement the text. - -![FWHMesl method figure](./fwhmesl.png) -*(a) Interpolated CT showing showing raycast lines from centre where the 1D profile is captured. (b) the equivalent interpolated segmentation. (c) the 1D intensity of the CT showing the bell shape. (d) the 1D intensity of the segmentation along the same cast ray in the CT. L_B and W_B mark the lumen and wall boundary respectively. P_A marks the wall peak, also known as the peak attenuation point.* - - -## Raycasting and stop points. -The centre of the airway is identified by `Return_centre_pt_image` and confirmed by the segmentation too `Check_centre_with_segmentation`. - -The rays are then cast by `Raycast` from this point outwards to get then 1D intensity profiles outwards from this point. The radial ray sampling is set by `obj.ray_interval` (0.2 degrees by default). - -The FWHM method is executed on each individual pair of ray 1D intensity profiles by `AirQuant.computeFWHM`. This identifies the stop points for each of the inner wall, wall peak-attenuation and outer wall. The FWHMesl method identifies the peak closest to the segmentation image edge. The peak is then travelled downwards along the negative gradient until it plateaus on both sides; it computes the mid point between these trough points and the peak to get the width at half maximum. Hence the name. - -Finally for each set of raycast stop points an ellipse is fit by `ComputeEllipses` in 2D. The results are saved in the private property `obj.FWHMesl`. diff --git a/mandocs/fwhmesl.png b/mandocs/fwhmesl.png deleted file mode 100644 index 25da486..0000000 Binary files a/mandocs/fwhmesl.png and /dev/null differ diff --git a/mandocs/github_demo_airway3.png b/mandocs/github_demo_airway3.png deleted file mode 100644 index 367a64f..0000000 Binary files a/mandocs/github_demo_airway3.png and /dev/null differ diff --git a/mandocs/github_demo_avgd.png b/mandocs/github_demo_avgd.png deleted file mode 100644 index 87aa156..0000000 Binary files a/mandocs/github_demo_avgd.png and /dev/null differ diff --git a/mandocs/github_demo_gen3d.png b/mandocs/github_demo_gen3d.png deleted file mode 100644 index 3cb7614..0000000 Binary files a/mandocs/github_demo_gen3d.png and /dev/null differ diff --git a/mandocs/github_demo_graph.png b/mandocs/github_demo_graph.png deleted file mode 100644 index 0aa4b3b..0000000 Binary files a/mandocs/github_demo_graph.png and /dev/null differ diff --git a/mandocs/github_demo_lobe3d.png b/mandocs/github_demo_lobe3d.png deleted file mode 100644 index 7e5d0b6..0000000 Binary files a/mandocs/github_demo_lobe3d.png and /dev/null differ diff --git a/mandocs/github_demo_ortho.png b/mandocs/github_demo_ortho.png deleted file mode 100644 index 2f1d866..0000000 Binary files a/mandocs/github_demo_ortho.png and /dev/null differ diff --git a/mandocs/github_demo_skel3d.png b/mandocs/github_demo_skel3d.png deleted file mode 100644 index 2d3f867..0000000 Binary files a/mandocs/github_demo_skel3d.png and /dev/null differ diff --git a/mandocs/github_demo_splinetree.png b/mandocs/github_demo_splinetree.png deleted file mode 100644 index 7d5c7aa..0000000 Binary files a/mandocs/github_demo_splinetree.png and /dev/null differ diff --git a/mandocs/github_demo_splinevec.png b/mandocs/github_demo_splinevec.png deleted file mode 100644 index 6697718..0000000 Binary files a/mandocs/github_demo_splinevec.png and /dev/null differ diff --git a/mandocs/github_demo_tree.png b/mandocs/github_demo_tree.png deleted file mode 100644 index 7d691ca..0000000 Binary files a/mandocs/github_demo_tree.png and /dev/null differ diff --git a/mandocs/github_demo_view3dgen.png b/mandocs/github_demo_view3dgen.png deleted file mode 100644 index ac52c85..0000000 Binary files a/mandocs/github_demo_view3dgen.png and /dev/null differ diff --git a/mandocs/github_demo_view3dlobe.png b/mandocs/github_demo_view3dlobe.png deleted file mode 100644 index 7bd140a..0000000 Binary files a/mandocs/github_demo_view3dlobe.png and /dev/null differ diff --git a/mandocs/interp.md b/mandocs/interp.md deleted file mode 100644 index 2294c19..0000000 --- a/mandocs/interp.md +++ /dev/null @@ -1,68 +0,0 @@ -# CT Airway interpolation - -The key step to being able to measure the airways at perpendicular angles lay in interpolating the CT image and even the segmentation. This section will cover the methods of AirQuant that completes this step and give an idea of the concept behind them. Strongly encouraged to read [K Quan et al.](https://doi.org/10.1117/12.2292306) on which this software is based. The two higher level functions are documented first, the lower level methods are mentioned to their purpose in a discussion of the theory below. - -### AirwayImageAll(obj) -Ultimately this whole process can be called to be executed on all airways by calling this function. The time it takes can vary from case to case. Smaller voxel size, larger volume, larger airways, longer airways (larger subjects), more airways segmented can all cause this step to take longer. - -It is recommended to execute this step on a high performance workstation or better yet a dedicated cluster. Ortherwise you may find that this process will take a long time (>24h) and easily run out of memory. On average, a case will take 6h (i7 processor, 32gb ram etc.). - -All results are stored in the AQ object itself. See the lower level functions that are called by this method below for more information but each step. - -*Notes* - -This step effectively calls `CreateAirwayImage` to run on every airway branch except the trachea. If the user insists on running this function of the trachea, they can call this function explicitly with the trachea index. - -The interpolation parameters are automatically derived from the CT Voxel Size and stored in the AQ object properties. The user can override these parameters by setting the properties themselves after initialisation and before running this method. `max_plane_sz` `plane_sample_sz`, `spline_sample_sz`, `plane_scaling_sz`, `min_tube_sz`. See the Theory section for more information. - -*Example* -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -AirwayImageAll(AQ) -``` - -### CreateAirwayImage(obj, link_index) -It is rare to request the interpolation for just one airway, but this function is available if so, it is intended as a lower level method that can be called by `AirwayImageAll`. It may be useful for testing an individual case and calling some segmental visualisations. - -Due to the fact that it can take a long time to process each airway, the AQ object is saved at the end of processing of each one by calling the `save` method. - -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -CreateAirwayImage(AQ, 41) -``` - - -# Theory -There are 5 listed steps to this process. Those concerning the spline are fast and those regarding the final interpolation step take much longer. - -* Compute the spline of an individual airway. -* Identify interpolation sample points along the spline -* Identify the normal vector to each sample, i.e. tangent along the spline. -* Compute the plane of each sampled point. -* Interpolate the CT (and seg) to the computed plane (The most computationally expensive step). - -## Spline computation -A spline is fitted to the skelton points of an individual airway in order to interpolate the centreline continuously at subvoxel points. Tangental vectors along the spline at these sampled spline points will be at right angles to the airway. - -`ComputeSpline` The skeleton points of the current airway and the parent (more central) airway are smoothed with a moving average filter, a spline is then fit to the current airway processed skeleton points only. The splines are stored in `obj.Splines` - -`ComputeSplinePoints` will then identify the sample points at a set interval `obj.spline_sample_sz` (half the smallest voxel dimension by default). The spline parametrized sample points are also stored in `obj.Splines` at second column, however the corresponding arclength (length along the spline) points are stored in `obj.arclength`. - -`AirQuant.ComputeNormal` MATLAB spline functions identify the real image point and the first derivative along the spline computes the normal. - -## Plane computation -`InterpolateCT` does the bulk of this step. The orthonormal vectors of the tangent are computed by the function `Orthonormal_basis_with_tangent_vector` and the plane grid at the sample point is generated by the function `Grids_coords_for_plane`. - -The maximum plane size is set by `obj.max_plane_sz` (40mm x 40mm by default) and the in-plane sample rate is set by `obj.plane_sample_sz` (half the smallest voxel dimension by default). - -AQ has the capability of adapting the actual plane size of a given point based on an indication of size by the airway segmentation. `ComputeDmapD` gets the value of the distance transform for a given CT point and the size of the plane is set by `obj.plane_scaling_sz` (5 by default) times this value. Note that the distance transform is only computed once and stored at `obj.Dmap`. - -The interpolation of the segmentation and CT is carried out by the MALTAB `interp3` function to populate the generated plane. - -The interpolated CTs and segmentation images are stored in `obj.TraversedImage` and `obj.TraversedSeg` respectively . diff --git a/mandocs/segvis.md b/mandocs/segvis.md deleted file mode 100644 index 4d3fac1..0000000 --- a/mandocs/segvis.md +++ /dev/null @@ -1,65 +0,0 @@ -# Segmental Visualisation - -These visualisation methods are intended for individual airway branches, hence the input argument `link_index`. Specify which branch to investigate by setting this value, it can be useful to see which airway corresponds to which by looking at the digraph first with `plot(AQ)`. - -### h = PlotSplineVecs(obj, subsamp, link_index) -Plot the normal vectors along the spline points of a particular airway branch. -The current branch will appear in blue and neighboring airways in yellow for context. - -`link_index` can be an individual airway or multiple provided in a single array. It will not be as detailed with a single airway if multiple are provided. If left blank, it will plot the splines of all airways which may take a while. - -One can provide a `subsamp` value to indicate how often to sample the spline points. e.g. greater the number, the fewer sampled vectors will be shown. - -This function requires the user to have directly or indirectly called the method `ComputeSpline` for all airways see [CT Airway Interpolation](/docs/interp.md). - -*Notes* - -It may become easier to interpret the resultant plot by increasing the `subsamp` value. - -Caution: This method has not been robustly tested. - -*Example* -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -figure; -h = PlotSplineVecs(AQ, 2, 41) -``` - -![Example result of graph plot](./github_demo_splinevec.png) - -### PlotAirway3(obj, link_index) -Interactively scroll through the interpolated airway slices of a particular branch. The computed inner and outer airway wall edge are overlaid on top with further information shown at the top in the yellow bar. In slices where measurements failed, no overlay is shown. - -This function requires airway measurements of `link_index` to first be processed, see [CT Airway Interpolation](/docs/interp.md) and [FWHMesl method](/docs/fwhm.md). - -*Example* -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -PlotAirway3(AQ, 41) -``` - -![Example result of graph plot](./github_demo_airway3.png) - -### s = OrthoViewAirway(obj, link_index) -Interactively view the interpolated airway slices of a particular branch. - -Uses MATLAB's built in orthosliceViewer. `s` the output is the slice viewer object see `help orthosliceViewer` for more information about it. - -This function requires airway measurements of `link_index` to first be processed, see [CT Airway Interpolation](/docs/interp.md) and [FWHMesl method](/docs/fwhm.md). - -*Example* -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -s = OrthoViewAirway(AQ, 41) -``` - -![Example result of graph plot](./github_demo_ortho.png) diff --git a/mandocs/skel.md b/mandocs/skel.md deleted file mode 100644 index c1c0987..0000000 --- a/mandocs/skel.md +++ /dev/null @@ -1,28 +0,0 @@ -# Skeletonisation - -The airway skeleton/centreline should be true to the natural centreline of the airway and contain no loops such that it can be broken down into a tree object. AQ parses the skeleton into a network digraph allowing analysis on individual airway segments. - -AQ will warn of any anomalies in the skeleton on initialisation and flag them on visualisation of the skeleton graph but make no effort to solve them. - -A robust skeleton is necessary for successful analysis with AirQuant, to this end a suitable algorithm, [library/PTKskel] is packaged with AirQuant. Example use can also be found in [scripts/example_skel] - -## PTKskel -Provide the filename of the segmentation to be skeletonised. The result is saved in the same folder and the same name with an appended "_PTKskel". The saved nifti file will be in alignment with the original segmentation file, allowing it to be used easily outside of PTK/AirQuant. - -The PTK library requires that the input segmentation be in the MATLAB current folder. - -*Notes* - -This algorithm utilises the library of the [PulmonaryToolKit (PTK) by Tom Doel](https://github.com/tomdoel/pulmonarytoolkit) to skeletonise an already complete segmentation. Original PTK plugins only allows skeletonisation of airway segmentations complete by its own algorithm, this algorithm essentially calls PTK region growing algorithm to re-segment the airways using a propagating wavefront method originating from the trachea, this algorithm parses the airways at the same time. The resultant PTKairways object can then be passed to the PTKskeletonisation library to employ the algorithm based on [Palágyi et al.](doi.org/10.1016/j.compbiomed.2005.05.004). It also checks for loops and removes 'offending' skeleton branches. - -AirQuant has methods to plot the segment and skeleton in one figure, see [Visualisation](/docs/vis.md). - -*Example* -``` -% must be in matlab current path -segname = 'github_demo_seg.nii.gz'; -% Ensure all AirQuant files are in matlab path -AirQuantAddPath(); - -PTKskel(segname); -``` diff --git a/mandocs/taper.md b/mandocs/taper.md deleted file mode 100644 index 08b2afc..0000000 --- a/mandocs/taper.md +++ /dev/null @@ -1,54 +0,0 @@ -# Tapering Metrics - -One of the ultimate goals of this intense data processing is to measure the gradient of tapering of airways. There are two 'schools of thought' built into AirQuant so far. One we will call **long tapering**, this is where you can look at the gradient of tapering from the carina to the distal point of the outermost airways, the other we call **segmental tapering**, where we consider each airway graph edge a segment as an individual unit to measure. - -## Long Tapering - -With this analysis we are looking at the tapering gradient from the carina to the most distal point of every airway end terminal. So we will have the same number of measurements as terminal airway end nodes (i.e nodes with one edge minus the top trachea node). - -### AllTaperResults = ComputeTaperAll(obj) - -We have one high level function which will compute the tapering gradient of every long airway path for all three sets of airway measurements and store it in a table. This is done relatively fast and can be saved in a new variable which the user can export to csv, xlsx etc. if they wish. Furthermore, the table is also stored within the AQ property structure `obj.Specs`. - -*Notes* - -This method calls on `ListTerminalNodes` to get a list of terminal nodes and process each of them through `ConstructTaperPath`. - -This is the measure that K Quan et al. considers and analysed extensively [Ref 1](https://doi.org/10.1117/12.2292306), [Ref 2](https://doi.org/10.1117/1.jmi.6.3.034003) and [Ref 3](http://arxiv.org/abs/1906.12225). Though their method is not exactly the same to process airway measurements. - -*Example* - -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -AllTaperResults = ComputeTaperAll(AQ) -``` - - -## Segmental Tapering - -This analysis is looking at the tapering of each individual airway branch as a single unit. Thus we will have the same number of measurements as there are branches minus the trachea branch. Specifically, there are two measures considered here **Intratapering** and **Intertapering**. The former considers percentage change in airway diameter along the segment relative to the first airway diameter measurement, the latter considers percentage change in average airway diameter relative to the previous segment. - -### SegmentTaperResults = SegmentTaperAll(obj, prunelength) - -We have one high level function which will compute the tapering gradient of every branch for all three sets of airway measurements and store it in a table. This is done relatively fast and can be saved in a new variable which the user can export to csv, xlsx etc. if they wish. Furthermore, the table is also stored within the AQ property structure `obj.Specs`. - -The second argument is optional and sets the length in mm to prune either end of the airway branches, it is a two element array. This is sometimes decided upon in order to avoid natural diameter changes due to bifurcations. For no pruning set the `prunelength` to `[0 0]`. - -*Notes* - -This method calls on `ComputeIntraTaper`, `ComputeIntraTaperAll` and `ComputeInterTaper`. - -This is a measure that [Kuo et al.](doi.org/10.1007/s00330-019-06606-w) considers extensively. - -*Example* - -``` -% reloading processed AQ object. -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% call function -SegmentTaperResults = SegmentTaperAll(AQ, [0 0]]) -``` diff --git a/mandocs/vis.md b/mandocs/vis.md deleted file mode 100644 index 81403bb..0000000 --- a/mandocs/vis.md +++ /dev/null @@ -1,139 +0,0 @@ -# Visualisation -AirQuant has a number of methods for visualising useful data. Often there are 3D visualisations where the airway segmentation is classified and 2D visualisations where the airway graph can be represented. - -## 3D -### PlotMap3D(obj, mode) -Visualise the airway segmentation in 3D and classify by either 'lobe' or 'generation' using MATLAB's standard `isosurface` library. - -*Notes* - -In cases where the result may appear buggy, it is recommended that the user use `View3D`. - -If the entire segmentation is not skeletonised properly i.e. the volume is not fully connected then those detached parts will most likely appear anomalous as they will be labelled by whichever lobe's airways they are physically nearest to. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% classify by generation -figure; -PlotMap3D(AQ, 'generation'); -% classify by lobes -figure; -PlotMap3D(AQ, 'lobe'); -``` - -![Example result of generation plot](./github_demo_gen3d.png) -![Example result of lobe plot](./github_demo_lobe3d.png) - - -### View3D(obj, mode) -An alternative to PlotMap3D, it uses MATLAB's VolumeViewer to visualise and label the airways by generation and lobe. Accepted modes either 'lobe' or 'generation'. - -*Notes* - -If possible, it is recommended that the user use `PlotMap3D` as it result in a better standard. - -If the entire segmentation is not skeletonised properly i.e. the volume is not fully connected then those detached parts will most likely appear anomalous as they will be labelled by whichever lobe's airways they are physically nearest to. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% classify by generation -PlotMap3D(AQ, 'generation'); -% classify by lobes -figure; -PlotMap3D(AQ, 'lobe'); -``` - -![Example result of view3d for generations](./github_demo_view3dgen.png) -![Example result of view3d for lobes](./github_demo_view3dlobe.png) - - -### PlotSegSkel(obj) -plots the airway skeleton within a translucent faced segmentation. This can be useful to identify any problems with the skeleton if results appear suspicious. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% plot airways and skeleton -figure; -PlotSegSkel(AQ); -``` - -![Example result of segskel](./github_demo_skel3d.png) - -## Graph - -### plot(obj, type) -Plot the digraph as a tree structure and specify edge label type as either 'index', 'lobes', 'generation' or 'none'. By default if the second param is not specified it will set the graph edges as the airway index number set within AirQuant. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% show digraph -figure; -plot(AQ); -``` - -![Example result of graph plot](./github_demo_graph.png) - -### PlotTree(obj) -Plot the skeleton with all graph information. Airways are numbered by their index. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% show tree plot -figure; -PlotTree(AQ); -``` - -![Example result of plot tree](./github_demo_tree.png) - -### PlotSplineTree(obj) -Plot all the splines that are fit to the airway skeleton. The colour of each spline is random, but useful for distinguishing. - -This function requires the user to have directly or indirectly called the method `ComputeSpline` for all airways see [CT Airway Interpolation](/docs/interp.md). - -*Notes* - -This resultant plot may not align with other 3D objects due to the way in which spline data is stored in the AQ object at a fundamental level. This will be addressed in future releases. - -*Example* -``` -% after initialising object -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% show tree plot -figure; -PlotSplineTree(AQ); -``` - -![Example result of graph diameter plot](./github_demo_splinetree.png) - -### GraphPlotDiameter(obj) -Plot the airway digraph where edges line thickness correspond proportionally to the average diameter of each airway. - -This function requires all airway measurements to first be processed, see [CT Airway Interpolation](/docs/interp.md) and [FWHMesl method](/docs/fwhm.md). - -*Example* -``` -% after initialising object and measuring all airways -savename = 'results/github_demo/github_demo_AQ.m' -AQ = AirQuant(savename); -% show average diameter plot -figure; -GraphPlotDiameter(AQ); -``` - -![Example result of graph diameter plot](./github_demo_avgd.png) diff --git a/measure/AQEllipse.m b/measure/AQEllipse.m index 4becaad..9e46bd0 100644 --- a/measure/AQEllipse.m +++ b/measure/AQEllipse.m @@ -41,9 +41,16 @@ end % compute values from ellipses - obj.Area(); - obj.NominalDiameter(); - obj.HydraulicDiameter(); + if obj.Rx < 0 || obj.Ry < 0 + % if ellipse RA or RB turns out to be negative + obj.area = nan; + obj.hydraulic_diameter = nan; + obj.diameter = nan; + else + obj.Area(); + obj.NominalDiameter(); + obj.HydraulicDiameter(); + end end function obj = FitEllipse(obj) diff --git a/measure/AirwayCNN.m b/measure/AirwayCNN.m deleted file mode 100644 index b6eb639..0000000 --- a/measure/AirwayCNN.m +++ /dev/null @@ -1,118 +0,0 @@ -classdef AirwayCNN < SuperMeasure - methods - function obj = AirwayCNN(varargin) - % See superclass constructor method SuperMeasure. - % - % long desc - % - % .. todo: add documentation to this function - % - % Args: - % x(type): - % - % Return: - % y(type): - % - obj@SuperMeasure(varargin{:}); - end - - function Measure(obj, modulepath, model_path, subtype) - % Measure airways using a deep learning method - % - % long desc - % - % .. todo: add documentation to this function - % - % Args: - % modulepath(char): absolute path to awyGAN parent directory. - % model_path(char): absolute path to model load checkpoint. - % - % - arguments - obj - modulepath char - model_path char - subtype char {mustBeMember(subtype,{'awyGAN','CNR'})} = 'CNR' - end - - % subtype - switch subtype - case 'awyGAN' - subtype = 'runawyGANmodel'; - case 'CNR' - subtype = 'runCNNmodel'; - end - - - % save patches to temp dir - obj.pixsize = 0.5; - temp_dir = tempname; - obj.tube.ExportOrthoPatches(temp_dir, 'temp') - - % compute image stats - M = cell2mat(obj.tube.source); - meanval = mean(M(:)); - varval = std(M(:)); - - % add module path to python paths and import - GANCNNpath = fullfile(modulepath,'GANCNN'); - pyrun(['import sys; sys.path.append("', modulepath,'");']); - pyrun(['import sys; sys.path.append("', GANCNNpath,'");']); - mod = py.importlib.import_module('AQ_CNR'); - py.importlib.reload(mod); - - pycode = ... - ['import sys; sys.path.append("', modulepath,'");', ... - ' from AQ_CNR import ',subtype, ';',... - ' modelpath = "', model_path, '";', ... - ' data_path = "', temp_dir, '";', ... - ' batch_size = ', num2str(256), ';', ... - ' mean = ', num2str(meanval), ';', ... - ' var = ', num2str(varval), ';', ... - ' out = ',subtype,'(modelpath, data_path, batch_size, mean, var)']; - - pymeasures = pyrun(pycode, 'out'); - measures = double(pymeasures); - - % set up memory - allslices = 1:length(obj.tube.source); - chosenslices = obj.tube.PruneMeasure(allslices); - - total_slices_sz = length(allslices); - ellipse_inner = cell(total_slices_sz, 1); - ellipse_outer = cell(total_slices_sz, 1); - - initslice = chosenslices(1); - - for k = 1:size(measures,1) - imcc = size(obj.tube.source{k}) / 2; - center = imcc + obj.pixsize + measures(k,3:4) / obj.pixsize; - - % process inner ellipses - innerstruct = struct; - innerstruct.center = center; - innerstruct.Rx = measures(k,1)/obj.pixsize(1); - innerstruct.Ry = measures(k,2)/obj.pixsize(1); - innerstruct.rotation = measures(k,7); - - % process outer ellipses - outerstruct = struct; - outerstruct.center = center; - outerstruct.Rx = measures(k,5)/obj.pixsize(1); - outerstruct.Ry = measures(k,6)/obj.pixsize(1); - outerstruct.rotation = measures(k,7); - - % compute ellipses - ellipse_inner{initslice,1} = AQEllipse(obj.pixsize, innerstruct); - ellipse_outer{initslice,1} = AQEllipse(obj.pixsize, outerstruct); - initslice = initslice + 1; - end - obj.measures = cell(2, length(ellipse_inner)); - obj.measures(1,:) = ellipse_inner; - obj.measures(2,:) = ellipse_outer; - - end - - - end -end diff --git a/network/ClinicalAirways.m b/network/ClinicalAirways.m index 1afa8bd..434e39c 100644 --- a/network/ClinicalAirways.m +++ b/network/ClinicalAirways.m @@ -29,8 +29,12 @@ obj.regioncategories.lobe = {'B','RUL','RML','RLL','LUL','LML','LLL','T'}; obj.IdentifyCarinaAndTrachea(); - -% obj.ClassifyLungLobes() + + % attempt to classify into lunglobes + try + obj.ClassifyLungLobes() + catch + end end @@ -121,12 +125,18 @@ % identify upper lobe and lingular MLULLML2 = MLlung(MLULLML).children; [~, ~, MLl2z] = obj.I2S(ClinicalAirways.SkelEnds(MLULLML2)); - [~,MLUL] = max(MLl2z); + % set LML by lowest z descendant [~,MLML] = min(MLl2z); - MLULLML2(MLUL).SetRegionDescendants('lobe','LUL'); MLULLML2(MLML).SetRegionDescendants('lobe','LML'); + + % set remaining descendants to LUL + MLULLML2(MLML) = []; + for idx = 1:length(MLULLML2) +% [~,MLUL] = max(MLl2z); + MLULLML2(idx).SetRegionDescendants('lobe','LUL'); + end - % % Identify right upper lobe + % identify right upper lobe MRlung = MB(rightMBI).children; [~, ~, MRz] = obj.I2S(ClinicalAirways.SkelEnds(MRlung)); [~,MRULi] = max(MRz); @@ -157,12 +167,14 @@ % get ancestor list and find intersection RML_endpath = [RML_eptube.Ancestors().ID]; RLL_endpath = [RLL_eptube.Ancestors().ID]; - [~, intersections] = intersect(RML_endpath, RLL_endpath); - RML_RLLNid = RML_endpath(min(intersections)); + [intersected_B, intersections] = intersect(RML_endpath, RLL_endpath); - % assign right bronchus intermedius - obj.tubes(RML_RLLNid).SetRegion('lobe','B'); - obj.tubes(RML_RLLNid).SetRegion('name','RightIntermedius'); + % assign right bronchus intermedius to branches between carina + % and RML except the major right bronchus + for id = intersected_B(2:end) + obj.tubes(id).SetRegion('lobe','B'); + obj.tubes(id).SetRegion('name','RightIntermedius'); + end % assign RML RML_id = RML_endpath(min(intersections)-1); @@ -184,9 +196,188 @@ end % set gen by lobe -% obj.RunAllTubes('SetRegionGeneration', 'lobe') + obj.RunAllTubes('SetRegionGeneration', 'lobe') + end + + function X = GraphLobarLayout(obj, h, g) + % Order lobes in graph plot to standard layout. + % + % :class:`network.TubeNetwork.Plot` usually orders the lobes randomly. + % this function will manipulate x positions of lobes so that + % they appear in the standard layout. This helps make plots + % more homogenous and therefore comparable to each other. + % + % .. warning: + % This may cause edges to cross. + % + % .. todo: + % Integrate a cross over correction mechanism. + % + % + % Args: + % h = graphical object of graph plot + % g = digraph object + % + % Example: + % >>> run CA_base.m; + % >>> figure; + % >>> [h, g] = AQnet.Plot(colour='lobe'); + % >>> AQnet.GraphLobarLayout(h, g) + % + + X = h.XData; + + % identify root nodes of lobes + lobe_ids = GetTubeValues(obj, 'lobe'); + + % find root nodes of lungs + lobes = {'RUL', 'RML', 'RLL', 'LLL', 'LML', 'LUL'}; + lobe_origins_id = zeros(length(lobes),1); + % tube index origins in same order as lobes + for ii = 1:length(lobes) + lobe = lobes{ii}; + % first occurance of that lobe idx will be the root of it. + lobe_origins_id(ii) = find(strcmp(lobe_ids,lobe),1); + end + + % check if lungs in the correct order. + % this is done by comparing the parent nodes of the RUL and + % LLL. + + % get x pos of the two nodes + RUL_edge = find(g.Edges.ID == lobe_origins_id(1)); + R_node = g.Edges.EndNodes(RUL_edge,1); + LLL_edge = find(g.Edges.ID == lobe_origins_id(4)); + L_node = g.Edges.EndNodes(LLL_edge,1); + % right must be less than left + if X(R_node) > X(L_node) + % swap the right and left lungs at these nodes. + swapnodes(R_node, L_node) + end + + % focus on ordering right lobes + + all_lobe_node_xpos = get_all_lobe_node_xpos(); + % if RUL not in 1st pos, swap RUL and right BI + % find the two decendents of the RUL node + [~,I] = min(all_lobe_node_xpos); + if I ~= 1 + % identify right BI + RUL_node = lobe_origins_id(1); + R_node_suc = successors(g,R_node); + RBI_node = R_node_suc(R_node_suc ~= RUL_node); + swapnodes(RBI_node, RUL_node) + end + % swap RML and RLL if in wrong order + if get_lobe_node_xpos(2) > get_lobe_node_xpos(3) + % swap LML node and LLL node + swapnodes(lobe_origins_id(2), lobe_origins_id(3)) + end + % + + % focus on ordering left lobes + all_lobe_node_xpos = get_all_lobe_node_xpos(); + left_lobe_node_xpos = all_lobe_node_xpos(4:6); + % if LLL not in 1st pos, swap LLL and left BI + [~,I] = min(left_lobe_node_xpos); + if I ~= 1 + % swap LLL node and left BI node + LLL_node = lobe_origins_id(4); + L_node_suc = successors(g,L_node); + LBI_node = L_node_suc(L_node_suc ~= LLL_node); + swapnodes(LBI_node, LLL_node) + end + + % swap LML and LUL if in wrong order + if get_lobe_node_xpos(5) > get_lobe_node_xpos(6) + % swap LML node and LLL node + swapnodes(lobe_origins_id(5), lobe_origins_id(6)) + end + + % set XData + h.XData = X; + + function out = get_lobe_node_xpos(idx) + % get a lobes x position given index in lobe cell array. + lobe_edge = find(g.Edges.ID == lobe_origins_id(idx)); + lobe_node = g.Edges.EndNodes(lobe_edge,2); + out = X(lobe_node); + end + + function out = get_all_lobe_node_xpos() + % get all lobes x postion + out = zeros(length(lobes),1); + for ii = 1:length(lobes) + out(ii) = get_lobe_node_xpos(ii); + end + end + + function swapnodes(node_idx_1, node_idx_2) + % swaps nodes positions in x of two nodes and all + % subsequent nodes. + + % identify all nodes on 1 + v1 = dfsearch(g,node_idx_1); + + % identify all nodes on 2 + v2 = dfsearch(g,node_idx_2); + + % identify relative shift + xshift = X(node_idx_1) - X(node_idx_2); + + % execute relative shift. + X(v1) = X(v1) - xshift; + X(v2) = X(v2) + xshift; + + end + + end + + function node_table = ExportGraph(obj, path) + % export the airway network as a graph. If lobe classifications + % are present, they will be included. + % + + parse_filename_extension(path, '.csv'); + + % construct node table + id = [obj.tubes.ID]'; + xyz = cell(length(id),1); + edge = cell(length(id),1); + for ii = id' + % get xyz position from midpoint of tube + skelpoints = obj.tubes(ii).skelpoints; + midpoint = floor(length(skelpoints)/2); + point_pix = obj.I2S(skelpoints(midpoint)); + point_mm = point_pix .* obj.voxdim; + xyz{ii,1} = num2str(point_mm,'% .2f'); + % convert children nodes to string + try + children = num2str([obj.tubes(ii).children.ID]); + catch + children = ''; + end + try + parent = num2str(obj.tubes(ii).parent.ID); + catch + parent = ''; + end + % save edge + edge{ii,1} = [parent, ' ', children]; + end + node_table = table(id, edge, xyz); + + try + % add lobe classifications if exist + lobe = AirQuant.list_property({obj.tubes.region},'lobe')'; + node_table = addvars(node_table, lobe); + catch + end + + % export to csv + writetable(node_table, path) + end - end methods(Static) diff --git a/network/TubeNetwork.m b/network/TubeNetwork.m index 1af3866..be5a86c 100644 --- a/network/TubeNetwork.m +++ b/network/TubeNetwork.m @@ -51,6 +51,7 @@ Dmap tubemat tubepath + isloops end methods @@ -283,7 +284,7 @@ options.type = 'both' options.usesegcrop logical = false options.method char = 'linear' - options.gpu logical = 1 + options.gpu logical = false end for ii = progress(1:length(obj.tubes), 'Title', 'Making tube patches') @@ -355,7 +356,7 @@ if nargin < 2 method = 'topnode'; end - [g, glink, gnode] = skel_2_digraph(obj.skel, method); + [g, glink, gnode, obj.isloops] = skel_2_digraph(obj.skel, method); end function ge = TubesAsEdges(obj) @@ -363,7 +364,6 @@ % init edgetable gn = TubesAsNodes(obj); - asedges = gn.Edges; % add incomming edge to nodes that have no incoming edges nin = find(indegree(gn) == 0); @@ -376,7 +376,7 @@ % asedges.EndNodes(height(asedges)+1,:) = [max(asedges.EndNodes(:))+1 nini]; end - gn.Edges.ID = gn.Edges.EndNodes(:,2); + gn.Edges.ID = str2double(gn.Edges.EndNodes(:,2)); ge = gn; @@ -390,29 +390,41 @@ % * Scope for improving efficiency, variables that change % size. % + + %%% identify origin node of all graphs + origintubes = []; + % loop through every tube + for ii = 1:length(obj.tubes) + currenttube = obj.tubes(ii); + % for current tube, find the origin tube (tube with no parent) + while ~isempty(currenttube.parent) + currenttube = currenttube.parent; + end + if ~any(origintubes == currenttube) + origintubes = [origintubes, currenttube]; + end + end - % get all tubes with generation 0 - zerogentubes = find([obj.tubes(:).generation] == 0); - - % set up initial nodes with zero gen tubes g = digraph(); % DFS from each zero gen tube, creating edgetable. - for ii = 1:length(zerogentubes) - edgestosearch = obj.tubes(zerogentubes(ii)); - if isempty(edgestosearch.parent) && isempty(edgestosearch.children) - g = addnode(g,1); + for origintube = origintubes + % case of isolated tube + if isempty(origintube.parent) && isempty(origintube.children) + g = addnode(g,num2str(origintube.ID)); + continue end + % independent DFS search from each origin tube. + edgestosearch = origintube; while ~isempty(edgestosearch) current_tube = edgestosearch(1); - parent = current_tube.ID; + edgestosearch(1) = []; if ~isempty(current_tube.children) for childtube = current_tube.children - g = addedge(g, parent, childtube.ID); + g = addedge(g, num2str(current_tube.ID), num2str(childtube.ID)); edgestosearch = [edgestosearch, childtube]; end end - edgestosearch(1) = []; end end end @@ -434,7 +446,13 @@ assert(isa(tubefunc,"char") || isa(tubefunc,"string"), ... 'tubefunction must be provided as char/string.') for ii = progress(1:length(obj.tubes), 'Title', strcat('RunAllTubes: ', tubefunc)) + try obj.tubes(ii).(tubefunc)(varargin{:}); + catch e + warning([char(tubefunc), ' failed for tube index ', num2str(ii), ' .']) + fprintf(2,'error identifier :\n%s\n',e.identifier); + fprintf(2,'There was an error! The message was:\n%s\n',e.message); + end end end @@ -496,7 +514,11 @@ function Measure(obj, varargin) % get values for ii = 1:numtubes - label_all = [eval(object_eval).(labelname)]; + try + label_all = [eval(object_eval).(labelname)]; + catch % incase tube doesn't have the label + label_all = NaN; + end if ischar(label_all) outlabel{ii} = label_all; else @@ -511,6 +533,25 @@ function Measure(obj, varargin) end + function [statsopt, regionopt] = VisualisationOptions(obj) + % get visualisation options from :attr:`tubes` by :attr:`tube.Tube.stats` and :attr:`tube.Tube.region`. + statsopt = []; + regionopt = []; + + % get stats options + for ii = 1:length(obj.tubes) + if ~isempty(obj.tubes(ii).stats) + current_stats = fieldnames(obj.tubes(ii).stats); + statsopt = [statsopt, setdiff(current_stats,statsopt)]; + end + if ~isempty(obj.tubes(ii).region) + current_region = fieldnames(obj.tubes(ii).region); + regionopt = [regionopt, setdiff(current_region,regionopt)]; + end + end + + end + function [regionkwarg,regionid] = ParseRegion(obj, regionkwarg) % visualisation utility method to interpret the `region` % keyword argument @@ -570,7 +611,11 @@ function Measure(obj, varargin) % convert labels into numbers if isempty(options.types) - types = unique(vals); + if isfield(obj.regioncategories,options.name) + types = obj.regioncategories.(options.name); + else + types = unique(vals); + end else types = options.types; end @@ -632,7 +677,7 @@ function vol3daxes(obj, ax) % VISUALISATION - function h = Plot(obj, options) + function [h, ge] = Plot(obj, options) % Very powerful graph visualisation tool of tubes. Can % manipulate by label, edge weight and colour. % @@ -845,17 +890,17 @@ function Plot3D(obj, options) % colour: *OPTIONAL* `default = 'c'` color of surface. Can be any % MATLAB accepted color format, e.g. RGB. % type(char): *OPTIONAL* `default = 'seg'` can be either - % `'seg'` or `'skel'`. which to plot surface. - % region + % `'seg'` or `'skel'`. which to plot surface. + % smooth_sz(int): *OPTIONAL* `default = 0` size of gaussian + % smoothing kernel. 0 for no smoothing. % % .. note: - % plotting multiple patches (>20) significantly reduces - % performance. This is why this function tries to collate everys - % patch that is designated a different colour and then plots - % them. + % Plotting multiple patches (>20) significantly reduces + % performance. This function attempts to group + % segments of the same colour to plot together. % % .. warning: - % Choosing a colour option that is continuous will cause this + % Choosing a colour option with several colours will cause this % method to take a very long time. % % @@ -878,6 +923,7 @@ function Plot3D(obj, options) options.colour = '' options.colouridx = 1 options.type {mustBeMember(options.type,{'seg','skel'})} = 'seg' + options.smooth_sz = 0 end @@ -910,6 +956,8 @@ function Plot3D(obj, options) V(tubeii.([options.type,'points'])) = 1; end end + + % plot each color vol as isosurface individually @@ -917,16 +965,23 @@ function Plot3D(obj, options) for ii = 1:max(cdata) U = zeros(size(V)); U(V==ii) = 1; + % smoothing + if options.smooth_sz > 0 + U = smooth3(U,'gaussian',options.smooth_sz); + end patch(isosurface(U),... 'FaceAlpha', options.alpha,... 'FaceColor', rgb(ii,:),... 'EdgeColor', 'none'); hold on end - else + else + if options.smooth_sz > 0 + V = smooth3(V,'gaussian',options.smooth_sz); + end patch(isosurface(V),... 'FaceAlpha', options.alpha,... - 'FaceColor', 'k',... + 'FaceColor', 'c',... 'EdgeColor', 'none'); hold on end @@ -1079,7 +1134,7 @@ function PlotSpline(obj, options) end - function [h, out] = Histogram(obj, options) + function [h, out] = Histogram2(obj, options) % Plot histogram of :attr:`tubes`. % % @@ -1091,6 +1146,9 @@ function PlotSpline(obj, options) % % Behaviour untested for label containing NaNs. % + % .. todo:: + % + % Remove redundant variables. % % Args: % label = *OPTIONAL* `default = 'generation'` set variable to count. @@ -1100,14 +1158,16 @@ function PlotSpline(obj, options) % in order. % labelidx(scalar) = *OPTIONAL* `default = 1`. Index of % chosen property in `label`. - % print(bool) = *OPTIONAL* `default = false`. print frequency + % print(bool) = *OPTIONAL* `default = false`. Print frequency % table. + % exact(bool) = *OPTIONAL* `default = false`. Bin values by + % their exact number if label is numerical. % % % Example: % >>> run CA_base.m; % >>> figure; - % >>> AQnet.PlotGen(); + % >>> AQnet.Histogram(); % % arguments @@ -1115,6 +1175,7 @@ function PlotSpline(obj, options) options.label = 'generation' options.labelidx = 1 options.print = false + options.exact = false end % get index for each branch @@ -1123,25 +1184,39 @@ function PlotSpline(obj, options) % get unique values outlabel_unique = unique(outlabel); - + categoricalinput = 0; % set up cell for each unique val if isnumeric(outlabel_unique) rownames = compose('%d', outlabel_unique); rownames = cellstr(string(rownames)); % 2nd variable defines the right side of each bin only - [count, ~] = histcounts(outlabel, [outlabel_unique; outlabel_unique(end)+1]); else % must be categorical, convert + categoricalinput = 1; rownames = outlabel_unique; outlabel = categorical(outlabel); outlabel_unique = categorical(outlabel_unique); - [count, ~] = histcounts(outlabel); end + if ~options.exact || categoricalinput + % determine counts by binning algorithm if numeric + [count, edges] = histcounts(outlabel); + else + % use exact bins + [count, edges] = histcounts(outlabel, [outlabel_unique; outlabel_unique(end)+1]); + end + + if categoricalinput + % convert to categoric + edges = categorical(edges); + else + % remove right most edge if numerical. + edges = edges(1:end-1); + end % count for all vals % show figure result as bar chart - h = bar(outlabel_unique, count'); + h = bar(edges, count'); title(['Number of tubes per ', options.label]) xlabel(options.label) ylabel('count') @@ -1161,7 +1236,80 @@ function PlotSpline(obj, options) end end - + + function h = Histogram(obj, options) + % Plot histogram of :attr:`tubes`. + % + % + % .. note:: + % + % Consider case where label contains NaNs. + % + % .. warning:: + % + % Behaviour untested for label containing NaNs. + % + % .. todo:: + % + % Remove redundant variables. + % + % Args: + % label = *OPTIONAL* `default = 'generation'` set variable to count. + % if `char` must be an :class:`tube` property, :class:`tube` `stats` or `region` field. + % e.g. `'generation'`, `'arclength'`, `'lobe'`. + % Can also be vector of length equal to number of tubes + % in order. + % labelidx(scalar) = *OPTIONAL* `default = 1`. Index of + % chosen property in `label`. + % region(char) = *OPTIONAL* `default = ''`. Region to + % group label values by. + % + % + % Example: + % >>> run CA_base.m; + % >>> figure; + % >>> AQnet.Histogram2(); + % + % + arguments + obj + options.label = 'generation' + options.labelidx = 1 + options.region = '' + end + + % get index for each branch + outlabel = GetTubeValues(obj, options.label, options.labelidx); + outlabel = outlabel'; + + if ~isempty(options.region) + regions = GetTubeValues(obj, options.region, 1); + else + regions = cell(1, length(outlabel)); + regions(:) = {'all'}; + end + + % make any possible nans into chars + [regions{~cellfun(@ischar,regions)}] = deal('none'); + % convert to categorical for histogram function + regions = categorical(regions); + regions_unique = unique(regions); + + % count for all vals + % show figure result as bar chart + tcl = tiledlayout('flow'); + for iregion = regions_unique + nexttile + h = histogram(outlabel(regions == iregion)); + title(iregion) + xlabel(options.label) + ylabel('count') + end + title(tcl,options.label) + grid on + + + end % Data IO function ExportOrthoPatches(obj, path, casename, options) diff --git a/test/Test_workflow.m b/test/Test_workflow.m new file mode 100644 index 0000000..c95c701 --- /dev/null +++ b/test/Test_workflow.m @@ -0,0 +1,41 @@ +%% Main function to generate tests +function tests = Test_workflow() +tests = functiontests(localfunctions); +end + +%% fixtures +function setupOnce(testCase) + AQpath = AirQuantAddPath; + + % download test data if necessary + casename='chestct'; + AQdownload_data(casename); + + % set names + data_dir = fullfile(AQpath,'data','airquant'); + results_dir = fullfile(AQpath,'data','results'); + + % clear results dir if anything stored + if exist(fullfile(results_dir,casename),'dir') + rmdir(fullfile(results_dir,casename),'s') + end + + % get file names + filesuffix = '.nii.gz'; + sourcef = fullfile(data_dir,[casename,'_source',filesuffix]); + segf = fullfile(data_dir,[casename,'_airway',filesuffix]); + skelf = fullfile(data_dir,[casename,'_airway_PTKskel',filesuffix]); + + % save vars + testCase.TestData = struct('casename',casename,'sourcef',sourcef,'segf',segf,'skelf',skelf,'results_dir',results_dir); +end + +function test_runcases(testCase) + vars=testCase.TestData; + wf_clinicalairways_fwhmesl(vars.casename, vars.sourcef, vars.segf, vars.skelf, vars.results_dir); + % check AQnet saves + verifyTrue(testCase,exist(fullfile(vars.results_dir,vars.casename,[vars.casename,'_AQnet.mat']),'file')>0) +end + +function teardownOnce(testCase) +end diff --git a/test/runalltests.m b/test/runalltests.m index 1a3b717..a960d71 100644 --- a/test/runalltests.m +++ b/test/runalltests.m @@ -1,4 +1,5 @@ % AirQuant https://github.com/ashkanpakzad/AirQuant % Run all tests -runtests('TestSkel2Digraph.m') \ No newline at end of file +runtests('TestSkel2Digraph.m') +runtests('Test_workflow.m') \ No newline at end of file diff --git a/tube/Tube.m b/tube/Tube.m index 02652c2..af481a0 100644 --- a/tube/Tube.m +++ b/tube/Tube.m @@ -255,9 +255,14 @@ genname = [regiontype, '_gen']; currentube = obj; count = 0; - while ~strcmp(currentube.region.(regiontype), ... - currentube.parent.region.(regiontype))||... - isempty(currentube.parent) + % check same region +% if length(currentube.parent) > 1 + while ~isempty(currentube.parent) + current_reg = currentube.region.(regiontype); + parent_reg = currentube.parent.region.(regiontype); + if ~strcmp(current_reg, parent_reg) + break + end if length(currentube.parent) > 1 obj.region.(genname) = NaN; warning('Multiple parents, not possible to set generation') @@ -512,9 +517,10 @@ % i.e. the spline interval as found by :meth:`tube.Tube.FindSplinePoints` and saves them to :attr:`tube.Tube.source`. % % .. note:: - % This function will use the GPU if available only when - % :param:`options.method` = 'linear' and the `parallel copmputing toolbox`_ - % is installed. + % * GPU only compatible when :param:`options.method` = 'linear' and the `parallel copmputing toolbox`_ + % is installed. + % * GPU is not compatible with :param:`options.usesegcrop` = true. + % * GPU is a legacy feature and will be removed in future versions. % % .. todo:: % * Heavily reliant on the network class structure. @@ -541,7 +547,7 @@ options.usesegcrop logical = false options.method char = 'linear' options.sample_sz = obj.network.plane_sample_sz - options.gpu logical = true + options.gpu logical = false end assert(~isempty(obj.spline), 'spline is empty, see method MakeSpline') @@ -1396,17 +1402,24 @@ function PlotSpline(obj,options) % Data IO - function ExportOrthoPatches(obj, path, casename) + function ExportOrthoPatches(obj, tarpath, casename) % export perpendicular slice patches of this tube. % % export the perpendicular slice patches of this tube stored in - % source as int16 tiff files. + % source as int16 tiff files to a tar file. + % + % + % .. note: + % This uses the GNU tar command from the system to create the + % final tarball. Therefore it is only expected to work on unix system. + % This was the only way to achieve the desired 'append' + % functionality that MATLAB `tar` doesn't offer. % % % Args: - % path(str): path to directory to save the exported patches. - % The directory will be created if it doesn't already - % exist. + % tarpath(str): path to tarfile save the exported patches. + % The tarfile will be be appended to or created if it + % doesn't already exist. % casename(char): casename to use as prefix to each patch % slice name. % @@ -1414,11 +1427,6 @@ function ExportOrthoPatches(obj, path, casename) % >>> run CA_base.m; % >>> AQnet.tubes(98).ExportOrthoPatches('patches','example') % - - % make directory - if ~exist(path, 'dir') - mkdir(path) - end % ensure char casename = char(casename); @@ -1427,21 +1435,28 @@ function ExportOrthoPatches(obj, path, casename) allslices = 1:length(obj.source); chosenslices = obj.PruneMeasure(allslices); + % save files to temp dir + atemp_dir = tempname; + mkdir(atemp_dir); + + % check tar path + tarpath = parse_filename_extension(tarpath,'.tar'); + % loop through slices for k = chosenslices + % save as int16 float = obj.source{k,1}; img = int16(float); - - paddedk = sprintf('%08d', k); - % save as int16 TIF - imgsavename = fullfile(path, [ ... - casename, ... + % patch name + paddedk = sprintf('%08d', k); + patchname = [casename, ... '_id_',num2str(obj.ID), ... '_gen_', num2str(obj.generation), ... - '_slice_',paddedk, ... - '.tif']); + '_slice_',paddedk]; + % save as tif to temp dir + imgsavename = fullfile(atemp_dir, [patchname, '.tif']); imgdata = img; t = Tiff(imgsavename,'w'); @@ -1458,6 +1473,18 @@ function ExportOrthoPatches(obj, path, casename) write(t,imgdata) close(t); end + + if length(chosenslices) > 0 + % append files to tar path + tarcommand = ['tar --append -f ', tarpath, ' --remove-files -C ', atemp_dir,... + ' $(cd ', atemp_dir, ' ; echo *.tif)']; + + % execute tar + system(tarcommand); + end + + % remove temp dir + rmdir(atemp_dir,'s'); end function toGif(obj, filename, options) diff --git a/workflow/batch_clinicalairways_fwhmesl.m b/workflow/batch_clinicalairways_fwhmesl.m new file mode 100644 index 0000000..1271dcb --- /dev/null +++ b/workflow/batch_clinicalairways_fwhmesl.m @@ -0,0 +1,93 @@ +% Parent function to set up and run AQ on several cases based on given +% config file + +function batch_clinicalairways_fwhmesl(source_dir, seg_dir, skel_dir, overwrite) + +if nargin < 1 + % default behaviour is not to overwrite. + source_dir = 'source'; +end + +if nargin < 2 + seg_dir = 'airway'; +end + +if nargin < 3 + skel_dir = 'airway_skel_pakzad'; +end + +if nargin < 4 + % default behaviour is not to overwrite. + overwrite = 0; +end + +%% parse input +results_dir = 'results'; + +% get list of casenames from source dir +alldir = dir(fullfile(source_dir)); +alldir = alldir(~[alldir.isdir]); +raw_names = string({(alldir.name)}); +casenames = strrep(raw_names,'.nii.gz',''); + + +%% Set-up directories +% names input files, give fullpaths if not in matlab path +% add AirQuant library to path +try + AirQuantAddPath(); +catch + warning('Could not find AirQuantAddPath, try running it manually first.') +end + +mkdir_existok(results_dir) + +%% run loop +skippedcases = []; + +for ii = 1:length(casenames) + % init + casename = char(casenames(ii)); + skip = 0; % do not skip by default + disp([num2str(ii),' of ', num2str(length(casenames))]) + + % check if it already exists + case_dir = fullfile(results_dir,casename); + if overwrite == 0 && isfolder(case_dir) + disp(['[',casename,'] already exists. overwrite == false, therefore skipping.']) + fskip; + continue + end + + % Get filenames + filesuffix = '.nii.gz'; + sourcef = fullfile(source_dir,[casename,filesuffix]); + segf = fullfile(seg_dir,[casename,filesuffix]); + skelf = fullfile(skel_dir,[casename,filesuffix]); + try + skip = wf_clinicalairways_fwhmesl(casename, sourcef, segf, skelf, results_dir); + catch e %e is an MException struct + fprintf(2,'error identifier :\n%s\n',e.identifier); + fprintf(2,'There was an error! The message was:\n%s\n',e.message); + skip = 1; + end + + if skip == 1 + fskip; + end +end + +disp('The following cases were skipped:') +disp(skippedcases) + + function fskip + % if informed to skip the current case, it throws a warning if not + % already for the current case and saves the name to provide a full + % list at the end. + if skip == 0 + skippedcases = [skippedcases, string(casename)]; + warning('skipping case...') + end + skip = 1; + end +end \ No newline at end of file diff --git a/workflow/make_lobe_dataset.m b/workflow/make_lobe_dataset.m new file mode 100644 index 0000000..0e825c6 --- /dev/null +++ b/workflow/make_lobe_dataset.m @@ -0,0 +1,19 @@ +%% parse input +top_results_dir = 'results'; + +% get list of casenames from source dir +alldir = dir(fullfile(top_results_dir)); +% remove any that begin with '.' (i.e. '.', '..', and '.DS_Store') +alldir = alldir(~startsWith({alldir.name}, '.')); +casenames = string({(alldir.name)}); + +for ii = 1:length(casenames) + casename = casenames(ii); + disp(casename) + disp([num2str(ii), ' of ', num2str(length(casenames))]) + results_dir = fullfile(top_results_dir, casename); + % load data + load(fullfile(results_dir, casename + '_AQnet.mat')); + AQnet.ExportGraph(fullfile(results_dir, casename + '_graph.csv')); + clear AQnet +end \ No newline at end of file diff --git a/workflow/wf_clinicalairways_fwhmesl.m b/workflow/wf_clinicalairways_fwhmesl.m new file mode 100644 index 0000000..25e1ca2 --- /dev/null +++ b/workflow/wf_clinicalairways_fwhmesl.m @@ -0,0 +1,205 @@ +% Parent function to set up and run AQ on several cases based on given +% config file + +function [skip, runinfo] = wf_clinicalairways_fwhmesl(casename, sourcef, segf, skelf, root_results_dir) + + % prune ends of airways + prune_ends = [2, 2]; + + % recommended fwhmesl parameters + num_rays = 60; + ray_interval = 0.2; + +%% run loop + % init + runinfo = struct('casename',casename); + skip = 0; % do not skip by default + disp(['[',casename,'] ','Starting']) + results_dir = fullfile(root_results_dir,casename); + mkdir_existok(results_dir) + + % set up log + logname = fullfile(results_dir,[casename, '_log.txt']); + diary(logname) + disp(datetime) + + % savename is given to automatically save/load results. + savename = fullfile(results_dir, [casename, '_AQnet.mat']); + + % check if results already exists + if isfile(savename) + warning('AQ object already exists.') + end + + % Check if each datafile exists, if not it throws a warning and skips + % the current case. + checkfile(sourcef, 'source') + checkfile(segf, 'segmentation') + checkfile(skelf, 'skeleton') + + if skip == 0 + + tic; % start timer + + % Load CT data as double + meta = niftiinfo(sourcef); + CT = double(niftiread(meta)); + + % Load Airway segmentation and its skeleton as logicals + S = logical(niftiread(segf)); + skel = logical(niftiread(skelf)); + + % Initialise AirQuant + % Parses CT, segmentation and skeleton to compute airway tree graph + disp(['[',casename,'] ','Init AirQuant.']) + AQnet = ClinicalAirways(skel, source=CT, header=meta, seg=S,fillholes=1, ... + spline_sample_sz=0.1, plane_sample_sz=0.2, ... + max_plane_sz=40, largestCC=1,originmethod='carina'); % removed spline_sample_sz=0.5, plane_sample_sz=0.5 as gsk scans are huge, max_plane_sz=40, + + % save graph + AQnet.ExportGraph(fullfile(results_dir, [casename, '_graph.csv'])); + + %%% Generate initial analysis figures and save + % segskel + f = figure(); AQnet.Plot3D(alpha=0.3); hold on; AQnet.Plot3D(type='skel',alpha=1); + fig_save(f,fullfile(results_dir, strcat(casename,"_skel"))); + close(f); + + % plot3 + f = figure(); AQnet.Plot3D(alpha=0.3); hold on; AQnet.Plot3(); + fig_save(f,fullfile(results_dir, strcat(casename,"_plot3"))); + close(f); + + % splines + f = figure(); AQnet.Plot3D(alpha=0.3); hold on; AQnet.PlotSpline(); + fig_save(f,fullfile(results_dir, strcat(casename,"_spline"))); + close(f); + + % generation per lobe histogram + f = figure(); AQnet.Histogram(label='lobe_gen',region='lobe') + fig_save(f,fullfile(results_dir, strcat(casename,"_hist_lobegen"))); + close(f); + + % plot 2d + f = figure(); [h, g] = AQnet.Plot(colour='lobe',weightfactor=3); + try + AQnet.GraphLobarLayout(h, g) + catch + end + fig_save(f,fullfile(results_dir, strcat(casename,"_plot"))) + close(f); + + % lobe 3d + plot 3 + f = figure(); AQnet.Plot3D(colour='lobe'); hold on; AQnet.Plot3(colour='lobe'); + fig_save(f,fullfile(results_dir, strcat(casename,"_plot3d_lobe"))); + close(f); + + % gen 3d + f = figure(); AQnet.Plot3D(colour='generation'); + fig_save(f,fullfile(results_dir, strcat(casename,"_plot3d_generation"))); + close(f); + + % tortuosity + f = figure(); AQnet.Histogram(label='tortuosity',region='lobe') + fig_save(f,fullfile(results_dir, strcat(casename,"_hist_tortuosity"))); + close(f); + + % set prune + AQnet.RunAllTubes('SetPruneLength',prune_ends); + + % get patches + AQnet.MakeTubePatches(method='linear'); + + % make FWHM measurement + AQnet.Measure('AirwayFWHMesl', num_rays, ray_interval, 0); + + % save measures + save(savename, "AQnet"); + + % save measurements + AQnet.ExportCSV(fullfile(results_dir,strcat(casename,"_FWHMesl"))); + + % save orthopatches + tarpath = fullfile(results_dir,[casename, '_patches.tar']); + AQnet.RunAllTubes('ExportOrthoPatches',tarpath, casename); + % save grid preview of patches + grid_preview(tarpath, 8, 8, fullfile(results_dir,[casename, '_patch_preview.png'])); + + % generate post analysis figures + % plot 2d - avg D + f = figure(); [h, g] = AQnet.Plot(colour='lobe', weight='diameter_mean', ..., + weightfactor=10, label='lobe_gen'); + try + AQnet.GraphLobarLayout(h, g); + catch + end + fig_save(f,fullfile(results_dir, strcat(casename,"_plot_diameter"))); + close(f); + + % intertapering + f = figure(); AQnet.Histogram(label='intertaper',region='lobe'); + fig_save(f,fullfile(results_dir, strcat(casename,"_hist_inter"))); + close(f); + + % diameter + f = figure(); AQnet.Histogram(label='diameter_mean',region='lobe'); + fig_save(f,fullfile(results_dir, strcat(casename,"_hist_diameter"))); + close(f); + + % reset + disp(['Case: ', casename, ' complete.']) + disp(['Total time: ', num2str(toc/60), ' minutes.']) + disp(datetime) + + % save info to summary table + % total number of tubes + runinfo.n = length(AQnet.tubes); + % lengths + lengths = cell2mat(AirQuant.list_property({AQnet.tubes.stats},'arclength')); + runinfo.arclen_mm = sum(lengths); + % lumen volume + vols = cell2mat(AirQuant.list_property({AQnet.tubes.stats},'volume'))/1e6; + vols = vols(1,:); + runinfo.lumen_vol_l = sum(vols,'omitnan'); + % max lobe generation + try + gens = cell2mat(AirQuant.list_property({AQnet.tubes.region},'lobe_gen')); + runinfo.maxlobegen = max(gens); + runinfo.runtime_m = toc/60; + catch + warning([casename,' ---------------------- cell2mat converrsion FAILED because of missing lobe labels? CHECK! ----------------']) + runinfo.mesfailures = casename; + end + % breakdown by lobe + try + lobes = AirQuant.list_property({AQnet.tubes.region},'lobe'); + lobecodes = {'RUL','RML','RLL','LUL','LML','LLL'}; + for i = 1:length(lobecodes) + idx = strcmp(lobes,lobecodes{i}); + runinfo.([lobecodes{i},'_n']) = sum(idx); + runinfo.([lobecodes{i},'_len']) = sum(lengths(idx)); + runinfo.([lobecodes{i},'_vol']) = sum(vols(idx),'omitnan'); + runinfo.([lobecodes{i}, '_maxlobegen']) = max(gens(idx)); + end + catch + warning([casename,' lobe info not available']) + end + % append runinfo to csv + runinfo_T = struct2table(runinfo); + writetable(runinfo_T, fullfile(root_results_dir, 'summary.csv'), 'WriteMode', 'append'); + end + + + + close all; + diary off; + + function checkfile(filepath, filetype) + % checks if the file exists, if not it will throw a warning and + % inform the script to skip the current case. + if ~exist(filepath, 'file') + warning([filetype,' file at path: ', filepath, ' does not exist.']) + end + end + +end \ No newline at end of file