Systematic Method for Morphological Reconstruction of the Semicircular Canals Using a Fully Automatic Skeletonization Process

: We present a novel method to characterize the morphology of semicircular canals of the inner ear. Previous experimental works have a common nexus, the human-operator subjectivity. Although these methods are mostly automatic, they rely on a human decision to determine some particular anatomical positions. We implement a systematic analysis where there is no human subjectivity. Our approach is based on a specific magnetic resonance study done in a group of 20 volunteers. From the raw data, the proposed method defines the centerline of all three semicircular canals through a skeletonization process and computes the angle of the functional pair and other geometrical parameters. This approach allows us to assess the inter-operator effect on other methods. From our results, we conclude that, although an average geometry can be defined, the inner ear anatomy cannot be reduced to a single geometry as seen in previous experimental works. We observed a relevant variability of the geometrical parameters in our cohort of volunteers that hinders this usual simplification.


Introduction
The spatial orientation of the three semicircular canals (SCs) has been a matter of interest since early anatomical studies of the inner ear and temporal bone. Recent clinical studies show that SCC spatial orientation plays a main role in specific pathologies (i.e., positional vertigoes [1]). In addition, the specific morphology of the vestibular system can explain the variability on the induced adverse symptoms (i.e., vertigoes and dizziness) on patients that are subjected to MRI high magnetic fields [2][3][4]. Hyrtl first, and then Sato, Siebenmann and Schoenemann later, in the transition from the 19th to the 20th centuries, gave the first measurements of the angles between the 3 SCs. The methodology was mainly based on corrosion casts and major differences were seen between the results given by different authors [5,6]. Since then, other methods have been developed to study the anatomy of the SC post-mortem.
With the advent of modern imaging technology and in particular with the ability to assess axial scans and perform reconstructions in several planes, more detailed in vivo evaluations in higher amounts of bones have been performed. However, even for computed-tomography (CT) and magnetic resonance imaging (MRI) studies, the differences in results are still significant, although of lower magnitude than before. CT scans [7,8] show some variability in results except for the angle between the horizontal SC and posterior SC, where results are more consistent. When compared to results in more recent anatomical measurements made in explanted temporal bones [9], the differences are smaller, and can, in general, be attributed to the exclusion of some parts of the canal in the analysis (common crus and ampullae). In the case of MRI, data also differ and the angles between canals (in the same ear) are close to 90 • in the case of the anterior and posterior SC [10] as well as in the case of the anterior and horizontal SC [11].
In general, the canals do not lie in specific planes in the space. Since the first measurements, it became clear that the canals follow a sinuous course and only their mid-portion is appropriately in a plane; that portion makes an obtuse angle with the ampullated and non-ampullated ends. In addition to this, other attributes of the canal (non-planarity, radius of curvature and length of the common crus), well documented with micro-CT scans [12], render the spatial orientation characterization task complex and elusive.
The relevance of spatial orientation to normal range of function dates back to the very first experiments in vestibular physiology. Since Crum-Brown, coplanarity (as documented and experimented in the horizontal canals) was extended to the six SCs. The planes for the left anterior canal with the right posterior (LARP) and for the right anterior with the left posterior (RALP) considered, as anatomical and functional units, became a matter of anatomical interest throughout the years. More recently, since their specific evaluation has become a standard in vestibular testing at bedside [13], the interest has grown.
Previous experimental works based on post-mortem dissections [9,14] or imaging [7,10,15] have a common nexus, the data processing protocol and extracted information always depend on subjective factors that are involved on how to delimit the section of the canal that must be considered in the analysis. In all the previous works, this is the decision of an expert operator. Recently published works by A.P. Bradshaw et al. proposed a mathematical method that allows 3D measuring and modeling of the semicircular ducts in living humans [12]. Their mathematical method automatically reconstructs the semicircular ducts geometry from computed tomography images validating the results with micro-computed tomography as ground truth. The final result is a centroid path of the cavities that allows extracting global characteristics of the vestibular cavity. Although the method proposed by A.P. Bradshaw et al. largely suppresses the arbitrariness, the initial step of the protocol still requires an operator to select the initial point of the ducts. In addition, this method accurately defines the centerlines of the ducts but ignores the SCs local diameter as a geometric parameter that can help to improve the characterization of the other parameters.
In this study, we were interested in defining in a group of normal subjects the planar orientation of the three SCs with the use of 3D high-resolution MRI images. We present an automated method that allows computing the SC duct morphological characteristics without subjectivity or operator variability. Our approach is based on an algorithm that computes the geometric characteristics of the ducts (3D orientation, diameter, and non-planarity) through a skeletonization process based on a fictitious repulsive force. In addition, the algorithm tests the influence of retaining only a fraction of the ducts (mimicking the potential inter-operator variability).

Subjects
The cohort of participants on this study are healthy volunteers recruited among the students and staff of the University of Navarra. The study was approved by the Ethics Research Committee of the University of Navarra. The inclusion criteria were healthy volunteers aged 18 years or older. The exclusion criteria included any potential illness related to the inner ear, such as dizziness, a medical history of hearing or visual problems and MRI incompatibilities such as claustrophobia, pregnancy, metal implants, etc.
Before inclusion in the study, subjects were evaluated at the Department of Othorinolaringology of the Clínica Universidad de Navarra by an expert otholaringologist in order to verify the absence of pathologies. The normality tests included normal Video Head Impulse Test (vHIT), Cervical Vestibular Evoked Myogenic Potential (cVEMP), Ocular Evoked Myogenic Potential (oVEMP) and Romberg test. After these evaluations, we retained a final group of 20 participants from an initial set of 30 volunteers. The volunteers received oral and written explanations and provided written informed consent.
The final volunteer group after the selection process was formed by 20 subjects (age span, 19-45 years, with an average age of 23.9 years). The group was divided into 9 men with an average age of 26.5 (±7.8) years and 11 women with an average age of 21.72 (±1.9) years.

Imaging Protocol
MRI studies were performed on a 3T whole-body scanner (Trio TIM, Siemens AG, Erlangen, Germany) using a 32-channel head array. First, a T2 weighted 2D fast spin-echo (FSE) sequence was employed to acquire anatomical images in coronal orientation, which were subsequently used to localize the high-resolution scan.
High-resolution images were acquired using a three-dimensional Fourier transformation constructive interference in steady state (3DFT-CISS) gradient-echo sequence [16] in axial oblique orientation with the imaging parameters shown in Table 1. The slab was positioned to include the inner ear, based on the FSE localizer ( Figure 1).

Labyrinth Volume Reconstruction
We restricted the analysis to a sphere of diameter 20 mm centered around each one of the semicircular ducts. From the 3DFT-CISS sequence, we recovered a matrix with 320 × 320 × 72 voxels that correspond with the x-, y-, and z-axes (usual Reference Coordinates System (RCS) reference frame) and a resolution of 0.4 mm in the head-foot direction and 0.6 mm in the two other directions. Our first step was to uniformize the resolution to 0.4 × 0.4 × 0.4 mm 3 using a bicubic interpolation in the x-y plane (transverse plane).
The MRI sequence used in this experiment provided high signal intensity in all spaces filled with intralabyrinthine fluid, while the signal outside was very small (5%) compared to the value inside the ducts. We could use this factor to split these volumes through a threshold value that binarized the original image. A typical problem of the binarization is the voxelization caused by the loss of information when a threshold is applied to a continuous variable. To prevent this, a common approach is to upscale the image and then binarize. We applied an upscale by a factor 3, thus the final binarized image used for the volume rendering was 0.13 × 0.13 × 0.13 mm 3 . In Figure 2, we present an example of this process, the original data, and two binarized sections, one using the original resolution and the other one with the upscaled version, can be observed.
In the last panel, we present an isosurface that corresponds to 20% of the maximum intensity value in the interpolated image. This value is used as a threshold for the binarization in the whole group of volunteers. The specific value used as threshold is not critical, as far as it remains far from the two extrema (the values inside and outside the canals).

Labyrinth Skeletonization Method
One of the possibilities for the analysis of the ducts' orientation is to assume that they follow a specific geometry, such as a circle, ellipse, or other 3D shape and then adjust the one that best describes the data. Here, we relaxed this restriction and obtained an arbitrary shape of the vestibular system, and then analyzed how the characteristics of these reconstructed ducts depend on the specific details of the analysis.
The shape of the ducts was obtained through a skeletonization [17] process of the binarized volume. This name encompasses a broad battery of different algorithms that generate a curve retaining the topology of the original 3D volume ( Figure 3). The difference among the various algorithms depends on how this curve is defined. Skeletonization methods have been studied for decades and their results have been tested by many authors in different experimental and applied topics, such as virtual navigation, computer graphics animation, medical applications, data analysis, etc. [18][19][20][21][22]. Some methods are based on the definition of this line as the set of points that are equidistant to the surface (using a geometric distance that can be Euclidean or not). Other methods erode the original volume reducing the thickness until it collapses to a line (thinning methods). However, these two families are very sensitive to noise and to the roughness of the volume boundary.
Here, we focused on a specific family called Potential Field Methods that are less sensitive to these effects. The idea is that a fictitious repulsive and conservative force is associated to each one of the surface elements, and then a scalar potential field is created in the interior of the considered volume. This potential field will decrease when we move away from the surface and when we approach the central part of the volume. For any possible cross section of each canal, there will be at least one point inside where this potential reaches its minimum value. The skeleton defined through this algorithm is the line d that connects the set of minima that will be recovered from any plane that chops a duct.
From a practical point of view, we can define a vector force field from this potential. Any fictitious particle under the action of this force will be pushed to the skeleton and once it reaches this region it will restrict its dynamics to displacements along this line.
We used an implemented MatLab (The MathWorks, Massachussets, USA) toolbox named Volume Skeleton Toolbox (VSK) [23]. In this toolbox, the specifically used method is called Generalized Potential Field method. This script splits the volume into boundary elements (voxels with at least one empty neighbor voxel) and interior voxels. Then, a potential field V is generated (we used a spatial dependence with the fifth power V ∼ 1/r 5 of the distance to the surface r ). A dynamical system is defined where the acceleration of the fictitious particles depends on Φ = −∇V. These particles will find (along the trajectory obtained connecting all these minima) some points where the dynamic stops (critical points). We can classify these critical points as local minima of the potential and saddle-points using the Jacobian of this function Φ. The skeleton can be recovered simulating the evolution of points close to the saddle-points and that propagate towards the local minima, with a dynamics always restricted to the skeleton itself (because in any other direction different from the one defined by the skeleton the potential will increase, and thus the particles will experience an acceleration pushing them back to the skeleton). These critical points are obtained through a Newton iterative method with a convergence criterion of 4 · 10 −3 (being 1 the voxel size). The step for trajectory tracing is 0.25 and the maximum number of steps is 200. A more detailed description can be obtained in [18,20].
Although VSK has been tested in numerous geometries, we noted the presence of several artifacts (extra lines) in our working volumes with the parameters specified above due to the ducts complexity and the cochlea presence in the working volume. We reunified redundant points when their mutual distance was smaller than one voxel. In addition, we deleted the cochlea points because they fall out of the scope of this analysis. With these new steps, the merging of all the trajectories generated on the semicircular ducts define the centerline (Figure 3a). Table 2 shows a schematic description of the centerlines definition code. Perform a binarization (Intensity-based, threshold 20%) ; 5 Perform a boundary/internal voxels labeling (Neighbourhood based); 6 Define saddle-points/local minima (Potential field); 7 Skeleton curves generation (Newton iterative method); 8 Reunify redundant points/curves; 9 Output: centerline curves; 10 end

Geometrical Parameters and Morphology Characterization
It has been clearly stated above that an accurate definition of the SCs morphology is important. One crucial point is that the cupulae are part of the vestibular cavity producing an increase in the diameter of the SCs and should be withdrawn from the calculations. Here, we explored a systematic way to restrict the geometrical analysis to the ducts themselves, reducing progressively the length of the canal considered and testing the effects on the geometrical parameters.
The procedure is as follows: once the original central lines of the ducts have been identified (Figure 3a), we characterized the minimum set of geometrical parameters that can describe their morphology. From these recovered centerlines, our first step was to obtain a plane that fits best (least-squares of the normal distance to the plane) of the set of points that define the trajectory. From the normal vector of this plane, we recovered the 3D spatial orientation of each canal, and we could obtain the relative orientation between each pair (the angle between these vectors). Other parameters that we obtained are the non-planarity (normal distance from the canal to the plane) and the ellipse that fits best (least-squares algebraic distance) to the projection of the canal on the plane.
Depending on the cupula section considered as part of the semicircular ducts, the planes that represent the canals show slight differences. For each point of the original trajectory, we computed the distance from the surface of the volume r( d). In Figure 4, we present these values for all the canals and volunteers. As the different canals have slightly different sizes, we normalized the position along each duct with its corresponding length. The variability of the canal radii between different volunteers (gray lines) and along the canal position is larger than the uncertainty due to the algorithm (i.e., the size of the binarized voxel).  These distances were retained as an indicator of the diameter of the duct along the centerline and we refer to them as radius, as each section can be assimilated to a cylinder. The cross-section of the canal defined by a plane perpendicular to the skeleton has an arbitrary shape, in most cases a convex shape different from a circle. Here, we used the generalization of the definition of diameter as the largest distance between two tangent straight lines at opposite sides of the shape boundary, and the radius is the shortest distance from the centerline to the surface. The local diameter T of the duct varies along the trajectory d and corresponds to T( d) 2 r( d). We defined a threshold for these radii to determine the ampullae region: all points above that threshold were assumed to be part of the cupula of the channel and were removed from the calculation. We repeated the geometrical analysis described above determining the planes, ellipses, non-planarities and angles. Once this step was finished, we slightly decreased the threshold and checked its effect on the new parameters. The main advantage of this process is that we could sequentially determine the effects of the excluded area in the SCs definition without subjectivity due to the action of the operator. Table 3 shows a schematic description of the SCC morphological characterization code. Define canals radii (Distance from volume surface); 4 Define ampullae threshold range ini and end (Local radii thresholding); 5 for (Ampullae ini ) to (Ampullae end ) 6 Remove ampullae points from centerlines; 7 Define plane fits best centerlines (Least-squares of normal distance to the plane); 8 Define normal vectors to planes; 9 Define relative orientation between canals (Angles between vectors); 10 Define non-planarity (Normal distance from canal to plane); 11 Define representative ellipse (Least-squares algebraic distance); 12 end 13 Output: Radii, Planes, Ellipses, Nonplanarities and Angles 14 end

Results and Discussion
Based on the proposed method, we computed the centerline for all three SCs of each volunteer's right ear. From the initial skeleton ( Figure 3, top), we removed the common crus structure. The initial SCs centerlines considered in the study are the remaining sections of the ducts: the anterior and posterior ducts paths are defined from the common intersection point to their ampullae. Figure 4 shows the local radius for all the studied volunteers and the mean value for each SCs. As expected from the data available in the literature [24], all the A, P and H semicircular ducts present a central region where the radius reaches a roughly constant value r ∼ 0.5 mm regardless of the volunteer.
A histogram of the whole data gives an accurate approximation of the SCs' radii ( Figure 5). There are no data for radius smaller than 0.3 mm, being the mode within 0.5-0.6 mm for all ducts, which corresponds to the value obtained on the central regions ( Figure 4). In contrast, the extreme values of the plots show an increase of the radius value that can be up to three times (1.5 mm) the central region radius. These last values correspond to the ampullae. As happens in other studies, in many volunteers, the SCs are far from simple toroids and are bent ( Figure 6). This bending is more relevant close to the ampulla region, and consequently the orientation of the plane that fits to the trajectory that defines the duct (Figure 7) is largely affected by the presence of the ampulla, be it complete or partial, in the calculations. In all previous works, the determination of the ampullae position was done by an expert operator. There is an inherent inter-operator variability due to this decision, and we sought to evaluate its role on the determination of the geometric parameters. We explored the effect of progressively removing the ampullae based on a threshold of the radius that defines the section to be removed. Once the distribution of radii r( d) in decreasing order was obtained for each canal, the threshold was defined using an increasing percentile. A zero value means that the whole trajectory was retained, and 100% that everything was rejected. For any intermediate value, we found the first and last points that fulfill the requirement. Due to the roughness of the channels, the diameter is not a monotonic function of the canal position d, and it may happen that the diameter increases again in some sections (see Posterior SC behavior in Figure 4, for a position d/d Max = 0.8). In those situations and for a given threshold, we could retain sections with a diameter larger than the rejection threshold value. This is why (see Figure 8, bottom right) we obtained retained lengths for the posterior duct that can be up to 10% larger than the value that corresponds to the rejection threshold.   The results show that in some cases there is a clear dependence of the computed SCs pair angle on the rejection threshold. PH angle presents the most regular behavior, with a plateau close to 90 • regardless of the discarded region. Nevertheless, there are two volunteers whose angles are around 80 • .
However, for the AP and AH angles, we have a large variability. The AH angle depends on the percentile threshold applied, although, for a certain value of canal discarded, the angle tends to an asymptotic angle around 90 • too. On the other hand, the AP is totally dependent on the duct retained, and the angle has a growing trend without saturation for the range of explored lengths. Even more important, on the AP angle and for small percentile thresholds, there is a region where the angle changes more than 5 • , and it becomes evident that we must expect a large inter-observer variability for this value when an operator has to decide the position of the ampullae. A similar result can appear on the AH angle. We do not expect a large inter-observer variability on the average PH angle when this value is computed for a large set of volunteers. This could explain the high variability on the results present in the literature with different techniques (skull sectioning, imaging techniques and mathematical model).
Concerning the non-planarity of the ducts, we present in Figure 9 the normal distance to the fitting plane for the trajectories that define the centerline of the SCs. We present two situations: the whole duct and a rejection threshold using the 40th percentile. In both cases, the distance is small, comparable to the diameter of the SCs, except when we are in the ampulla regions. This result is even more evident in the PSC, where the ampulla induces a deviation of 1.5 mm. When the rejection threshold is applied removing the ampullae, these distances decrease to values smaller than the radius of the canal (0.5 mm), and we can conclude from this analysis that the central regions (that retain more than 60% of the whole length) of the SC lie in the fitting plane. Figure 9. A, P and H non-planarity (defined as normal distance) referred to the reference plane: (left) this distance when the whole duct is retained; and (right) the evolution of this distance when the 40th percentile is discarded. Grey lines represent different volunteers, the black line is the average and the red lines represent the limits for one standard deviation for each SCs.
An important statistical indicator about the influence of the ampulla presence in the non-planarity is the standard deviation along the trajectory of the distances between the SCs ducts and the reference planes. Here, we present the average value of this statistical computed for all the volunteers ( Figure 10). From this analysis, it is clear that for the PSC and HSC a rejection threshold larger than 10% does not produce a noticeable effect, as the global distance does not substantially change. However, on the ASC, this non-planarity remains important even using a rejection threshold of 40%. This result indicates that the Anterior SC is the most curved canal, and also explains the behavior of the AP and AH angles, which strongly depend on the applied threshold. Even though the separations we obtained for a rejection threshold of 40% are s H = 0.08 mm, s A = 0.15 mm, and s P = 0.10 mm, where · denotes the average computed over the whole group of volunteers, individual separations can be as large as s H = 0.36 mm, s A = 0.34 mm and s P = 0.36 mm for the whole canal and s H = 0.17 mm, s A = 0.24 mm and s P = 0.17 mm for 40% rejection, respectively. Bradshaw et al. [12] described the bending of the SCs ducts assuming that they lie on a curved surface (a generalized cylinder) instead of a plane. From the data published in the Appendix of that paper, we could estimate that the non-planarity for their model assigns a standard deviation of the distance to the respective planes of s A = 0.29 mm, s P = 0.19 mm and s H = 0.16 mm, values which are compatible with the values we obtained using a rejection of roughly 5% (see the thick black lines on Figure 10). From our data for posterior and horizontal SCs, we observe that a slight variation on the decision about where are the ampullae located could modify the parameters of this model. Again, the decision on how the ampullae is considered can affect the parameter values of any model.
As in previous works, we verified that the centerlines are very close to ellipsoidal trajectories, thus we extracted the parameters that can be obtained from the individual paths. A least square fitting was done (based on an algebraic distance) for each volunteer and canal. Figure 11 shows in red the ellipse whose semi-axes are the average of the whole set of ellipses that better describe each of the SCs projected on its corresponding plane. We present two situations, retaining all the duct or with a rejection threshold of 40%. The individual trajectories appear in gray in the same figure. It must be noted that the ampullae position, in most of the cases, is located at one end of the major axis. Figure 11. Centerlines projected on the respective reference planes for the A, P and H ducts: (top) the behavior for the whole ducts; and (bottom) when 40% is discarded. The red thick line correspond to the ellipse that best fits all the data from the different volunteers.
We can appreciate small differences of the fitting parameters of the ellipses in both 0% and 40% rejection thresholds. In Table 4, we present the variability of the semi-axes. In both cases, our data include the average and their corresponding standard deviation for the whole set of volunteers. We include for reference the values published in two previous works, by Bradshaw et al. [12] and Daocai et al. [24]. The values from this last study were inferred from their published values, as they only included the size of the circular cavity defined inside each canal. As our characterization is done using the central path along the ducts, we enlarged their values with the canal radii values of our work. We can appreciate that there is no real difference between these three studies, as the variability (defined from the standard deviation in our cohort) of these parameters is a 10%. However, it should be noted how the rejection of the ampullae region can affect in a similar degree those parameters and the operator action that removes that part is a source of artifact that increases the data variability. Table 4. Mean value and standard deviation of the fitting ellipse control parameters (major and minor ellipse axis) for the three different SCs. The values presented correspond to two cases: retaining the whole duct and using as rejection threshold the percentile 40%. We include also the results of Bradshaw et al. [12] and Daocai et al. [24] for comparison.

Anterior
Posterior

Conclusions
A new method to extract the geometrical parameters of the semicircular canals of the inner ear is presented. Our approach relies on the skeletonization of a binary image obtained from a high-resolution magnetic resonance image acquired from a cohort of 20 volunteers. This process avoids subjective decisions made by the observer, and systematically analyzes the quantitative effect of the manual definition of the canals.
Other previous experimental works have analyzed this problem using similar techniques (i.e., X-ray TC, MRI, or skull analysis) but all of them need the subjective input of a human operator at some point during the analysis. This introduces an inter-operator variability that our method completely removes. Even more interesting, we could mimic the role of this operator by systematically reducing the canal length taken into account on the calculations and observe how this affects the retrieved geometrical parameters.
The whole set of geometrical parameters can be obtained: width, radius, angle of the functional pairs, and non-planarity. Our results reveal that there is a large variability across the parameters, which makes it difficult to restrict all the geometries to a single one that represents the whole group.
Author Contributions: All authors contributed to the conceptualization, methodology, software, validation, investigation and writing; and I.C.-D. and J.B. contributed to the formal analysis and visualization.

Funding:
The authors would like to thank the support of the Spanish Government through grants FIS2011-24642, FIS2014-54101-P and FIS2017-83401-P and the Caja Navarra Foundation (2014 call). One of the authors (I.C.) acknowledges Asociación de Amigos de la Universidad de Navarra for a research grant.