Validated reconstructions of geometries of nasal cavities from CT scans

Developing validated objective measures of nasal airflow is paramount for improving the management of nasal obstruction. Nasal airflow can be studied objectively by simulating the flow in the computer. This requires a faithful reconstruction of the nasal geometry since the simulated airflow is sensitive to small variations of the complex shape of the nasal cavity. We here show that altering the geometry by less than 1 mm can change the airflow two-fold. We also show that a faithful reconstruction of the nasal geometry is possible with a threshold-based segmentation. Utilizing the known geometry of a CT phantom, we determine an optimal segmentation threshold of −450 HU. Changing this threshold by 100 HU alters the geometry by only about 0.1 mm. We use this verified segmentation to extract nasal geometries of three patients and simulate the respective airflows using the Lattice Boltzmann method. Using a simple model, we can predict how the reconstruction threshold affects the resistance to airflow. Since the segmentation and the simulation can be automated completely, this is an important step toward an objective analysis of nasal airflow based on CT scans.


Introduction
The geometry of the nasal cavity critically affects nasal airflow and thus nasal function (Elad et al 2008). This is evident in cases of nasal obstruction, which have a prevalence of around 30% and have been shown to have a negative impact on both disease-specific and global quality of life (Fuller et al 2017). However, the accuracy of the diagnosis and treatment of nasal obstruction is compromised by the lack of an accurate objective measure of nasal airflow (Pawar et al 2010). Measuring the patient's airflow directly often only provides gross values and depends on the patient's effort (Ottaviano and Fokkens 2016). Alternatively, the nasal cavity can be reconstructed from CT scans in the computer and the airflow can be simulated The study of nasal airflow using computational fluid dynamics (CFD) can be broken down into several steps: (i) Imaging: CT or MRI scans containing the complete nasal cavity need to be acquired. (ii) Segmentation: the surface outlining the nasal cavity needs to be extracted from the image data. (iii) Meshing: a computation mesh needs to be created, which defines the region and the resolution of the computer simulation. (iv) Simulation: the equations describing the airflow need to be solved on the mesh. (v) Analysis: the resulting flow data needs to be analyzed to extract useful information. Steps (iii)-(v) are necessary for every CFD simulations and have thus been studied extensively. Similarly, CT scans are routinely acquired and standard protocols have been established. Consequently, the crucial step in analyzing nasal cavities is the segmentation, which needs to be tailored to this particular geometry (Kim et al 2013). Moreover, to be useful in clinical settings, it is important that all steps combined work robustly and automatically, also in the hands of inexperienced users.
Computational fluid dynamics (CFD) simulations of the airflow in the nasal cavity have been carried out using various algorithms (Doorly et al 2008, Quadrio et al 2014. Most of these simulations solve the Navier-Stokes equations, which describe the dynamics of the velocity field of the airflow, using approximate methods, like laminar simulations, Reynolds averaged simulations, Large-Eddy simulation, or various turbulence models (e.g., k-ò or k-ω). In these methods, the computation mesh needs to conform to the boundary, which often requires additional smoothing of the surface ( (Gambaruto et al 2009). This smoothing alters the surface, degrades the effective resolution, and might thus introduce additional errors in the simulation (Gambaruto et al 2009). Methods that do not rely on solving the Navier-Stokes equations include Multiparticle collision dynamics (Gompper et al 2009) and the Lattice Boltzmann method (LBM) (Succi 2001). The advantage of these methods is the straight-forward implementation of boundary conditions, which do not require additional smoothing of the surface. For these reasons, the LBM has been used to simulate nasal airflow before (Finck et al 2007, Lintermann et al 2013a, 2013b. The reliability of the CFD simulations critically depends on a faithful reconstruction of the nasal geometry. Here, the key step is the image segmentation, i.e., the identification of the regions belonging to the nasal cavity in the grayscale 3D images obtained from CT or MRI scans (Sharma and Aggarwal 2010). The 3D data could in principle be analyzed by first segmenting 2D slices and then combining these into a single 3D geometry. The segmentation of 2D slices can be controlled manually, but the subsequent combination step is challenging, in particular in complex geometries like the nasal cavity. Consequently, it is advantageous to perform the segmentation using the 3D data directly. Typical segmentation techniques are similar for 2D and 3D images (Smistad et al 2015): threshold-based methods discriminate regions based on a single grayscale value. These methods, typically used in conjunction with the wide-spread marching cubes algorithm (Lorensen and Cline 1987), are computationally robust, but the choice of the threshold is crucial. Other methods, like region growing (Adams and Bischof 1994), clustering, watershed methods, and graph-based methods, are not based on thresholds but instead often require manually supplied seed positions of the regions that need to be segmented. Conversely, threshold-based methods can be automated and are thus used more frequently. However, it is not clear what threshold provides the most faithful reconstruction of the geometry (Sezgin and Sankur 2004). One possibility is to determine the segmentation threshold automatically, e.g., using Otsu's method (Otsu 1979), but such methods have often been developed for general segmentation problems and might thus not provide the best result in the specific case of nasal cavities. Alternatively, specific algorithms for finding an optimal threshold have been designed (Nakano et al 2013), but it is not clear how robust they are and how such algorithms behave in challenging situations.
Geometries of nasal cavities have been obtained using various segmentation techniques in the literature. Typical procedures vary widely, although many authors do not provide any details on the segmentation step beyond the name of the program that was The aim of the present paper is two-fold. First, we seek an objective protocol to derive the best segmentation parameters based on the known geometry of a CT phantom. Second, we want to understand how sensitively the flow simulations in realistic nasal cavities depend on the choice of these parameters. While CT phantoms have been used in the past to study reconstruction algorithms (Zeiberg et al 1996, Nakano et al 2013, the fine geometric details that appear in the nasal cavity have been neglected. We put an emphasis on these structures and will thus be able to determine how well such fine structures can be reconstructed. Moreover, we will be able to predict how the reconstruction changes when the parameters are modified.

CT phantom
Two modules from an AAPM CT Performance Phantom (CIRS Model 610) were scanned in this study. The first module (outside diameter: 7.5″ length: 2.5″ is an acrylic test object (radiodensity 110 K 130 HU) with eight rows of air-filled channels of square shape with side lengths 1.75, 1.5, 1.25, 1.00, 0.75, 0.61, 0.50, and 0.40 mm. The rows are 5 mm apart from each other and contain each 5 channels separated by a distance equal to their side length. The second module (outside diameter: 8″; length: 1.18″) is made of a proprietary epoxy with radiodensity 6 K 10 HU above water. It contains multiple air-filled, circular channels with diameters ranging from 2.5 to 7.5 mm, in 0.5 mm steps, and channel separations equal to the channel diameter.

Image acquisition
The phantom was scanned on a GE Discovery 750 HD scanner with our standard sinus stealth protocol (120 kV, rotation time: 0.8 s, pitch: 0.531, collimation: 20 mm, noise index: 3.85). The tube current was set to 90 mA, corresponding to a Volume CT Dose Index (CTDIvol) of 27.06 mGy (16 cm phantom). Images were reconstructed with 20% Adaptive Statistical Iterative Reconstruction (ASIR) and three different reconstruction protocols: 'bw' (bone kernel, slice thickness (ST) 2.5 mm), 'soft' (soft kernel, ST 1.25 mm), and 'skinny' (bone kernel, ST 0.625 mm). The Display Field of View (DFOV) was set to 20 cm, resulting in an in-plane resolution of h≈0.4 mm.
The patient scans were obtained from the Rhinoplasty Data Repository at the Massachusetts Eye and Ear Infirmary (MEEI). This repository is approved by the institutional review board (IRB) at MEEI and all patients in the data repository completed an electronic consent form. The same scanner as for the phantom was used for patient 3, while patient 1 and 2 were scanned on a Siemens Sensation 40 scanner. All scans were obtained using our standard sinus stealth protocol. Images were reconstructed with a slice thickness of 2 mm and two different reconstruction protocols: 'bw' (Siemens scanner: H70h kernel; GE scanner: bone kernel, 20% ASIR) and 'soft' (Siemens scanner: H30s kernel; GE scanner: soft kernel, 20% ASIR). The DFOV was set to about 18 cm, resulting in an in-plane resolution of h≈0.35 mm.

Geometry reconstruction
3D surface meshes of the scanned geometries were reconstructed in the computer using the marching cubes algorithm implemented by the 'GrayscaleModel-Maker' plugin of Slicer 3D (Fedorov et al 2012). This algorithm uses linear interpolation to extract a mesh representing the isosurface of a given radiodensityΘ (Lorensen and Cline 1987). The geometry was isolated and augmented with boundary conditions using custom python code calling paraview (Ayachit 2015) and Meshlab (Cignoni et al 2008). The segmentation, geometry processing, and CFD simulation were automated using custom python code.

Computational fluid dynamics
We simulate the airflow in nasal cavities using the Lattice Boltzmann method (LBM) implemented by Muphy (Bernaschi et al 2009). The LBM utilizes a uniform Cartesian grid, which we obtain by determining which cells of a regular grid lie inside the surface mesh. Muphy utilizes a compact representation of the simulation domain, such that only the mesh nodes that represent the fluid and the surrounding boundary are explicitly taken into account. This strategy results in a large saving in computational time and computer memory. In addition, Muphy is very flexible in simulating the airflow by leveraging either CPU or GPU-based hardware, thus being an ideal tool for high through-put simulation of the nasal airflow. Since this paper focuses on reconstructing the geometry, we only considered a grid resolution of 0.25 mm, which allowed us to simulate steady exhalation on a single GPU by keeping the box placed around the external nose at ambient pressure and imposing a volumetric flux of 25 cm 3 /s at the throat. This low flow rate was chosen since it simplifies the simulation. However, this is not a principle limitation of the LBM, which can also be used to simulate turbulent and oscillatory flow

Reconstructing CT phantoms
We determine the optimal reconstruction algorithm using an object of known geometry and radiodensity, known as a CT phantom. Here we chose a phantom with radiodensity and geometric properties similar to those of nasal cavities. We scanned the CT phantom and reconstructed its geometry as described in the Methods section. Figure 1 shows the reconstructed shape in an axial cross-section. The size of the channels is mainly affected by the threshold levelΘ used for reconstruction, but also by the parameters of the CT scan and the subsequent image reconstruction. The aim of this section is to determine what parameters yield the most faithful reconstruction of the phantom geometry. Figure 1 shows that the reconstructed shapes are often too irregular to measure the side lengths or diameters directly. However, we can always determine the crosssectional areaA and the perimeterΓ for each channel in each cross-section. From this, we obtain two estimates of the side lengtha s of the small, square channels

Quantifying reconstruction quality
Similarly, we get two estimates for the diametera l of the large, circular channels, Using these definitions, we measure the sizes a for the small and the large channels in many slices. This allows us to determine a mean channel sizeā as well as the associated standard deviationσ a . Figure 1 shows these values for a reconstruction with the 'bw' protocol at the threshold Θ=−425 HU. The dotted lines in the background indicate the channel sizesâ i s andâ i l for the small and the large channels that are given in the Methods section. The figure reveals that the typical error of the size estimate is smaller when it is based on the area. This is not surprising since small fluctuations of the channel shape will lead to larger fluctuations in perimeter than in area. Consequently, we exclusively use the size measurement based on the channel area in the following.
To estimate the overall error of the reconstruction, we sum the squared deviations of all channels and normalize them by the total numberN of channels, This value denotes the average deviation of the predicted channel sizeā from the true valueâ. The deviationδa depends on the reconstruction algorithm and the thresholdΘ.

Optimal reconstruction algorithm
The quality of a reconstruction of the CT phantom can be quantified by the associated deviationδa, which is affected by two contributions: first, there is the systematic error of choosing the wrong threshold for reconstructing the geometry. Second, there are inherent fluctuations due to variability in the production of the phantom, physical limitations of the image acquisition, and numerical uncertainties in the subsequent analysis. By minimizing both errors, we can determine the optimal parameters that lead to minimalδa and thus the most faithful reconstruction of the geometry.
The two error sources have different effects on the deviationδa. While the inherent fluctuations can be modeled as random noiseη, the systematic error depends on the reconstruction method. In the marching cubes algorithm that we use (Lorensen and Cline 1987), the location of the isocontour is determined using a linear interpolation between the values given at the discrete pixel locations of the CT image, see figure 2. Consequently, the deviation in the reconstructed position due to a sub-optimal threshold Θ is given by a simple linear function. This linear approximation is valid if the variation in Θ is smaller than the intensity difference between adjacent pixels, which is certainly the case for the features with high contrast that we are considering. Taken together, the measured channel sizeā i can thus be modeled as whereâ i is the true channel size, the second term models the systematic error, and the last term captures fluctuations. Here, χ denotes the sensitivity ofā i to changes in the threshold away from the optimal value * Q . Moreover, we consider uncorrelated white noiseη i of strengthω with h h w d á ñ = i j ij 2 and h á ñ = 0 i . Since the deviationδa is given by an average of -â a i i , see equation (3), the associated expected value d á ñ a follows from an ensemble average over the noise, where c, * Q , and ω are unknown parameters that can be determined by fitting equation (5) to measured values of δa as a function of Q.
We determinedδa for several reconstruction thresholdsQ and the three different reconstructions protocols 'bw', 'soft', and 'skinny' described in the Methods section. Figure 3 shows that the theoretical model explains the data very well. We thus fit equation (5) to the data to determine the optimal threshold * Q for each algorithm, see table 1. Moreover, figure 3 shows that the 'soft' reconstruction protocol is worst and leads to the largest noise strength ω. This is likely caused by too much smoothing of the data. Conversely, the 'bw' and 'skinny' protocol perform virtually identically for the large channels and show only small differences for the small channels. For  both protocols, the optimal segmentation threshold is close to * Q » -440 HU. The fit parameters also indicate that the intrinsic noiseω is on the order of 0.02 mm and thus much smaller than the CT resolution of h≈0.4 mm. Moreover, we find that the sensitivity of the reconstruction is roughlyc » -10 3 mm HU . Consequently, the position of the wall moves by 0.1 mm when the threshold is changed by 100 HU. Note that this value depends on the in-plane resolu-tionh of the CT scan. Since h is the only relevant length scale in this problem, dimensional analysis implies that the sensitivityc is proportional to h. The wall position would thus move by 0.25 mm after chan-gingQ by 100 HU at the lower resolution of h≈1 mm. Taken together, we thus identified an optimal threshold of * Q » -HU 440 and determined the sensitivity of the reconstruction to this choice.
The optimal threshold can be compared to typical intensities of the solid material of the CT phantom and the air inside it. Table 2 shows the corresponding radiodensities, which we measured manually using Horos (Horos Project 2018). Note that the arithmetic mean of the densities in the solid and in the air is close to the optimal threshold value obtained above. In fact, the reconstruction procedure shown schematically in figure 2 already suggests that a threshold value of the arithmetic mean between the intensity values of the two plateaus will lead to a good result and compensate for partial volume averaging. Our optimized reconstruction provides evidence that this heuristic works well.
In summary, we propose that two parameters are important for a faithful reconstruction of thin geometries with large contrast, like nasal cavities. First, the reconstruction kernel needs to preserve the edges. In our case, the 'bone' kernel, which is used in the 'bw' and 'skinny' protocol, works best, while the 'soft' kernel seems to introduce too much smoothing, thus leading to significant noise. Second, the intensity threshold on which the segmentation algorithm is based should be chosen such that it lies close to the arithmetic mean of the intensities of the two regions that need to be separated.

Reconstructing realistic nasal cavities
We next use the principles obtained from the analysis of the CT phantom to reconstruct geometries of realistic nasal cavities and simulate the airflow within them.

Automatic extraction of the nasal geometry from CT scans
To check whether our algorithm works on realistic data, we reconstructed nasal geometries from CT scans that were already present in the Massachusetts Eye and Ear Infirmary (MEEI) Rhinoplasty Data Repository. These CT scans contained image stacks created with the 'soft' and the 'bw' protocol along the axial and coronal axes of the patient. Based on the analysis presented in the previous section, we only considered the series using the 'bw' protocol, assuming that the 'soft' protocol would work less well, similar to the phantom reconstruction.
To see whether the axial or coronal CT scans are more suitable for reconstruction, we analyzed the resolution of the scans. All scans contain most of the head and since the images were limited to 512 columns and 512 rows, the in-plane resolution was in the range of = ¼ ( ) h 0.34 0.36 mm. Conversely, the slice thickness was set to 2 mm, so the in-plane resolution of the CT scans is higher than the resolution along the CT image stack. To maximize the quality of the reconstruction, it is thus important to choose the axis of the stack adequately. In essence, the resolution in each direction needs to be good enough to resolve the finest structural variations in this direction. The geometry of the nose exhibits the smallest variation in the coronal direction, i.e., the predominant direction of the flow. Conversely, the finest details, and thus the larger variations, can be seen in the coronal plane, perpendicular to the coronal axis. Consequently, the axis of the CT stack, where the resolution is worst, should be aligned with the coronal axis. Based on these qualitative arguments, we focused on reconstructing coronal scans exclusively.
We determine the optimal threshold valueQ for the reconstruction based on the principle that it should be close to the arithmetic mean of the radiodensity of the airspace defining the nasal cavity and that of the surrounding tissue. Estimating the respective intensities using Horos (Horos Project 2018), we find significant variations of both intensities, see table 3. The variations are likely patient-specific, although a different scanner was used for patient 3 than for the other two. Note also, that these measurements were done manually and might thus be subjective. Moreover, the intensities varied from one location to the next, which is captured in the uncertainty presented in table 3. These few examples show that the expected threshold might vary significantly from patient to patient. However, manual determination of the tissue intensities for each patient is not practical for an automated analysis, so we here propose an average threshold of Q = -450 HU as a compromise that should yield good results in most cases. We will analyze the sensitivity of our results on the choice of this threshold later. The reconstruction of the nasal geometry from CT scans can be performed completely automatically. To do this, we developed python code that uses the marching cubes algorithm implemented in Slicer 3D (Fedorov et al 2012) to turn the images of the CT scan into a surface mesh. The tip of the nose is then located automatically and used as a reference point to remove surface elements attributed to the ears and the mouth, see figure 4(a). The same reference point is also used to remove parts of the face so that only the external nose remains, see figure 4(b). We then place a box around this external nose, where ambient pressure will be imposed. Additionally, the pharynx is elongated artificially to model part of the larynx and prevent boundary artifacts in the fluid dynamics simulation (Taylor et al 2010). The geometry reconstruction took (80±20) s per mesh on a standard laptop. The automatic addition of the box and elongation of the larynx can fail when non-manifold vertices and edges are created during geometric manipulations. To circumvent such problems, we could randomly vary the size and the position of the boundary regions slightly until the reconstruction completes. However, this was not necessary for the geometries discussed in this paper. Taken together, this procedure automatically reads a series of CT images and creates a surface mesh together with suitable boundary regions, which can be used for the airflow simulation discussed below.
We performed the automatic geometry extraction for three different patients using various thresh-oldsQ. Figure 5(a) shows a direct comparison of the CT image of one patient and the reconstructed geometries for three values ofQ. Clearly, larger thresholds lead to more expanded geometries because they imply that more regions are interpreted as air-filled cavities. The figure indicates that changing the threshold by 100 HU moves the boundaries of the geometry by less than a pixel. To see whether such a change is significant, we first consider the total volume of the geometries as a function ofQ. Figure 5(b) shows that increasing the threshold by 100 HU increases the respective volume by less than 10%. While this might be negligible, the change in the geometry might have a strong effect on the fluid dynamics.

Impact of reconstruction threshold on flow resistance
We simulated the flow in the geometries of the three patients as described in the Methods. Figure 6(a) shows the streamlines of a single simulation. We here do not analyze the simulation in detail, but rather measure the pressure dropΔP between the nostrils and the throat together with the volumetric fluxQ through the cavity, as indicated in figure 4(d). The resistanceK to the flow is then given by = D K P Q. We measured K as a function of the thresholdQ for three examples. Figure 6(b) shows that changing the threshold by 100 HU changes the resistanceK by 40 % in the vicinity of the threshold Θ=−450 HU. This indicates that choosing the correct threshold is important.
The strong dependence of the measured resis-tanceK on the thresholdQ used for the reconstruction can be understood from a simple theoretical model of the flow in the nasal cavity. It has been shown that the narrow central region of the nasal cavity behaves similar to two parallel plates separated by the gap widthℓ (Zwicker et al 2018). Since the flow in the nose is laminar, the associated resistance is given by where h = -· · 1.8 10 Pa s 5 is the dynamic viscosity of air, L is the length of the cavity, and W is its width, so that ℓW is the cross-sectional area. Typical values for the geometric parameters are ℓ≈3 mm, L≈7 cm, and W≈10 cm (Zwicker et al 2018).
Reconstructing the geometry with different thresholdsQ influences the small scales of the geometry much more strongly than the overall shape. Consequently, the small gap widthℓ will be affected significantly, while the widthW and the lengthL will hardly change. The gap width ℓ of a reconstructed geometry is a function of the physical gap widthℓ * and the chosen thresholdQ, similarly to the channels sizes of the CT phantom discussed in section 3.1. Consequently, we can use the average of equation (4) to express ℓ as a linear function of Q, where the factor of 2 stems from the fact that the gap is defined by two walls. Here, c denotes the sensitivity of the reconstruction on the choice of threshold, which we determined in section 3.1.2. The resistanceK changes if the thresholdQ, and thus the gap widthℓ, is varied. Combining equations (6) and (7), we obtain to linear order in Q around a reference threshold * Q with associated resistance * * = Q ( ) K K . Note that this linear dependence of K on Q is observed in the data shown in figure 6(b). For typical valuesℓ * ≈3 mm and c » -10 3 mm HU , changing the threshold by 100 HU thus implies that the resistance changes by about 20%, comparable with what we observe in figure 6(b).
The simple model presented here allows us to extract information about the nasal geometry from the dependence of the resistanceK on the reconstruction thresholdQ. Using the sensitivityc » -10 3 mm HU and * Q = -450 HU, we fit equation (8) to determine * K and * ℓ . The resulting gap widths * ℓ indicated in figure 6(b) are very close to the expected value of 3 mm. Note that ℓ * represents an average gap width, which is related to the hydraulic radius of the nasal cavity. Taken together, this analysis indicates that the influence of the threshold on the flow simulation is adequately described by the linear model given in equation (7). This linear model combined with equation (6) can then be used to predict how changing the threshold affects the resistance without performing costly numerical simulations.

Discussion
We presented a validated protocol to reconstruct geometries of nasal cavities. The segmentation step of the protocol is threshold-based to provide a robust implementation that can be automated. We demonstrated that choosing the right threshold value is crucial since small variations influence the geometry and lead to significant changes in the resulting airflow. Using the known geometry of the CT phantom, we established that the optimal threshold is Q » -450 HU. This value is close to the optimal thresholds determined in other phantom studies (Nakano et al 2013, Zeiberg et al 1996. Moreover, this value and our theoretical analysis verify the heuristic of choosing a threshold halfway between the radiodensity of the soft tissue and air, e.g. employed by Wootton et al (2014). Note that this value is considered to be low by other studies based on the fact that the maxillary sinuses might not be connected to the main nasal cavity in the reconstructed geometry (Croce et al 2006, Quadrio et al 2014. However, we argue that this observation is a poor criterium for selecting reconstruction thresholds since the maxillary sinuses do not affect the flow through the main nasal cavity. In fact, the narrow connection between a sinus and the main nasal cavity could be clogged, which would be captured by a reconstruction that separates the two regions. While the reconstruction of the sinuses is unimportant, we have shown that the reconstruction of the narrow regions of the main nasal cavity is critically affected by the choice of the threshold. Consequently, we verified thresholds based on their capability to reconstruct the narrow channels of a CT phantom of known geometry. Note that this also implies a faithful reconstruction of coarser structures. We also determined how varying the threshold affects the reconstruction of the geometry and the subsequent simulation of the airflow. We established that varying the threshold by 100 HU would move the walls of the nasal cavity by about 0.1 mm when the CT resolution is » h 0.4 mm. This modification of the geometry will change the overall resistance of the nasal cavity by roughly 20%. However, this resistance is a global quantity and local quantities, like the wall shear stress (Inthavong et al 2014), might thus change even more strongly. This large effect of minute changes in geometry on the flow are due to the non-linear dependence of the flow resistance on geometric parameters, see equation (6). This sensitive dependence needs to be considered in all reconstructions of nasal geometry. For instance, typical manual and semimanual segmentations can introduce errors of up to 1 mm (Tingelhoff et al 2007), which implies a more than two-fold change in resistance for a characteristic gap width of ℓ≈3 mm. Even worse, some fully-automatic segmentation methods used in image analysis often lead to even larger deviations (Huang et al 2016) and are thus unsuitable for reconstructing nasal cavities. Moreover, smoothing procedures, which are often employed to prepare the surface for the generation of the computation mesh, can move boundaries by up to one pixel (Gambaruto et al 2009). In the case of a pixel size of h≈0.4 mm and ℓ≈3 mm, this implies that the resistance changes by a factor of 2. Even if smoothing moves the boundary by only half a pixel, we predict the resistance to change by 50%. Consequently, it is important to reconstruct the geometry faithfully and process it as little as possible. We here showed that smoothing is not necessary when the airflow is simulated using the Lattice Boltzmann method. Additionally, this method has the advantage that it employs a regular mesh, whose quality is controlled by a single parameter, the mesh resolution.
The combination of a verified threshold-based segmentation with a simulation based on the Lattice Boltzmann method has the advantage that it is very robust and can thus be automated efficiently. In fact, we already showed in this paper that a completely automated simulation, including segmenting the CT scans, creating the mesh, and simulating the flow, is possible. In the future, it will become important to also automate the analysis step to provide useful visualizations to the provider. However, a better understanding of the relationship between nasal geometry, the airflow, and the symptoms is required for this. Our automated workflow could provide the necessary data for such future investigations.
Further improvements of our algorithm could lead to geometric reconstructions with less variability and thus more reliable flow simulations. For instance, the resolution of the CT scans could be increased by using more slices, a windowed reconstruction focusing only on the nasal cavity, larger image matrix sizes, or new technological developments, like photoncounting CT (Shikhaliev et al 2005). Moreover, advances in the removal of typical CT artifacts, like beam hardening due to bony structures, could further improve the image quality and thus the geometry reconstruction. Our proposed reconstruction algorithm using marching cubes with a fixed threshold could be improved by considering variable thresholds in the future. Here, our results could serve as a baseline for quantitative comparisons. The airflow simulations could be improved by increasing the resolution of the Lattice Boltzmann grid, which should be possible by using more compute power, employing adaptive meshes (Lintermann et al 2013a), or removing negligible parts, like the sinuses. It has also been proposed that the time point of taking the CT scans could be important (Kim et al 2013). This is because the airflow in the two nasal cavities varies with time on the order of hours, a process known as the nasal cycle (Kahana-Zweig et al 2016).
Finally, our segmentation method based on a fixed threshold is simple enough that it might also be applicable to other image modalities, e.g., Magnetic resonance imaging (MRI). In this case, we would also expect that the reconstruction accuracy scales with the image resolution and that thresholds close to the mean intensity between the air-filled cavity and the surrounding tissue are optimal. Other image modalities might reduce the radiation patients are exposed to and prevent scattering artifacts from metal implants, but a quantitative analysis of the reconstruction quality is necessary to assess their performance.

Conclusion
We showed that the simulated airflow in the nose is very sensitive to the details of the reconstruction of the nasal geometry from underlying CT scans. It is thus paramount that geometric variations are kept to a minimum and all details are spelled out in studies on nasal airflow. We here examined a threshold-based reconstruction for which we determined an optimal threshold parameter of Q » -450 HU based on the comparison with a CT phantom. More generally, the threshold should be close to the mean radiodensity of the two regions that should be separated, so that our algorithm will also be applicable to other CT scanner settings and other geometries with fine structures and high contrast. We also showed that the Lattice Boltzmann method is suitable for simulating nasal airflow and can be used in a completely automatic workflow. Taken together, our work thus provides an important step towards the automatic analysis of nasal airflow based on CT scans alone and it establishes a quantitative baseline for future improvements.