Royal Academy of Sciences New Zealand Open Science
Open Science

High quality statistical shape modelling of the human nasal cavity and applications

Published:

The human nose is a complex organ that shows large morphological variations and has many important functions. However, the relation between shape and function is not yet fully understood. In this work, we present a high quality statistical shape model of the human nose based on clinical CT data of 46 patients. A technique based on cylindrical parametrization was used to create a correspondence between the nasal shapes of the population. Applying principal component analysis on these corresponded nasal cavities resulted in an average nasal geometry and geometrical variations, known as principal components, present in the population with a high precision. The analysis led to 46 principal components, which account for 95% of the total geometrical variation captured. These variations are first discussed qualitatively, and the effect on the average nasal shape of the first five principal components is visualized. Hereafter, by using this statistical shape model, two application examples that lead to quantitative data are shown: nasal shape in function of age and gender, and a morphometric analysis of different anatomical regions. Shape models, as the one presented here, can help to get a better understanding of nasal shape and variation, and their relationship with demographic data.

1. Introduction

The human nose is an important and complex organ having many functions, including ventilation, filtration and olfaction. People who suffer from a nasal function impairment have a reduced overall quality of life [1]. Improving our knowledge of the key structures and key features of the nasal cavity is therefore important. The human nose shows large variation in shape between different individuals, for example variations of the bone structure [2,3] or temporal variations in the nasal cycle [4]. These variations alter the airflow through the nasal cavity, and can thereby affect the functioning of the nose [5]. To get a better understanding of normal and pathological functioning of the nasal cavity, it is important to study the effects of these variations. Nowadays, surface models of the nasal cavity are typically based on tomographic data coming from computed tomography (CT) scans. A whole range of studies have been conducted in the past, using data coming from a single or few patients [6,7], to studies where data of a much larger population is used [813]. The approach in this work is different from the previous in that it is not directly limited to the shapes obtained from the tomographic data. From this point of view, it is much more related to the work done by Gambaruto et al. [14] about nasal airflow modelling using the representation of a nasal passage based on a Fourier descriptor.

We propose to use the statistical shape modelling (SSM) technique to capture the shape variations [15]. In this way, we avoid a priori assumptions about morphology. SSM is a widely used tool to capture the morphological variability present in a population. A main application is their use as prior knowledge in model-based automated image segmentation. They were popularized for this purpose by Cootes et al. [16]. An overview of their use in medical image segmentation can also be found in Heimann & Meinzer [17]. In the present work, clinical CT data is used to capture the nose in 46 patients, and from this data a high quality statistical shape model is created. In literature, different techniques (e.g. level sets, Fourier descriptors) have already been used to create a shape model of the nose [1822]. In Liu et al. [20] and Sun et al. [22] image alignment and image processing of cross-sections are used to create a standardized geometrical model. In the latter work, a full upper respiratory airway model is created, containing less details of the nasal complex. A similar approach is followed in Nejati et al. [21], where cross-sections of the segmented mesh are taken and processed. Thinned representations of the nasal cavity are then computed, and a reference template is deformed such that it matches this thinned representation. Finally, all deformations are averaged to obtain an average geometry. They reported an increase in tolerance to a wider variety of nasal geometries than earlier methods. However, none of the methods discussed above use the direct three-dimensional information coming from training shapes. In Huang et al. [19] level sets are used to construct a shape model. Their model and the shape instances it produces for different positions along the PC axes show artefacts, i.e. left and right nasal channels are fused at the top of the nasal cavity, and different turbinates are merged together.

The technique used in the present work, based on cylindrical parametrization, is highly suitable for creating a nasal shape model because of the shape of the nasal cavity. More precisely, all nasal shapes of the population were closed at the post-nasal region in such a way that the resulting nasal cavity can be seen as a U-bended cylinder. For more information, the reader is referred to §2.2. We will show that this approach not only allows capturing and qualitatively reporting the natural anatomical variations present, but also more importantly, enables to relate shape variations with other features. In the following, this is called the quantitative analysis of the nasal cavity. In this paper, we show two quantitative applications of the SSM: we look at how age and gender are related with nasal shape variations, and we apply the shape model to study and discuss the influence of shape variation on the volume, surface and surface-to-volume ratio (SVR) of certain anatomical regions. These regions are the nasal valve and nasal vestibule region, the respiratory region, the olfactory region and the post-nasal region. The anatomy of these regions is important for the well-functioning of the nose, e.g. respiratory region for heating and humidifying the inhaled air. For sake of clarity, it is defined what is meant with ‘shape’ in this work. Strictly speaking, when two objects fall perfectly together after alignment (due to rotation and translation), they are said to have the same ‘form’. When these two objects are additionally allowed to have a different relative scale, they are said to have the same ‘shape’. In this work, however, we consistently use the term ‘shape’ instead of ‘form’, due to its connotation with ‘statistical shape model’, although we did not correct for scale in the alignment.

2. Material and methods

2.1. Data acquisition

Clinical cone-beam CT (CBCT) scans of six patients and helical CT scans of 40 patients were analysed. Data were obtained from patients with very different nasal or sinus related complaints, where a CT scan was requested by the treating ear nose throat (ENT) physician. The spatial resolution of the scans ranges from 250 to 430 µm perpendicular to the axial scan direction, and from 250 to 600 µm in the axial direction. The group included 29 females and 17 males. The average male subject age was 37 years (minimum age is 17, maximum is 60), and the average female subject age was 46 years (minimum age is 20, maximum is 84).

2.2. Individual model generation

First, the different steps are given that were followed to generate 46 different surfaces of the nasal cavity based on CT data. A classical and well-known approach was used in this work. For each CT scan a three-dimensional geometrical model was built by segmentation using Amira v. 6.0.1. On CT slices of the nose one can observe air, soft tissue and bone. Greyscale thresholding was used to determine the nasal passages, because it is the simplest way of segmentation and very suitable to find air spaces in the skull. This, however, led to inclusion of the sinuses, which are connected to the nasal cavity. Manual separation of these sinuses was overseen by an ENT physician. This separation, and by extension all steps described below, was done by the same operator for all 46 CT scans to avoid inter-operator generated variations between segmentations. At this stage and where necessary, some small corrections were made in the file containing the result of the segmentation, from now on called ‘label file’. Eventually, the goal is to represent each nasal shape with the same amount of points, where those points are located at the same anatomical position. This is done in the surface correspondence step, discussed below, but demands that all the training surfaces have the topology of a sphere (genus-0). To meet this requirement, holes in the label files were therefore filled such that the eventual geometrical surface extracted from the label file is genus-0. The field of view also varies for each scan, so that the end of the post-nasal region is different for each nasal shape. For this reason, the labels cannot be used as is. This would introduce a large non-anatomical variation in the shape model, as variation would be the extent of the region over which the scan was taken. Therefore, the labels were cut perpendicular to the axial scan direction at the floor of the inferior turbinates (figure 1).

Figure 1.

Figure 1. Segmentation result for a nose (side-view). Left top and bottom figures show the raw and smoothed nasal shape, respectively. The position of the four coronal slices (ad) on the right are indicated on the shape. The applied mask is shown in blue, and the segmentation labels in red.

The generated labels are used to build a triangulated surface geometry by applying the marching cubes algorithm [23]. As a next step, these surface models were smoothed, taking care that potential surface shrinkage was minimal. Algorithms exist that can be used for this (for example, the surface smoothing algorithm in Amira v. 6.0.1 is based on the paper ‘Curve and surface smoothing without shrinkage’ by Taubin [24]). At the end, the surface models were still quite rough when, for example, comparing them to those used in computational fluid dynamics. Nonetheless, further smoothing was unnecessary because of the smoothing effect of averaging when deriving shape model instances. The number of triangles was not decimated, because this is not necessary for the correspondence procedure.

All surfaces were exported in the stereo lithographic (STL) file format, and tiny holes were cut in both nostrils at the same anatomical position for all shapes. This was the final prerequisite of the correspondence algorithm. To increase the amount of morphological variation captured by the statistical shape model, mirrored versions of all nasal cavities are generated. Paraview v. 5.4.0 was used to generate these mirrored shapes, hereby doubling the size of the training data for the statistical shape model. The average shape obtained from principal component analysis (PCA) will be symmetrical because of this step.

2.3. Surface correspondence

Surface correspondence defines a dense one-to-one map between the surfaces in the population. In this paper, the correspondence method of Huysmans et al.[15] was modified to guarantee a uniform surface sampling. The method proceeds in four main steps. In the first step, each surface was cylindrically parametrized, using the method of Huysmans et al. [25], equipping the vertices of each surface with -coordinates. These coordinates define a low-distortion one-to-one mapping between each surface and the cylinder, where the small holes in the nostrils correspond to the two boundaries (ends) of the cylinder and where distortion is measured using the stretch metric of Sander et al. [26]. A correspondence between any two nasal cavities in the population can be derived by composition of these maps. In the second step, the surfaces were rigidly transformed and the parametrizations were rotated (-translation) for optimal alignment according to the minimum description length (MDL) objective [27]. The MDL is a measure based on the minimum description length principle: the sampled surfaces are coded in a message where the encoding is determined by the PCA model built from the correspondence. The length of the total message together with the encoded model, determines the quality of the model, and in this way the quality of the correspondence. Because of this, a trade-off exists between goodness-of-fit and model complexity. In this work the simplified MDL measure that was introduced by Thodberg [28], is used,

2.1

It is a function of the shape mode variances . The parameter  is set to be the expected variance of the noise in the data, such that the variation captured by all modes with a variation smaller than  is considered noise. The MDL favours compact models, this can be seen from the fact that goes to zero when all eigenvalues go to zero, where a lower value of indicates a better quality of correspondence. In the third step, the correspondences were further fine-tuned by additional non-rigid deformations in the -domain for each of the surfaces, measuring optimality again with the MDL objective. This optimization was done at three levels of detail. In the fourth and final step, the optimal parametrizations were used to generate a vertex correspondence for the population, where each shape is represented with a fixed number of landmarks that are at corresponding locations across the population. The term ‘landmarks’ here refers to Type III landmarks, defined by Bookstein [29] as mathematically derived points, which are different from Type I and II landmarks that are points at clearly identifiable features of the anatomy. The obtained correspondence is termed an anatomical correspondence in literature on image registration. This term is debatable, especially in inter-subject correspondences, where large differences in anatomy could be present. The obtained correspondence in this paper is the one that is optimal with respect to the MDL objective, commonly used in statistical shape modelling. As pointed out, the MDL objective measures the compactness of the PCA model calculated from the correspondences. Lower MDL values are assigned to correspondences that result in models with less variance. In this sense, models that only model variation in shape and do not include artificial variations like sliding of vertices along the surface are favoured. As it turns out, correspondences with optimal MDL tend to match up similarly shaped surface regions across the population, which are often also the same anatomical structures. See figure 2 for a visualization of the correspondence and cylinder domain for two subjects.

Figure 2.

Figure 2. A visualization of the correspondence for two subjects (A and B). The top two rows show the sagittal and axial view of the nasal cavity surface where the red and blue lines reveal the (u,v) coordinates of the map to the cylinder. The red lines are at constant v (around the cylinder) and the blue lines at constant u (along the cylinder). The lines are at corresponding locations between A and B. Third row: the cylinder domain is colour coded with the mean curvature of the corresponding surface points, revealing the correspondence to the surface on the cylinder. The red lines correspond to the ridges of the post-nasal region (centre of the cylinder) and the turbinates (towards both sides of the cylinder). A very similar map is obtained for A and B indicating a good correspondence. Bottom row: the cylinder domain with u and v iso-parameter lines. These lines correspond to the lines on the nasal cavity surfaces.

2.4. Statistical shape modelling

Using the correspondence calculated in the previous step, each shape is now represented by a set of 100 K corresponding landmark points distributed over the boundary surface of the nasal cavity and part of the post-nasal cavity. A SSM was built from the 92 corresponded nasal shapes by applying PCA on the vertices of the shapes. From a mathematical standpoint, a three-dimensional shape with n landmarks can be represented as a vector x with length 3n,

2.2

For all 46 training examples and 46 mirrored counterparts, such a vector was generated. An important step is to ensure that all the shapes are optimally aligned. A nasal shape was considered to be independent of position and orientation, but not of scale. To determine the optimal poses, the Procrustes analysis was used. In this analysis, all shapes are aligned with a reference shape. In the first iteration no average shape is yet calculated, and cannot be used as reference shape. For this reason, a shape from the training set is chosen as reference shape to start the analysis. Because no scaling had to be used, the analysis came down to two steps. In the first part, differences due to translation were removed by calculating the centroid of each shape and translating them to the origin. In the second part, quaternions were used to compute the rotation matrix between each shape of the training set and the reference shape.

In a next step, PCA was applied, delivering an average nasal shape (by averaging the corresponding landmarks) and an orthogonal set of shape variations (the so-called shape modes). Existing and new nasal shapes can then be described as the sum of the average nasal cavity and a specific linear combination of the shape modes. More formally, this is done by adapting the shape parameters , so that the point coordinates will be displaced and a new surface is created,

2.3

and

2.4

P contains the eigenvectors of the covariance matrix and b is a vector given by equation (2.4). It is common to limit the possible value of to the range ±3 times the standard deviation , so that the model generates shapes similar to those in the training set. As most of the variance is concentrated in the first few shape modes, in general SSM can provide a compact parametrization of the space of possible shapes. This is dependent, among other things, on the shape under consideration and the amount of training data.

In figure 3, the normalized cumulative sum of all  is visualized. In the beginning the curve shows a strong increase, which indicates that these shape modes are important to include in the model. More to the right, the curve flattens and adding new shape modes to the model will add very little new variance, or unwanted variance caused by noise on the original data. A (95%) cut-off has to be chosen, and the variance that is captured by the resulting model is called ‘explained variance’ in figure 3.

Figure 3.

Figure 3. Cumulative variance of the principal components of variation (compactness). The dashed line indicates the 95% limit and intersects with the cumulative curve when the number of shape modes is equal to 46.

2.5. Statistical analysis of shape

Because 95% of the total variation captured by the shape model is described by 46 shape modes, only these modes are considered in the following. This allows expressing the shapes of the training population in a useful way for further statistical analysis. More specifically, there are 46 values for each, with i running from 1 to 46. It is purely coincidence that the amount of shapes of the training population is equal to the cut-off number for the parameters ; there is no mathematical constrict enforcing this.

The total population was divided according to gender. Next, for each parameter a simple linear regression model to age was created using R v. 3.4.2. This meant that in the end 46 simple regression models are obtained relating a specific with age. Normality was checked using a graphical technique called ‘normal probability plot’. In such a plot, the values for a parameter are plotted against age, and deviations from a straight line suggest deviation from normality. In case of deviation, a Box–Cox power transformation was applied. This is a procedure to find an appropriate value for the exponent λ to use to transform the original data into normal distributed data. The value of λ represents the power to which all data should be raised. To accomplish this, the Box–Cox power transformation searches for the best value for λ in a certain predefined region. The Box–Cox transformation function that was used has the form:

2.5

where and are the transformed and the original response variable, respectively. Outlier detection was based on Cook's distance, which measures for each data point the change of the regression fit with and without the presence of that specific point. In this way, the impact of each individual data point can be assessed and outliers can be detected. Robust regression is applied in the case of an outlier where the effect could not be neglected. The t-value coming from this robust regression is used to determine if the null hypothesis can be rejected or not (correlation or not). All regression models were tested at a 5% significance level.

This allows creating a visual map of how the shape of the nasal cavity is correlated with age. For the sake of clarity only the regression models of the female nasal shapes will be considered; the procedure is identical for the male shapes. As input, a lower age of 20, and an upper age of 80, is given to the regression model. The model outputs . For the models with p < 0.05, the corresponding entries are subtracted, giving . The values that showed no significant correlation with age were set to zero. As a final step, one can calculate by using equation (2.3).

Other approaches similar to the one above exist, e.g. linear discriminant analysis, but for the sake of brevity we will not go deeper into them.

For the relation between nasal shape and gender a similar approach to the one above was used. For each corresponding male and female set of , a two-sample t-test or a Wilcoxon test was used; the latter in case of non-normally distributed data. Significance was tested at a 5% level. Those values that showed no significant correlation with gender were set to zero, and was calculated.

2.6. Volume partitioning

A simple approach was followed to partition the total volume of each shape. The average nasal cavity was divided into four anatomical regions (figure 4ac) based on [3033]. The division was done using a bounding box in Paraview, and by selecting those triangular faces on the surface that comprise the different volumes of interest. This puts a constraint on the way the anatomical regions can be defined. A trade-off exists; lowering the box will at the same time add unwanted parts to the region. Therefore, the regions definition was overseen by an ENT physician and the current division made an optimal compromise. The nasal valve and nasal vestibule region compose a first region (green, NVVR). The respiratory region is shown in (blue, RR), the olfactory region on top in (yellow, OR) and posteriorly the post-nasal region in (pink, PNR). Because of the way the different regions are defined, the nasal valve (separating the NVVR and RR) is implicitly represented by a plane. In reality, the nasal valve has a non-planar shape. However, in literature it is not uncommon to make such an approximation because it allows a precise definition, e.g. for modelling purposes [32]. A fairly good agreement exists when comparing the anatomical regions with definitions in existing literature [3133]. Small differences are not unusual because of a difference in nasal shape, and variation also exists between literature. For example, the definition of the olfactory region used in this work, has a better agreement with Croce [31] than with Netter [33]. As a last step the NVVR, OR and PNR were closed by defining a plane for each region connecting the boundary triangles, i.e. those triangles with neighbours at only one side. The tiny holes made for the correspondence step were also closed at this point. These steps made it possible to define a volume for all anatomical regions.

Figure 4.

Figure 4. (a) The different anatomical regions: nasal valve and vestibule region (NVVR), olfactory region (OR), respiratory region (RR), post-nasal region (PNR). (b) Birds-eye view of nasal shape, where the RR is not visualized, showing the opening of the other regions. (c) Bottom view of the nose. In light grey, the regions that were subtracted from the NVVR and PNR are clearly visible. This subtraction is needed for a correct surface area calculation. (d,e) The result of the first principal component (PC) on the average shape with equal to  and , respectively. It is clear that the regions defined before, vary accordingly.

Finally, the volume of the respiratory region was calculated by subtracting the volume of the total nasal shape from that of the other three previously mentioned. The ten first shape modes are consecutively applied to the average shape with ranging from  to  in steps of  (see figure 4d,e). The effect of each shape mode on the surface area and on the volume of the different anatomical regions is captured and compared.

3. Results

3.1. Qualitative visualization

The different shape modes can be visualized by applying them to the average nasal shape using equation (2.3). In figure 5 this is shown for the first five modes, with values for ranging from  to . For each PC, the colour map represents the magnitude of variation. Regular PCA is applied, which results in global modes, meaning that these modes act on the entire mean nasal surface.

Figure 5.

Figure 5. Results of applying the different PCs to the average shape, with different values for . Each time the colour map represents the magnitude (in mm) of variation for .

The first eigenmode applies a general scaling to average nasal shape. It is not a perfect isotropic scaling, but nonetheless there is a clear size increase/decrease in all three dimensions. This could be expected because in the alignment of the shapes only rotation and translation were used and no scaling. The second mode predominantly results in the contraction or expansion of the nasal shape along the anterior–posterior axis, and can be clearly observed in the region behind the choana. The expansion of the shape along the anterior–posterior axis goes hand in hand with a, be it smaller, expansion in the perpendicular plane, making the nose wider. The third eigenmode has a large effect on the curvature of the line connecting the anterior tip of the nose and the end of the post-nasal region. More precisely: the anterior tip of the nose drops and the region behind the choana flattens with decreasing values for . The fourth mode affects the width of the nasal shape. In addition, the posterior part of the olfactory region increases in height. It seems that the length of the olfactory region relative to that of the total nasal length is changed by this mode. The fifth eigenmode is the first one in line that has a non-negligible effect on the symmetry aspect of the nasal model. Applying this mode to the average nasal shape, which is symmetrical, results in bending the anterior part of the nose aside. One could say it reflects a specific type of nasal septum deviation. The fifth mode of the nasal shape is visualized from a front view, as the side view reveals little shape variation.

3.2. Quantitative applications

Figure 6 shows the result of the statistical analysis described in the ‘Material and methods’ section. The analysis gave and as a result. This is visualized by projecting the Euclidean distance map on the average nasal shape. This Euclidean distance is calculated from the vectors using equation (2.3). This allows creating a visual map of how the shape of the nasal cavity is correlated with age and gender.

Figure 6.

Figure 6. (a,b) Euclidean distance map between young and old nose, projected on the average nose, for males and females. These two maps show the magnitude and location of shape variation in function of age. (c) Euclidean distance map between average male and female nose, projected on the average nose. This map shows the magnitude and location of shape variation in function of gender.

3.3. Morphometric analysis of anatomical regions

3.3.1. Nasal volume

As a second application example for the nasal shape model, the surface, volume and their ratio was studied for the different anatomical regions under the influence of each of the first 10 shape modes. These first 10 modes represent 69% of the total variation captured by the shape model. Not more are studied because of the extensive manual work involved per eigenmode.