Skip to content

CT Airway interpolation

Jump to bottom
Ashkan Pakzad edited this page Nov 22, 2021 · 1 revision

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. 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 .

Clone this wiki locally