 Research Article
 Open Access
 Published:
A deformable template method for describing and averaging the anatomical variation of the human nasal cavity
BMC Medical Imaging volume 16, Article number: 55 (2016)
Abstract
Background
Understanding airflow through human airways is of importance in drug delivery and development of assisted breathing methods. In this work, we focus on development of a new method to obtain an averaged upper airway geometry from computed tomography (CT) scans of many individuals. This geometry can be used for air flow simulation. We examine the geometry resulting from a data set consisting of 26 airway scans. The methods used to achieve this include nasal cavity segmentation and a deformable template matching procedure.
Methods
The method uses CT scans of the nasal cavity of individuals to obtain a segmented mesh, and coronal crosssections of this segmented mesh are taken. The crosssections are processed to extract the nasal cavity, and then thinned (‘skeletonized’) representations of the airways are computed. A reference template is then deformed such that it lies on this thinned representation. The average of these deformations is used to obtain the average geometry. Our procedure tolerates a wider variety of nasal cavity geometries than earlier methods.
Results
To assess the averaging method, key landmark points on each of the input scans as well as the output average geometry are located and compared with one another, showing good agreement. In addition, the crosssectional area (CSA) profile of the nasal cavities of the input scans and average geometry are also computed, showing that the CSA of the average model falls within the variation of the population.
Conclusions
The use of a deformable template method for aligning and averaging the nasal cavity provides an improved, detailed geometry that is unavailable without using deformation.
Background
Simulation of airflow through human upper airways is of importance in a number of medical applications, including the improvement of artificial respiratory devices [1] for conditions such as obstructive sleep apnea (OSA), and also for delivery of drugs through aerosol deposition [2]. The first step in simulating the flow is producing a geometrical shape. The nasal cavity has two roughly symmetric passages (divided by a thin plate of bone and cartilage called the nasal septum) with each comprised of at least three distinct meatuses (Fig. 1) [3]. Significant variation between individuals is present, such as developmental variations of the underlying bone structure [4, 5] and also the temporal variations that may occur due to the nasal cycle [6]. Deviation of the nasal septum is common [7]. Typically, studies have considered airway shapes obtained from cadaver casts or tomography scans of a single individual [2, 8, 9] or otherwise have studied separately the airflow through the geometries obtained from of a small number (<10) of individuals [7, 10]. They lack universality and may produce quite different airflow results due to the observed variation in nasal cavities of individuals [7].
In a study by Liu et al. [11], nasal cavity geometries from 30 healthy individuals were obtained and used to produce a ‘standardized’ nasal cavity. The procedure involved first segmenting sequences of xray computed tomography (CT) scan data from a number of individuals from the transverse plane, then aligning and scaling the individual geometries appropriately, and finally superimposing the geometries and taking mean/medians of each image. This method lacks the ability to work with the various shape deviations that are typically present; for instance it cannot work with deviated nasal septa so only nasal cavities with straight septa were considered. We observe that this does not adequately capture the variations in the population. The method used by Liu et al. also had difficulties coping with congestion (the blockage of the nasal passages due to membranes lining the nose becoming swollen from inflamed blood vessels, which is quite common and difficult to distinguish from tissue/bone in CT scans), and due to the nature of the approach much of the shape information in the superior nasal cavity section was lost.
A somewhat different approach was presented in [12], where the average CT images were used to produce a single geometry. That is, the result is the segmentation of the average scan images, rather than the average of the segmentations of the scan images. However, this approach is also limited to nasal cavities with similar shapes (taken from ethnically uniform females of similar age) and it is not evident how to extend this method to general nasal cavity scans.
In [13, 14], airway shapes were decomposed into a small set of morphological feature coefficients which are then directly averaged. Two feature descriptions were studied: Fourier descriptors and descriptions based on structural decomposition of the medial axis or ‘skeleton’ of coronal nasal cavity crosssections. The Fourier descriptor method relies on representing finitelength single closed curves as periodic coordinatevalued functions. Applying a lowpass filter on such a function produces a ‘simplified’ shape. There are numerous limitations to applying the Fourier descriptor method. Firstly, for obtaining a mean geometry, a common rotational reference frame must be specified, or else crosssections must be aligned rotationally. The method also requires modification to be applicable to crosssections that are not single closed curves i.e. where congestion is present, making a single closed curve insufficient to represent the boundary of the crosssection (see Fig. 2). It is also not known how well the Fourier descriptor method works with shapes that show significant variation.
In this work, we study the problem of producing a more representative airway geometry. By this we mean an airway geometry where the size and shape are averaged across airway geometries of a sample of the population. The method presented can accept more shape variation than previouslypublished methods (see Fig. 3). The contribution of this paper is thus a new shape registration method for the human nasal cavity, with the following advantages:

1.
It can succeed with a wider variety of nasal cavity shapes compared to previous methods, including those with different numbers of conchae (3 or 4), deformed septa, congestion, or conchae with differing forms of curvature.

2.
It provides a simple way to describe these geometrical variations, allowing for statistical analysis.

3.
It is insensitive to variations in nasal passageway crosssectional area. For instance, during the nasal cycle, or during periods of congestion.
We use a method based on the medial axis description, which consists of finding the ‘skeletons’ or ‘thinned’ representations of airway crosssections. For instance, a thick vertical slit would be represented by a vertical line passing through the center of the slit. We then use a different type of averaging procedure that is robust and can process geometries with large variations. The method is based on pointset registration [15], in which the resulting skeletons are subdivided into a large set of points and the correspondences between these points and the points from a ‘reference’ or ‘template’ skeleton are discovered.
The most similar previous work was that of [11]. The primary difference is that no deformable registration was used for alignment of nasal cavities to the templates. We show that the use of deformable registration allows better matching between nasal cavities, and results in an average geometry that includes finer details (see Fig. 4).
Methods
Patient CT data sets
Anonymized head and neck CT scans from 26 patients were used to create patient specific airway models. The data set included a mix of male and female (9 male, 17 female) subjects aged from 17 to 70 years, with an average age of 52.7 years and a standard deviation of 13 years. The obtained scans were retrospective and anonymized. No pathologies were observed in the airways. Patients were imaged awake, in a supine position. Scans were comprised of stacks of horizontal slices with 0.43 ×0.43 mm planar resolution; slices were spaced 0.6 mm apart vertically.
Airway segmentation
The segmentation procedure was as follows: First, the airway from the tip of the nose to the trachea was segmented from other structures in the CT images using 3D Slicer software (v4.3.1) [16]. A thresholding procedure was used for the segmentation of air with the threshold varied from 1000 to 400 Hounsfield units (HU); see Fig. 5. The optimal threshold values varied from scan to scan and within a scan for areas with poor contrast or resolution. The latter was used for the interface between the nasal cavity and paranasal sinuses interface due to the small scale and complex geometry of the nasal cavity passages and sinus drainage openings. Manual slicebyslice editing was used for these problematic areas to obtain anatomically accurate airway segmentation.
In the next step, the segmented image data from all three image planes (axial, coronal and sagittal) was used to create a 3D surface model of the airway using 3D Slicer’s builtin marching cubes algorithmbased ModelMaker module [17]. The procedure was as follows: We applied a default ‘ThresholdEffect’, followed by a ‘PaintEffect’, and then the ‘Unsmooth’ ModelMaker function to generate the surface mesh. The generated model was then exported in an triangular mesh format for editing. MeshLab software (v1.3.3) [18] was then used to remove the sinuses and to smooth the surface of the model (see Fig. 6), using the holefilling ‘Replace’ module (Smooth mean value coordinates (MVC) method), and the localised robust smoothing, refinement submodule (default parameters: refinement: 0, reduction: 5.0, smoothing: 5.0). The frontal, maxillary, ethmoidal and sphenoidal sinuses were removed from the nasal cavity by manual editing in MeshMixer software (v10.6.53)^{1}. To smooth the airway model surface, while preserving its topology and size, a twostep MeshLab Laplacian Smoothing was applied (1D Boundary Smoothing, cotangent weight method). This procedure resulted in smooth 3D representations of the airway of each patient.
There is a large ‘air’ volume in front of the face which is topologically connected to the air inside the nasal cavity. We remove this volume via a morphological opening of the segmented volume [19], in the anterior volume, by a spherical structuring element of radius r. An additional erosion of radius r _{ ε } is used to capture ‘outlier’ areas that spherical structuring elements do not fit e.g. the narrow channel between the nose and the upper lip. The radius r must be set to a value that is larger than the internal walltowall dimensions of the nasal cavity but smaller than the volume outside the face. We used a value of r = 10 mm.
3D mesh model alignment
The absolute coordinates of the origin and relative rotation of each individual varies when CT scans are acquired. The first step in the analysis is the alignment of the airways of scans of different subjects so that they are all approximately in the same spatial position and have the same alignment and scale. This is because we ‘slice’ the scan into crosssections later, and we require the crosssections to approximately contain the same features in the structure. For alignment, the common approach is to identify a set of ‘landmarks’ on the shape. Landmarks must be “homologous anatomical loci that provide adequate coverage of the morphology, and can be found repeatedly and reliably” [20]. Specifically, landmarks are zerodimensional points, and we use the term ‘landmark’ as opposed to ‘feature’ for these points, both to make this distinction and for consistency with the field of geometric morphometrics.
Alignment of meshes according to shape similarity requires specification of seven degrees of freedom corresponding to three translation, three rotation, and one scale degree of freedom [11]. We locate a pair of landmarks in each scan. Bonebased landmarks are preferable to tissuebased landmarks as these are less subject to variation based on experimental conditions and nasal cycle phase. In [11], the following landmarks were selected for this purpose:

The tip of the anterior maxillary spine (AMS), where the nasal septum and maxilla bone intersect.

The choana; the beginning of the nasopharynx marked by the point where the two interior nostrils (choanae) meet. The two choanae are the openings between the nasal cavity and the nasopharynx, where the left and right nasal cavity passageways meet.
These landmarks are illustrated in Fig. 7 and specify six degrees of freedom. The final degree of freedom is rotation about the line between the landmarks. In [11], an automatic alignment procedure was used to specify this but there is no need to do this for our method since it is insensitive to rotation about this axis. We simply align the airways such that the AMS is at the origin and the line between the two landmarks is parallel to the yaxis (the zaxis is in the upwards direction towards the top of the head, and the xaxis is in the sideways direction). Additionally, once this alignment has been done, the geometry is uniformly scaled such that the choana lies at a distance of 60 mm from the AMS. This is based on the typical average length of the AMSChoana distance; see, for example [11].
Anatomical description of the individual nasal cavities
Coronal crosssections of the nasal cavity are extracted. This converts the 3D geometry to a set of 2D crosssections. These crosssections are processed independently. This is a common approach to studying nasal cavity shape since twodimensional shape analysis techniques are highly developed [11, 12, 21]. For instance, we can robustly identify the inferior and middle sections of the airway [14]. Sections are shown in Fig. 8. The major sources of variability in the airway shapes are:

Overall dimensions e.g. size of the nasal cavity (which is dependent on the size of the individual as well as other factors).

Crosssectional thickness of the passages due to congestion, the nasal cycle, or inflammation.

Height of the supreme meatus due to a lack of CT scan data near the orbital plane.

The curvature of the nasal septum. Variation is common even in healthy individuals [22].
We make the observation that, despite the variability in individuals, many airway shapes can be viewed as ‘warped’ or smoothly deformed versions of one another. Such variations are common in biological contexts [23–25]. Even though large variations are present, the various structures (meatuses, etc.) have approximately the same position and alignment with respect to one another. Thus we describe the crosssections by first obtaining the medial axis (e.g. Fig. 9).
The medial axis of a shape is the set of all points having more than one closest point on the object’s boundary. In practical terms, the medial axis provides a ‘thinning’ or ‘skeleton’ for a shape. Here the medial axis is computed according to the ZhangSuen thinning algorithm [26], which provides several desirable properties, for example each pixel on an edge has exactly two neighbors, and each pixel at a bifurcation point has exactly three. The branching of the airway into various meatuses can be seen consistently with this medial axis transformation [14].
The use of the medial axis transform has the benefit of eliminating the variability due to the crosssectional area of the airway. The result of the medial axis transformation is a set of curves. We describe each curve by choosing a set of points that approximate the curve. These points are known as semilandmarks. The distinction between semilandmarks and landmarks is that the relative spatial position of landmarks important, but the relative spatial position of semilandmarks is not  we can slide semilandmarks along curves without changing the curve that the semilandmarks define. In other words, if two sets of semilandmarks have different coordinates but define similar curves, we consider the sets of landmarks to be aligned with respect to one another.
Now we envision some ‘reference’ shape template (set of curves) that undergoes a smooth nonintersecting deformation to superimpose on the curves of a desired ‘target’ shape, based on a similarity metric. The deformation process and similarity metric are given in section ‘Matching Procedure’.
In what follows we address the problem of how to deform the reference template to fit target shapes. We desire deformations that are as ‘simple’ as possible in the sense of not requiring complex, convoluted deformations. This is to preserve the spatial relationship of adjacent structures. Also, as we are interested in the position of landmarks, it is desirable to have a deformation formulation that is based on altering the position of (semi)landmarks [20]. That is, we prefer a deformation that can be represented by deforming a small set of ’sample’ points and ’extrapolating’ this deformation to the rest of the space.
To summarize, we require a deformation formulation that has the following features:

Produces a continuous and smooth deformation that matches the control points

Invariant to rotation and translation

Nonaffine; this allows us to take into account nasal septum deviations

Based on landmark positions
Thinplate splines (TPS) [27] satisfy all of these requirements, and are formulated as follows. Let x _{ i } be a set of control points which are mapped to the points y _{ i }. The thinplate spline is the unique deformation f that minimizes the following function,
The second term in this equation is the squared curvature of f. λ is a realvalued stiffness parameter that controls the relative weight of stiffness vs. accurate point matching. This deformation f can be found exactly via radial basis functions [28].
Matching Procedure
The first step of our method is the medial axis transformation, as described. The second step is nonrigid registration of the thinned crosssections to a ‘template’ consisting of an airway with demarcated subregions (using thinplate splines), using thinplate splines (see Fig. 10). The template itself is based on the CarletonCivic airway geometry [11] and is shown in Fig. 11. Even though the CarletonCivic geometry lacks precision in terms of capturing airway features, for this stage of the process only an approximate skeleton of the airway is required, since we later discard the template’s control point coordinates and use coordinates from the input data.
For finding the optimal deformation, we define the following sets of points:

The reference point set, which in this case is the set of points comprising the medial axis of the CarletonCivic geometry.

The target point set, which is the set of points comprising the medial axis of an individual airway section.

The control points, which are a small set of points which represent our TPS deformation.
For each pair of slices, we deform the set of control points, use TPS to extrapolate this deformation to the other points, and try to find the deformation that will optimize the alignment of the target and reference points (see Fig. 10).
Let the set of control points in the reference template be {X _{ i,j }}, where j indexes the crosssectional image number and i indexes each control point in a crosssectional image. After nonrigid registration, we obtain the points {Y _{ i,k,j }} where k indexes the target (not reference) crosssection. We compute the average position of each control point \(\bar {Y}_{i,j}=\frac {1}{K}\sum _{k=1}^{K}Y_{i,k,l}\) to obtain a new set of control points. We then align all of the crosssectional images to these coordinates, to obtain a superimposed image where most of the airways overlap. Finally, the median intensity of this image is determined. This completes the description of our algorithm, however two details remain: The pairwise matching algorithm, and the method used for selecting control points. For the pairwise matching procedure mentioned above, we use the method of deformable robust point matching (RPM), see [28] for further details.
We select control points as follows. There must be enough control points to compactly represent the natural range of deformations that are observed. In theory, if a curve in the reference template can be approximated by a piecewiselinear curve, then the nodes of this curve can be taken as control points and would provide a nearlycomplete representation of any possible deformation of the curve that preserves the piecewiselinear structure. Thus, we use the RamerDouglasPuecker (RDP) algorithm [29, 30] on the skeleton to produce a set of control points. The RDP algorithm requires a parameter ε representing the characteristic length scale of the resulting approximation. By varying ε, the number of control points can be varied, with larger ε leading to fewer control points, as shown in Fig. 12.
Our final procedure for obtaining the mean geometry is as follows:

1.
Align to landmarks (AMS & Choana)

2.
Slice into crosssections

3.
Medial axis transformation on crosssectional images

4.
For each image:

1.
Find the thinplate spline deformation that matches reference skeleton to target skeleton

2.
Discard coordinate information of reference control points and use average coordinates of target control points

3.
Deform all target images (not skeletons) to average control point coordinates

4.
Average all deformed images, filter, and take the median to result in a median image mask

5.
Optionally, use the skeleton of the median image mask as a reference skeleton, and repeat steps 1)4).
To visualize this procedure, a flowchart is presented in Fig. 13.
Perpassage CSA calculation
We calculate the crosssectional area (CSA) along the crosssectional images of the airway, and we also calculate the CSA for the various parts (inferior, middle, superior) of the airway. To do this, the different parts are labelled separately and the number of pixels with each label counted. See Fig. 14.
For labeling the airway sections separately, an automated procedure is used. The traditional anatomical classification of the various sections is illdefined. For instance, the regions encompassed by the inferior turbinate and main passage overlap. Here, labeling is done based on the skeleton. After localization of the main branch point (the point where the middle turbinate branches from the main passage), the three skeletal branches from this point are cut off at distance r = 3 mm (to remove any ambiguity of the skeleton in the neighborhood of the branch point), and then the ‘inferior’ branch is taken to be the branch with most negative minimum in the topbottom direction. The ‘middle’ branch point is taken to be the branch with the most negative minimum in the leftright direction (assuming the right passageway; for the left passageway this is reversed).
After this, the procedure is to label all other parts of the skeleton. We must allow for ‘broken’ skeletons as well as multiple extraneous branches. This is done via the following process: consider all nonlabeled curves in the skeleton, and pick the one with an endpoint closest to any labeled curve in the skeleton. Then label this curve with the same label. This process is repeated until all curves in the skeleton have been labeled.
The final part of the procedure consists of labeling all pixels in the binary mask that are not in the skeleton. Care must be taken here as we cannot simply assign labels based on proximity to nearest curve in the skeleton. This would give incorrect labels in situations where, for instance, a ‘thick’ curve is in close proximity to a ‘thin’ one. The labeling procedure is similar to the one for curves except it is done based on pixels, i.e. at each step we pick the nonlabeled pixel that is closest to the set of labeled pixels and we assign this pixel the same label. This process is repeated until all pixels have been labeled. For the purposes of computational efficiency, an optimized nearestneighbor search is employed (Fast Library for Approximate Nearest Neighbors (FLANN)) [31].
Geometry results
The final result of this procedure is a standardized geometry. This geometry will be referred to as the UoAUC (University of Auckland  University of Canterbury) standardized nasal model in the figures. The resulting shapes are shown in Fig. 15. The full set of crosssectional slices is available for download (Additional file 1).
To objectively measure the performance of the deformation procedure, we varied λ, the ‘stiffness’ parameter of the deformation procedure. For ε, the characteristic length of the RDP algorithm which controls the number of control points, the base value of ε=1 was used. After deforming all crosssectional images, we then calculated how well the images overlapped with each other. We did this using the following metric:
where S _{ i } is the deformed image for the scan indexed by i.
In Fig. 16, from the AMS to the choana, the crosssectional image shape becomes more complicated and the overlap decreases. Also, the lowest λ value performs poorly for the first crosssections, but performs well for the latter crosssections. This may be because the anterior crosssections have fewer control points (being of simpler shape), with better skeletal alignment at the expense of crosssection alignment. Thus different values of λ seem to be optimal for different crosssectional positions. However, the plot illustrates that a very high λ value, corresponding to a ’stiff’, nondeforming transformation, is generally suboptimal. This indicates that the method performs better when deformable registration is used.
To measure the performance of the deformation procedure in matching features for each of the individual scans we manually determine if the TPSbased matching procedure has correctly aligned the position of various features of the airway with the template. To do this, we check if the deformed point coordinates are superimposed on the reference point coordinates. For curves, we check if their deformed control points are superimposed on the reference points. The features that were considered are:

The inferior/middle meatus branch point i.e. the point on each crosssection where the middle meatus and the inferior meatus join (points indicated ’1’ in Fig. 12).

The curve of the inferior meatus i.e. whether the inferior meatus in the template is properly matched to the inferior meatus in the test slice (black region in Fig. 14).

The middle meatus curve (green region in Fig. 14).

The passageway above the middle meatus, including the superior meatus (blue region in Fig. 14).
Each of these four features were considered for both the left and right passageways, resulting in 8 features that could either be properly or improperly matched in each individual scan. We confirmed the proper matching of these features at the control points along the curves; in the nasal valve region the deformation procedure achieved 100 % matching accuracy. For the nasopharynx region, we also found that the outline of the airway is matched in every scan. For the turbinate region, results varied in more posterior sections (no superior meatus) and more anterior sections, thus we show the results separately in Table 1.
Some scans have less variation and thus match better than other scans. In Fig. 17 we show histograms of matching performance, this time on a perindividual basis. For the majority of individual scans, the method successfully matches all 8 features; however for others some or all of the features fail to properly match. The ‘peak’ near 0 is higher than would be expected if each of the individual features failed to match at a rate independent of the others fail to match any features. This suggests that those features that aren’t matched belong to scans that have outlying variation.
To assess the new geometry, we compare the crosssectional area of each ‘part’ of the geometry (inferior/middle/posterior meatuses) with the individual scans (Fig. 14; procedure detailed in §2.5). Each of the separate passageways (inferior, middle, superior), in turn, start at a low value and then gradually peak before decreasing again (Fig. 18). Near the posterior of the cavity, the inferior turbinate section once again increases in crosssectional area. The measurements from the mean geometry are well within the natural range of variation and, most importantly, reflect the trends of changing crosssectional area (CSA) at different positions along the cavity.
Results and discussion
The values of both total CSA and persection CSA show variation among subjects. A possible reason for this is the nasal cycle. The CSA of the average geometry closely tracks the CSA of the test subjects. In other studies [11, 21] only total CSAs are compared, not persection CSAs. The mean geometry shows less percrosssection variation of CSA, especially in the middle meatus section, than the individual subjects. In addition, the mean geometry shows a ‘dip’ in the CSA of the superior meatus near the posterior of the geometry. A possible explanation could be the lack of overlap of the different scans in this section.
We also compare the CSAs of the individual sub jects with CSA results reported in the literature (Fig. 19) [32–34]. It is important to note that due to different alignment systems used in the literature, the CSA values may be offset by some amount. In general, there is less variation in the middle section and the variation increases towards the anterior. An increase in total CSA is apparent towards the anterior. A large increase happens in the nasal valve area (up to 3 cm away from the AMS landmark) and then the CSA remains approximately constant, before a dramatic increase again in the nasopharynx area. CSA variation is more apparent than the rigid nasal cavity as the nasopharynx area is surrounded by more tissue and muscle.
We anticipate the use of the methods outlined in this paper for analysis of differences in nasal cavity geometries in different populations. We expect that our model will lead to more accurate simulations of e.g. assisted breathing using nasal cannulae, drug deposition studies, and also simplify studies of the dynamics of natural breathing, due to the finer precision in reflecting the CSAs of the various sections of the nasal cavity geometries of the population.
Conclusions
A novel method has been developed for obtaining mean nasal cavity geometries from individual planar CT scans. A skeletonization procedure was used to find the general shape of the airway of each crosssectional image (and to correct for nasal cycle variations) and then a deformation procedure was used to align all the images to a common reference. The use of deformation is the unique aspect of this method that has not been used in previous studies in the literature. We have compared the similarity coefficient of the method with and without deformation, which shows a marked improvement in overlap of the images when deformation is used. This validates use of deformation for obtaining mean geometries. The geometry developed in this work is of finer resolution and higher quality than other, similar geometries that are currently available.
Endnote
^{1} Autodesk MeshMixer 3.0. 2015.
Abbreviations
 AMS:

Anterior maxillary spine
 CSA:

Crosssectional area
 CT:

Computed Tomography
 HU:

Hounsfield units
 MRI:

(Nuclear) Magnetic resonance imaging
 MVC:

Mean value coordinates
 OSA:

Obstructive sleep apnea
 RDP:

RamerDouglasPuecker algorithm
 RPM:

Robust point matching
 TPS:

Thinplate spline
References
 1
Spence CJT, Buchmann NA, Jermy MC. Unsteady flow in the nasal cavity with high flow therapy measured by stereoscopic PIV. Exp Fluids. 2012; 52(3):569–79.
 2
Xi J, Si X, Kim JW, Berlinski A. Simulation of airflow and aerosol deposition in the nasal cavity of a 5yearold child. J Aerosol Sci. 2011; 42(3):156–73.
 3
Lang J. Clinical anatomy of the nose, nasal cavity and paranasal sinuses. Stuttgart: Thieme; 1989.
 4
Noback ML, Harvati K, Spoor F. Climaterelated variation of the human nasal cavity. Am J Phys Anthropol. 2011; 145:599–614.
 5
Franciscus RG, Long JC. Variation in human nasal height and breadth. Am J Phys Anthropol. 1991; 85:419–27.
 6
Eccles R. A role for the nasal cycle in respiratory defence. Eur Respir J. 1996; 9:371–6.
 7
Churchill SE, Shackelford LL, Georgi JN, Black MT. Morphological variation and airflow dynamics in the human nose. Am J Hum Biol. 2004; 16(6):625–38.
 8
Kim SK, Chung SK. An investigation on airflow in disordered nasal cavity and its corrected models by tomographic PIV. Meas Sci Technol. 2004; 15:1090–6.
 9
Stringer NM, Cater JE, EatonEvans J, White C. Numerical Comparison of Air Flow Patterns in the Upper Airways of Adults and Neonates. In: Proceedings of the 17th Australasian Fluid Mechanics Conference. Auckland, New Zealand: The Faculty of Engineering in association with the Centre for Continuing Education The University of Auckland: 2010.
 10
Schroeter JD, Garcia GJM, Kimbell JS. A computational fluid dynamics approach to assess interhuman variability in hydrogen sulfide nasal dosimetry. Inhal Toxicol. 2010; 22(4):277–86.
 11
Liu Y, Johnson MR, Matida Ea, Kherani S, Marsan J. Creation of a standardized geometry of the human nasal cavity. J Appl Physiol. 2009; 106:784–95.
 12
Lee CF, Abdullah MZ, Ahmad KA, Lutfi Shuaib I. Standardization of Malaysian adult female nasal cavity. Comput Math Methods Med. 2013. article ID:519071.
 13
Gambaruto AM, Taylor DJ, Doorly DJ. Modelling nasal air flow using a Fourier descriptor representation of geometry. Int J Numer Methods Fluids. 2009; 59:1259–83.
 14
Gambaruto AM, Taylor DJ, Doorly DJ. Decomposition and description of the nasal cavity form. Ann Biomed Eng. 2012; 40(5):1142–59.
 15
Myronenko A, Song X. Point set registration: Coherent point drift. IEEE Trans Pattern Anal Mach Intell. 2010; 32(12):2262–75.
 16
Fedorov A, Beichel R, KalpathyCramer J, Finet J, FillionRobin JC, Pujol S, et al.3D Slicer as an image computing platform for the Quantitative Imaging Network. Magn Reson Imaging. 2012; 30:1323–41.
 17
Lorensen WE, Cline HE. Marching cubes: A high resolution 3D surface construction algorithm. In: ACM SIGGRAPH Computer Graphics. New York: ACM: 1987. p. 163–9.
 18
Cignoni P, Corsini M, Ranzuglia G. Meshlab: an opensource 3d mesh processing system. Ercim News. 2008; 73:45–46.
 19
Haralick RM, Shapiro LG, Vol. 1. Computer and Robot Vision. Boston: AddisonWesley; 1992.
 20
Zelditch ML, Swiderski DL, Sheets HD. Geometric Morphometrics for Biologists: A Primer, 2nd ed. San Diego: Academic Press (Elsevier); 2012.
 21
Lee CF, Abdullah MZ, Ahmad KA, Shuaib IL. Analytical Comparisons of Standardized Nasal Cavity. J Med Imaging Health Inf. 2014; 4(1):14–20.
 22
Holton NE, Yokley TR, Figueroa A. Nasal septal and craniofacial form in European and Africanderived populations. J Anat. 2012; 221:263–74.
 23
Swiderski DL. Morphological evolution of the scapula in tree squirrels, chipmunks, and ground squirrels (Sciuridae): an analysis using thinplate splines. Evolution. 1993; 47(6):1854–73.
 24
Adams DC, Rohlf FJ, Slice DE. Geometric morphometrics: Ten years of progress following the ‘revolution’. Ital J Zool. 2004; 71:5–16.
 25
Kleisner K, Chvȧtalovȧ V, Flegr J. Perceived intelligence is associated with measured intelligence in men but not women. PloS ONE. 2014; e81237:9.
 26
Zhang TY, Suen CY. A fast parallel algorithm for thinning digital patterns. Image Process Comput Vision. 1984; 27(3):236–9.
 27
Duchon J. Splines minimizing rotationinvariant seminorms in Sobolev spaces. In: Constructive Theory of Functions of Several Variables. Berlin, Heidelberg: Springer: 1977. p. 85–100.
 28
Chui H, Rangarajan A. A new point matching algorithm for nonrigid registration. Comp Vision Image Underst. 2003; 89:114–41.
 29
Ramer U. An iterative procedure for the polygonal approximation of plane curves. Comput Graphics Image Process. 1972; 1:244–56.
 30
Douglas DH, Peucker TK. Algorithms for the Reduction of the Number of Points Required to Represent a Digitized Line or its Caricature. Can Cartographer. 1973;10(11222).
 31
Muja M, Lowe DG. Fast Approximate Nearest Neighbors with Automatic Algorithm Configuration. In: International Conference on Computer Vision Theory and Applications (VISAPP ’09). Lisboa: VISIGRAPP: 2009. p. 1–10.
 32
Terheyden H, Maune S, Mertens J, Hilberg O. Acoustic rhinometry: validation by threedimensionally reconstructed computer tomographic scans. J Appl Physiol. 2000; 89:1013–21.
 33
Cheng YS, Yeh HC, Guilmette RA, Simpson QS, Cheng KH, Swift DL. Nasal deposition of ultrafine particles in human volunteers and its relationship to airway geometry. Aerosol Sci Technol. 1996; 25:274–91.
 34
Corey JP, Gunfor A, Nelson R, Fredberg J, Lai V. A comparison of the nasal crosssectional areas and volumes obtained with acoustic rhinometry and magnetic resonance imaging. Otolaryngol Head Neck Surg. 1997; 117:349–54.
Funding
This research was supported by a New Zealand Ministry of Business, Innovation & Employment (MBIE) grant (project number 3706026).
Availability of data and materials
A public release of the source data (CT scans of individuals) has not been made available at this time.
Authors’ contributions
AN developed the method for geometry standardization used in the paper and drafted the manuscript. NK carried out the segmentation of the CT scans and contributed to the ‘Methods’ section of the manuscript. JEC participated in the design and coordination of the study, collection of the CT scan data, and helped to draft the manuscript. MCJ participated in the collection of the CT scan data and segmentation of the scans. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
In accordance with the 2014 New Zealand Health and Disability Ethics Committees (HDEC) standard operation procedures (SOPs) guide, we did not require ethics approval for use of anonymized scans.
Author information
Additional file
Additional file 1
Crosssections for final UoAUC geometry. File is a HDF5 data set containing crosssectional slices for computed mean geometry. (H5 347kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Nejati, A., Kabaliuk, N., Jermy, M.C. et al. A deformable template method for describing and averaging the anatomical variation of the human nasal cavity. BMC Med Imaging 16, 55 (2016). https://doi.org/10.1186/s1288001601548
Received:
Accepted:
Published:
Keywords
 Human airways
 Image registration
 Thinplate splines