# Evaluation Algorithm of Root Canal Shape Based on Steklov Spectrum Analysis

Authors: Dongqing Wu; Jian Gao; Xiaoli Hu; Zhengtao Xiao; Zhuwei Huang; Lanyu Zhang; Xin Chen; Yunbo He
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

In recent years, we have seen more and more interest in the field of medical images and shape comparison motivated by the latest advances in microcomputed tomography (μCT) acquisition, modelling, and visualization technologies. Usually, biologists need to evaluate the effect of different root canal preparation systems. Current root canal preparation evaluation methods are based on the volume difference, area difference, and transportation of two root canals before and after treatment. The purpose of root canal preparation is to minimize the volume difference and ensure the complete removal of the smear layer. Previous methods can reflect some general geometric differences, but they are not enough to evaluate the quality of root canal shape. To solve this problem, we proposed a novel root canal evaluation method based on spectrum and eigenfunctions of Steklov operators, which can be served as a better alternative to current methods in root canal preparation evaluation. Firstly, the ideal root canal model was simulated according to the root canal model before and after preparation. Secondly, the Steklov spectrum of the two models was calculated. Thirdly, based on the spectrum and the histogram of the Gaussian curvature on the surface, the weight of each eigenvalue was computed. Therefore, the Steklov spectrum distance (SSD), which measures shape difference between the root canals, was defined. Finally, the calculation method that quantifies the root canal preparation effect of root canals was obtained. Through experiments, our method manifested high robustness and accuracy compared with existing state-of-the-art approaches. It also demonstrates the significance of our algorithm’s advantages on a variety of challenging root canals through result comparison with counterpart methods.


---

## Body

## 1. Introduction

In recent years, imaging technology represented by 3D digital geometry technology has widely been used in the medical field. The 3D model can display critical information such as the three-dimensional topological structure, geometric information, and anatomical structure of the lesion more intuitively. It lays a solid foundation for successive operations, such as simulating the operation, human-computer interaction, quantitative morphometry, and telemedicine [1]. Researchers have paid more and more attention to the issues of segmentation, registration, difference evaluation, and classification retrieval of 3D medical models. With the popularity of measurement tools represented by μCT, the evaluation of anatomical structure before and after root canal preparation has become a hot topic in medical imaging [2–9].The root canal is the part of the cavity in the middle of the tooth. The upper part of the cavity is broad, called the pulp cavity, and the lower part has a tubular structure from which the blood vessels of the dental nerve and the nutrient nerve are derived, as shown in Figure1. In [10], Vertucci investigated 2,400 permanent teeth and categorized them into eight types. He studied the morphology and root canal of adult and established the evaluation criteria of treatment difficulty. From the point of view of dental anatomy, root canal cleansing and complete obturation with an inert filling material are very necessary when the root canal morphology of teeth is affected. Root canal therapy generally consists of three stages: root canal preparation, root canal disinfection, and root canal obturation. Root canal preparation is the critical step for root canal obturation therapy, including root canal cleaning and formation. Before root canal preparation, the location, length, shape, and number of root canal orifices need to be determined. Under the guidance of the dental operating microscope, ultrasonic root canal files and trumpet hand files are used to probe the root canal to the apical orifice. Root canals with small or blocked root canals should be cleaned and dredged with ethylenediaminetetraacetic acid (EDTA) or other solutions. Infected dentin should be removed and sterilized by root canal instruments to facilitate root canal obturation. The prepared root canals should contain the original root canal shape, retain the narrow part of the apex, and form the taper from the coronal to apical [3].Figure 1
Root canal and dentin in the first premolar. (a) Root canal and dentin. (b) Left view. (c) Right view. (d) Four Dicoms for top to bottom.
(a)
(b)
(c)
(d)To improve the effect of root canal therapy, various root canal preparation systems have been developed to maintain the original morphology of root canal anatomy (including using different instruments and cooperating with different solutions to eliminate the smear layer of the root canal wall). The results before and after the operation among groups were tested by the significance test, such as Mann–WhitneyU test [3, 6, 7, 9]. With the development of medical imaging measurement technology such as CT, more and higher resolution image data can be collected, which boost more accurate instrument evaluation. μCT is a nondestructive 3D imaging technology, which can clearly understand the internal microscopic structure of the sample. Hammad et al. [11] compared various instruments and considered that the root canal measurement results of the μCT are consistent with the actual results. Kierklo et al. [6] evaluated the root canal preparation operation based on 2D data; that is, the cross-sectional images collected by the μCT were compared before and after preparation. In [12], 3D models were reconstructed from images, and the volume and area of the model were used to evaluate the effect of preparation among different systems.In the field of medical imaging, many studies of therapy evaluation are based on the idea of registration before comparison. As far as root canal models are concerned, there were many changes in the models before and after treatment, so the rigid registration method is unsuitable. The nonrigid registration method mainly includes two approaches based on geometric features or image statistics. Both of the two ideas need to extract features or statistical information from the model for registration. However, in most cases, it is difficult to obtain features accurately when the feature extraction is carried out between different models, or when the feature points are not distinct, and always accompanied by outliers. The outliners probably present immense barriers to the subsequent registration algorithm. Moreover, the geometric characteristics and the topology of the root canal model may change before and after preparation. It is also challenging to apply the idea of nonrigid registration to root canal preparation assessment.Ørstavik [3] summarized the evaluation methods of root canal shaping ability of various instruments, including X-ray parallel projection, Bramante model, cross-sectional section method, CT, and μCT methods. The ideas of these methods are to obtain the volume, surface area, and transportation of root canal axis before and after preparation by integrating the slice images [13, 14]. Specific methods are as follows:(1)
Divide teeth with highh into N equal parts at the interval of δ=h/N, then get images Ii1N(2)
Measure the cavity areaSi and the perimeter of the cavity of the root canal li in Ii(3)
Calculate the total surface area withS=∑i=1N−1δli(4)
Calculate volume withV=∑i=1NδSi(5)
The deviation of the root canal center and the average distance of each cross-sectional cavity are calculated, respectively, after registration, as shown in Figure2Figure 2
Calculating volume and area from slice images. (a) Before canal preparation. (b) After canal preparation.
(a)
(b)Existing methods of root canal preparation evaluation have obvious drawbacks. First, the differences in permeability and taper of different root canals were submerged behind a single geometric quantity, which cannot distinguish the root canals with similar volume and significant shape difference (shown in Figure3). Second, the overall geometry and complex topological structure of root canals are not taken into account. Third, the cavity morphology before and after treatment is often not considered. A large extent of transformation (irregular to ellipsoid-like) occurs, which may cause intolerable errors in both rigid and nonrigid registration algorithms.Figure 3
Five shapes with the same volume (intuitively, if shapes (a) to (d) were sorted from small to large in order of shape differences from the shape (e), the result is (c), (b), (a), and (d)).
(a)
(b)
(c)
(d)
(e)Every root canal has its own individual shape. Accurate knowledge of the canal anatomy is a prerequisite for effective endodontic treatment [15, 16]. At present, the evaluation of the various instruments in root canal therapy is generally based on reexamining the volume percentage of the void (POV) of the root canal. This method typically requires a more extended period and cost. Therefore, by measuring the model of root canal preparation, defining the “distance” between the shapes and calculating the distance between the current root canal and the idealized root canal according to the criteria conducive to the subsequent root canal obturation, it will be helpful to improve the efficiency of this work.For the first time in [17], the method of studying geometric shape based on the spectrum was proposed. The paper answered the question of how much about the shape can be inferred from the knowledge of its geometry spectrum. By calculating the shape spectrum λi and its eigenfunction ϕi as shown in Figure 4, the shape surface can be uniquely determined and vice versa in ordinary cases.Figure 4
Root canal and its eigenfunctions from the Laplace–Beltrami spectrum.In recent years, with the development of 4G/5G technology and CT technology, as well as the unprecedented expansion of the accuracy and breadth of medical 3D data acquisition, more and more researchers pay attention to the massive data of medical geometry. On the contrary, the shape analysis paradigm, theory, and practice based on the geometric spectrum have gradually become a research hotspot [18–21]. Based on the works, the author proposed a novel method based on the geometric spectrum for root canal preparation treatment.The spectrum of differential operator contains rich geometric and topological information of the shape, including curvature, area, and genus. The LBO spectrum and Steklov spectrum have the characteristics of independent of the coordinate system and scale. Although in rare high-dimensional manifolds, there may be different shapes sharing the same spectrum. However, in low-dimensional manifold cases (such as root canals), the geometry of the shape can be uniquely determined by its spectrum. The traditional shape spectrum is based on the intrinsic description of the Gauss curvature, which has an insufficient perception of concavity, convexity, and volume of shape, while the application of shape-based extrinsic geometric spectrum, such as Steklov spectrum, is increasing. It has gradually become a research hotspot [22–26].In this work, we consider the problem of the root canal shape evaluation and investigate its relevance in a selection of literature from medical image and graphics. We propose a shape-aware evaluation method which is a precise evaluator for comparing different root canal preparation systems. Our key contributions can be summarized as follows:Firstly, despite being highly nonlinear and hard to compute, we indicate for the first time that the mapping between a 3D model of the root canal and its Steklov spectrum is addressable with modern mathematical toolsSecondly, for the first time, an evaluation method of root canal shape based on an ideal root canal model is proposed, and an algorithm for generating an ideal root canal model is givenAt last, we evaluate the root canal preparation effect by defining the reasonable distance depending on the histogram of the Gaussian curvature and the Steklov spectrum of root canal models


## 2. Related Research

In recent years, more and more attention has been paid to the evaluation method of therapy based on 3D model reconstruction from data collected byμCT. The high-resolution characteristics of μCT make it widely used in small-scale medical in vivo research. In [27], synchrotron radiation microcomputed tomography (SRμCT) was evaluated concerning the 3D morphology of the rat spinal cord microvasculature. In [28], Wang et al. used parallel programming techniques to accelerate existing image reconstruction algorithms for 3D optoacoustic tomography. In [29], Kobayashi-Velasco et al. compared cone-beam computed tomography (CBCT) with periapical radiograph (PR) in the diagnosis of alveolar and root fractures.Also, more and more applications ofμCT in the evaluation of root canal therapy have been made. The purpose of root canal preparation is to clear the infected root canal anatomical space. Under ideal conditions, on the premise of keeping the centerline unchanged, the inner wall of the root canal should all be processed, and a better taper can be obtained to facilitate subsequent obturation. According to the different geometry of root canal preparation, the evaluation methods of root canal preparation can be divided into two categories.The first category focuses on volume and area changes before and after root canal preparation. In [10], the root canal was segmented into the coronal third, middle third, and apical third. Scanning electron microscopy (SEM) was used to observe the removal of the smear layer; then the images before and after root canal preparation were calculated and compared, which was used to evaluate the therapy. In [4], the geometry of the root canals was described by three different mathematical models: the elliptical model, the 1-ring model, and the 3-ring model, and then a method for analyzing the geometric characteristics of root canals was realized. In [30], Berutti et al. compared the canal curvature and axis transportation after instrumentation with three different systems. In [31], the critical section images of the coronal third, middle third, and apical third before and after root canal preparation were analyzed and compared by several nickel-titanium systems. Anbu et al. [32] observed the efficacy of various technique to fill root canals using spiral computed tomography (SCT). Orhan et al. [2] used both μCT and nano-CT to measure and analyze the volume ratio of filling voids at coronal third, middle third, and apical third and recommend that a voxel size of 11.2 μm can be used as a reliable cutoff value in μCT imaging.The second category of the evaluation method focuses on the transportation of spatial centerline and the balance of periodontal dentin thickness. In [9, 16], root canal deviation and central positioning (or transportation) ability were considered to be essential parameters for evaluating the quality of canal preparation. In [33, 34], the thickness of the mesial and distal root canal wall was scanned and measured, respectively, and then the displacement and axial centrality of root canals in each group were calculated by using Gambill’s formula. In in vivo cases, the coordinate systems of canal measured before and after the operation are not consistent in general. How to register the two root canals remains a problem to be solved. Because of the particularity of teeth, it is not feasible to use registration points in industrial scenes for registration.There are two shortcomings in the above two methods: it significantly relies on knowledge and experience manipulators on understanding the geometry of the root canal. The accuracy and automation of the process are not high; the other is that the evaluation based on local sampling lacks a global shape assessment method suitable for subsequent obturation, so it is not conducive to accurately evaluate the quality of root canal preparation of various instruments.3D shape comparison is a hotspot in the field of geometric analysis. In the field of manufacturing, it is often necessary to evaluate the accuracy of parts processing by shape comparison, such as blades and special-shaped bolts. In the medical field, the manufacture of porcelain fused to metal teeth and artificial bones also need to evaluate the difference of shape. In online massive shape retrieval scenarios, it is necessary to quickly match one or some kinds of shape according to certain retrieval conditions. In the interactive game industry, it is required to cluster and understand facial expressions. In [35], the mapping between two surfaces is expressed by the linear function map, and the difference between different shapes is evaluated by defining reference shapes. Fish et al. [36] introduced a meta-representation that represents the essence of a family of shapes. Handels and Hacker [37] extended the volume-based representation of organs to enable the 3D visualization of organs’ shape variability in the atlas. In [38], a functional map between two shapes based on the multiscale Laplace–Beltrami operator (LBO) shape descriptor was proposed, and the difference between shapes was evaluated by the basis function family of the spectrum to avoid the NP-hard problem of the registration algorithm. In [39], the comparative shapes were embedded in canonical domains and evaluated.Shape comparison methods can be categorized into three types: first, the corresponding relationship between shapes is established, and then the geometric comparison is made; the other is to compress the geometric and topological information of three-dimensional image into a scalar, i.e., shape similarity, by using some global or local shape descriptors; and the third is to evaluate the shape difference by defining a specific norm from the global perspective, such as [40]. These shape descriptors are all shape descriptors based on intrinsic geometry or geodesic equidistant transformation invariant. However, the intrinsic description of shape, such as the LBO spectrum, is unable to distinguish concave and convex shapes (Figure 5). The conicity, concavity, and convexity of root canal shape and the overall permeability are the key factors to evaluate the effect of treatment. Therefore, the extrinsic shape descriptor, which can capture volume difference, is suitable in the scenario. According to [41, 42], the spectrum of the Steklov operator mapping the Dirichlet boundary to Neumann boundary is precisely extrinsic.Figure 5
Inward and outward bumps on the cylinders.The problem of calculating the spectrum of the Steklov operator of shape is generally transformed to solving the boundary value problem of the partial differential equation (PDE). Because of the complexity and discontinuity of the boundary, it is almost impossible to obtain the analytic solution. Therefore, finding the weak solution of the equation is the only feasible choice. Solving the PDE on the discrete point cloud model by the spectral method is an effective way to improve the accuracy and convergence speed. The basic idea is to map the original PDE to an integral equation. By using the orthogonal and recursive formulas of the basis function, the weak solution problem of PDE is simplified to an integral equation and then converted to the eigenvalue problem of linear equations [43]. Therefore, the main obstacle is to select a fast convergent of the basis function family to calculate the shape spectrum for the root canal as a specific boundary condition.The main contribution of this study is to overcome the shortcomings of the previous methods of evaluation of root canal preparation by surface area, volume, or eccentricity thickness of the central line. Based on strict mathematical theories, the proposed methods for comparing shape differences are given in terms of conicity, permeability, and smoothness of the whole root canal. This method does not require the establishment of calibration points for the shape of the root canal. This thought avoids the precondition that two shapes must share the same topological structure in previous algorithms and is very suitable for the shape comparison of the root canal. It also breaks through the problems of different CT sampling rates and inconsistent coordinates of in vivo data acquisition and evaluates root canal preparation more accurately and comprehensively.


## 3. Theoretical Analysis

### 3.1. Main Idea

Thanks to the excellent characteristics of the Steklov spectrum in describing concave and convex shapes and volume changes and its robustness in the calculation, the authors propose an algorithm of root canal treatment evaluation based on the shape difference of the Steklov spectrum. The overall idea is as follows:(1)
Firstly, the 3D models before and after root canal preparationM0 and M1 were acquired by μCT.(2)
To synthesize the changes of a cross section ofM0 and M1 at the same height, the key cross section was selected, and the ideal root canal model M˜1 was generated by lofting.(3)
Calculate the Steklov shape spectrum ofM1 and M˜1, and then get the eigenvalues λi and μj.(4)
Define the spectral distance Steklov spectrum discrete (SSD) to evaluate the difference. The distance is used to measure the shape difference betweenM1 and M˜1 so as to evaluate the quality of root canal preparation.


### 3.2. Mathematical Model of the Steklov Spectrum

In the three-dimensional Euclidean space, the Steklov eigenvalue problem is defined as follows.LetΩ∈ℝ3 be a compact Riemannian manifold, whose boundary Σ=∂Ω satisfies the Lipschitz condition. The Steklov equation is used to find the harmonic function ux and real number λ≠0. The Steklov problem on Ω is(1)Δux=0,x∈Ω,∂ux∂n⟶=λux,x∈Σ,where Δ=∑i=13∂2/∂xi2 is the Laplace operator on functions on Ω and ∂ux/∂n⟶ is the outward derivative along Σ. In the second equation in equation (1), it can be interpreted as an eigenvalue problem that transforms the Dirichlet boundary condition into the Neumann boundary condition. It maps function f∈C∞Σ to Df=∂nHf∈C∞Ω:(2)D:f⟶∂nHf,where D is the elliptic pseudodifferential operator of order one. Hf is the harmonic extension of f on Ω and extends smoothly from the boundary to the interior of Ω [44, 45]. The eigenfunctions ϕ1,ϕ2,…, corresponding to the eigenvalues constitutes the orthogonal basis in L2∑ space. The trace of the eigen matrix formed by ϕi is closely related to the intrinsic geometry and the extrinsic geometry of the shape itself.LetTre−tD be the heat kernel of the operator e−tD:f⟶u for the heat equation ∂tu+Pu=0, when t⟶0+, and there is an asymptotic expansion of Tre−tD [46]:(3)hx,t≔Tre−tD=∫Σe−tDx,xdx∼∑k=0∞akxt−n+1+k+∑l=1∞blxtllogt.The coefficientsa0,a1,…,an−1, as well as b0,b1,…,bn−1, are the Steklov heat invariants which are determined by the Steklov spectrum. According to [26], when n=3 and Ω∈ℝ3,(4)a0x=12π,(5)a1x=KHx4π,where KHx is the mean curvature at x.(6)a2x=116πKHx+RMx6,where RMx is the Gaussian curvature of boundary Σ.From equations (3) to (6), it can be seen that the mean curvatures reflecting the intrinsic geometry and the Gaussian curvature of the extrinsic geometry are determined by the Steklov spectrum.In summary, by calculating the eigenvalues and eigenfunctions of the Steklov operators, the intrinsic geometry and extrinsic geometry of shape can be transformed into the spectrum space for research, and then shape difference evaluation can be carried out. Our thought is based on PDE, which obtains the whole geometric structure from the boundary, which is a powerful mathematical tool for medical CT imaging and geological exploration imaging [47].According to the literature [48], the Steklov spectrum and LBO spectrum of shape are invariant in isometric transformation. Moreover, it has strong robustness for model downsampling, and even nonuniform sampling has better immune effect for topological noise. What is more, the conjugate gradient descent method can be used to solve the above spectrum problem, and the running speed of the algorithm does not depend significantly on the number of edges and faces. The Steklov spectrum has good accuracy even for degenerated irregular meshes. Therefore, it is very beneficial to realize the fast evaluation of root canal shape difference.


### 3.3. Discussion on the Numerical Solution of the Steklov Spectrum

The numerical solution of PDE boundary value problems for high-dimensional manifolds, and two-manifold with higher genus often falls into local minima or diverges, leading to calculation failure. A series of thresholds or carefully adjusted penalty functions were often needed to intervene manually. The 3D topology and geometry of root canals are relatively fixed (tree-shaped cavity), gradually thinning from top to bottom. Therefore, it is possible to avoid the shortcomings of the general algorithm and improve the automaticity and accuracy of the algorithm by thoroughly combining this priori knowledge. Thus, the design of a particular shape difference evaluation algorithm based on relatively fixed root canal topology will be based on the conicity, permeability, and smoothness of the whole root canal shape. The shortcomings of the previous local and single geometric evaluation method were broken through to provide a more scientific and comprehensive assessment method for root canal preparation effect. The calculation method is deduced as follows.The numerical solution is based on potential theory. The potential theory is a robust basis for solving mathematical and physical equations by using the integral equation method. By using the linear superposition principle and fundamental solutions, various forms of integral expressions are obtained. Then, by solving the boundary, the purpose of solving the definite solution problem of the differential equation corresponding to the physical problem is achieved.In this paper, the boundary element method is used to solve the Steklov eigenvalue problem. This method is different from the finite element method, the variational method, and the radial basis function method. It has the characteristics of simplicity, fast convergence, and high accuracy. It is very suitable for the case of complex internal structure and gridding [49]. The basic idea of the boundary element method is that first, the differential equation is transformed into a boundary integral equation. After that, the discrete quadrature method is used for numerical calculation. In solving the discrete matrix elements of the integral equation, the solution needs to be obtained by a singular integral. In this paper, after derivation, the asymptotic expansion of the odd order can be obtained. When the accuracy can be achieved, the convergence speed of the eigenvalue problem can be significantly improved. At the same time, a posterior error estimate can be obtained based on residual estimation [50].According to the potential theory of the differential equation, wheny is a nonsingular point on Σ, the Steklov eigenvalue problem can be converted to the following integral equation with the Green formula [51]:(7)12ux=−∫Σ∂hx,y∂nyuydsy+∫Σhx,yvydsy,where vx=∂ux/∂nx and U∗y,x=1/4πx−y is the fundamental solution of the following equation:(8)LyGx,y=δ0y−x,forx,y∈ℝ3.Thus, the Steklov eigenvalue problem is transformed into the boundary element numerical calculation problem. The first term on the right-hand side of equation (7) is the double-layer potential, the second term is the single-layer potential, and the corresponding density functions are −uy and vx.According to the continuity of single-layer potential on the boundary, equation (7) still holds for x∈Σ. According to the discontinuity of the single-layer potential derivative and the nature of double-layer potential, the following relations can be obtained:(9)∂u∂n=V−1λI+Ku,where V:H1/2Σ⟶H1/2Σ is the single-layer potential operator.(10)Vwx≔∫ΣU∗x,ywydsy,where U∗y,x=1/4πx−y is the fundamental solution of Δux=0. The integral operator V maps the density function wx to the potential distribution on Σ.K is the double-layer potential as follows:(11)ux=Kvx≔∫y∈Σ:y−x≥ε∂U∗x,y∂nyvydsy.Since the right-hand side of equation (9) is not necessarily a positive definite matrix after being discretized into a matrix, the conjugate gradient method cannot be used to calculate its eigenvalues. According to [51], equation (9) can be transformed into another eigenvalue problem with the same solution through the Calderón projection:(12)∂ux∂n=D+σxI+K∗V−1σxI+Kux=λux,forx∈Σ,(13)σx≔limε⟶∞14πε2∫y∈Ω:y−x=εdsy,forx∈Σ,for almost x∈Σ, σx=0.5 [51].In equation (12), D is a hypersingular operator which is defined as(14)Dvx≔−∫Σ∂2U∗x,y∂nx∂nyvydsy.K
∗
:
H
−
1
/
2
Σ
⟶
H
−
1
/
2
Σ is the adjoint double-layer potential, which is defined as the conormal derivative of V:(15)K∗wx=∫Σ∂U∗x,y∂nxwydsy.The potential functions discussed above are integral operators, which are defined in both continuous space and infinite-dimensional function space. However, the root canal model is discretized, so these boundary operators should be further discretized into corresponding matrices which is discussed in Section4.


### 3.4. Constructing an Ideal Root Canal Model as a Benchmark

The ideal model serves as a benchmark for the evaluation of preparatory effect. For clearing the inner wall of the root canal, facilitating subsequent obturation and minimizing the transportation of centerline of the root canal, the prepared root canals should have the following criteria:(1)
The central axis should be consistent with the original one(2)
The new cross section should be enlarged evenly on the basis of the old one to ensure that all the inner walls are treated evenly(3)
The volume of the cut inner cavity should be as small as possible on the premise of satisfying Step (2)(4)
The transition between sections should be as smooth as possibleBased on the four criteria, the second result is better than the first result in Figures6(d) and 6(e).Figure 6
Comparison of two results before and after root canal preparation. (a) Before preparation. (b) After preparation. (c) Compound root canal. (d) Result I at A-A (5x). (e) Result II at A-A (5).
(a)
(b)
(c)
(d)
(e)Since most of the sections are nearly circular, we propose an ideal algorithm of root canal model generation based on the ideal inner diameter. The constructing steps are as follows:(1)
Assuming that the height of the root canal ish, a root canal is divided into three segments: coronal, middle, and apical from top to bottom with heights hc,hm,andha, respectively. For the coronal segment, three sections are intercepted at three heights of hc/3−a,hc/3 and hc/3+a (a=0.05hc). Also in the middle and apical segments, three sections each are intercepted in the same way as in coronal’s. A total of three groups of 9 sections Sici=13,Sjmj=13,andSkak=13 are obtained, and the set of boundary points of the sections are pici=13,pjmj=13,andpkak=13, respectively.(2)
To calculate the ideal inner diameter of the section, the following formulas are used:(16)rideal⋅≔A⋅π,where A⋅ represents the area of the given section.(3)
For each section, the average distanced¯⋅ between the actual boundary and the ideal boundary is calculated by using the following formula:(17)d¯⋅≔∑i=1Npm−c⋅−rideal⋅N,where pi denotes the ith boundary point of the current section and c⋅ is the barycentric coordinates of the boundary points of each section. N represents the number of boundary points of the current section. For each segment of the root canal, the cross section with the minimum d¯⋅ is selected as the representative section of the segment. Obviously, when the section is circular, its d¯⋅ equals zero exactly.(4)
In Step (3), the critical cross sectionsSc,Sm,andSa, of the three segments were obtained. The top section St and the bottom section Sb of the root canal can be achieved from the original model. The ideal root canal model can be generated by lofting in the order of St, Sc,Sm,Sa, and Sb.


## 3.1. Main Idea

Thanks to the excellent characteristics of the Steklov spectrum in describing concave and convex shapes and volume changes and its robustness in the calculation, the authors propose an algorithm of root canal treatment evaluation based on the shape difference of the Steklov spectrum. The overall idea is as follows:(1)
Firstly, the 3D models before and after root canal preparationM0 and M1 were acquired by μCT.(2)
To synthesize the changes of a cross section ofM0 and M1 at the same height, the key cross section was selected, and the ideal root canal model M˜1 was generated by lofting.(3)
Calculate the Steklov shape spectrum ofM1 and M˜1, and then get the eigenvalues λi and μj.(4)
Define the spectral distance Steklov spectrum discrete (SSD) to evaluate the difference. The distance is used to measure the shape difference betweenM1 and M˜1 so as to evaluate the quality of root canal preparation.


## 3.2. Mathematical Model of the Steklov Spectrum

In the three-dimensional Euclidean space, the Steklov eigenvalue problem is defined as follows.LetΩ∈ℝ3 be a compact Riemannian manifold, whose boundary Σ=∂Ω satisfies the Lipschitz condition. The Steklov equation is used to find the harmonic function ux and real number λ≠0. The Steklov problem on Ω is(1)Δux=0,x∈Ω,∂ux∂n⟶=λux,x∈Σ,where Δ=∑i=13∂2/∂xi2 is the Laplace operator on functions on Ω and ∂ux/∂n⟶ is the outward derivative along Σ. In the second equation in equation (1), it can be interpreted as an eigenvalue problem that transforms the Dirichlet boundary condition into the Neumann boundary condition. It maps function f∈C∞Σ to Df=∂nHf∈C∞Ω:(2)D:f⟶∂nHf,where D is the elliptic pseudodifferential operator of order one. Hf is the harmonic extension of f on Ω and extends smoothly from the boundary to the interior of Ω [44, 45]. The eigenfunctions ϕ1,ϕ2,…, corresponding to the eigenvalues constitutes the orthogonal basis in L2∑ space. The trace of the eigen matrix formed by ϕi is closely related to the intrinsic geometry and the extrinsic geometry of the shape itself.LetTre−tD be the heat kernel of the operator e−tD:f⟶u for the heat equation ∂tu+Pu=0, when t⟶0+, and there is an asymptotic expansion of Tre−tD [46]:(3)hx,t≔Tre−tD=∫Σe−tDx,xdx∼∑k=0∞akxt−n+1+k+∑l=1∞blxtllogt.The coefficientsa0,a1,…,an−1, as well as b0,b1,…,bn−1, are the Steklov heat invariants which are determined by the Steklov spectrum. According to [26], when n=3 and Ω∈ℝ3,(4)a0x=12π,(5)a1x=KHx4π,where KHx is the mean curvature at x.(6)a2x=116πKHx+RMx6,where RMx is the Gaussian curvature of boundary Σ.From equations (3) to (6), it can be seen that the mean curvatures reflecting the intrinsic geometry and the Gaussian curvature of the extrinsic geometry are determined by the Steklov spectrum.In summary, by calculating the eigenvalues and eigenfunctions of the Steklov operators, the intrinsic geometry and extrinsic geometry of shape can be transformed into the spectrum space for research, and then shape difference evaluation can be carried out. Our thought is based on PDE, which obtains the whole geometric structure from the boundary, which is a powerful mathematical tool for medical CT imaging and geological exploration imaging [47].According to the literature [48], the Steklov spectrum and LBO spectrum of shape are invariant in isometric transformation. Moreover, it has strong robustness for model downsampling, and even nonuniform sampling has better immune effect for topological noise. What is more, the conjugate gradient descent method can be used to solve the above spectrum problem, and the running speed of the algorithm does not depend significantly on the number of edges and faces. The Steklov spectrum has good accuracy even for degenerated irregular meshes. Therefore, it is very beneficial to realize the fast evaluation of root canal shape difference.


## 3.3. Discussion on the Numerical Solution of the Steklov Spectrum

The numerical solution of PDE boundary value problems for high-dimensional manifolds, and two-manifold with higher genus often falls into local minima or diverges, leading to calculation failure. A series of thresholds or carefully adjusted penalty functions were often needed to intervene manually. The 3D topology and geometry of root canals are relatively fixed (tree-shaped cavity), gradually thinning from top to bottom. Therefore, it is possible to avoid the shortcomings of the general algorithm and improve the automaticity and accuracy of the algorithm by thoroughly combining this priori knowledge. Thus, the design of a particular shape difference evaluation algorithm based on relatively fixed root canal topology will be based on the conicity, permeability, and smoothness of the whole root canal shape. The shortcomings of the previous local and single geometric evaluation method were broken through to provide a more scientific and comprehensive assessment method for root canal preparation effect. The calculation method is deduced as follows.The numerical solution is based on potential theory. The potential theory is a robust basis for solving mathematical and physical equations by using the integral equation method. By using the linear superposition principle and fundamental solutions, various forms of integral expressions are obtained. Then, by solving the boundary, the purpose of solving the definite solution problem of the differential equation corresponding to the physical problem is achieved.In this paper, the boundary element method is used to solve the Steklov eigenvalue problem. This method is different from the finite element method, the variational method, and the radial basis function method. It has the characteristics of simplicity, fast convergence, and high accuracy. It is very suitable for the case of complex internal structure and gridding [49]. The basic idea of the boundary element method is that first, the differential equation is transformed into a boundary integral equation. After that, the discrete quadrature method is used for numerical calculation. In solving the discrete matrix elements of the integral equation, the solution needs to be obtained by a singular integral. In this paper, after derivation, the asymptotic expansion of the odd order can be obtained. When the accuracy can be achieved, the convergence speed of the eigenvalue problem can be significantly improved. At the same time, a posterior error estimate can be obtained based on residual estimation [50].According to the potential theory of the differential equation, wheny is a nonsingular point on Σ, the Steklov eigenvalue problem can be converted to the following integral equation with the Green formula [51]:(7)12ux=−∫Σ∂hx,y∂nyuydsy+∫Σhx,yvydsy,where vx=∂ux/∂nx and U∗y,x=1/4πx−y is the fundamental solution of the following equation:(8)LyGx,y=δ0y−x,forx,y∈ℝ3.Thus, the Steklov eigenvalue problem is transformed into the boundary element numerical calculation problem. The first term on the right-hand side of equation (7) is the double-layer potential, the second term is the single-layer potential, and the corresponding density functions are −uy and vx.According to the continuity of single-layer potential on the boundary, equation (7) still holds for x∈Σ. According to the discontinuity of the single-layer potential derivative and the nature of double-layer potential, the following relations can be obtained:(9)∂u∂n=V−1λI+Ku,where V:H1/2Σ⟶H1/2Σ is the single-layer potential operator.(10)Vwx≔∫ΣU∗x,ywydsy,where U∗y,x=1/4πx−y is the fundamental solution of Δux=0. The integral operator V maps the density function wx to the potential distribution on Σ.K is the double-layer potential as follows:(11)ux=Kvx≔∫y∈Σ:y−x≥ε∂U∗x,y∂nyvydsy.Since the right-hand side of equation (9) is not necessarily a positive definite matrix after being discretized into a matrix, the conjugate gradient method cannot be used to calculate its eigenvalues. According to [51], equation (9) can be transformed into another eigenvalue problem with the same solution through the Calderón projection:(12)∂ux∂n=D+σxI+K∗V−1σxI+Kux=λux,forx∈Σ,(13)σx≔limε⟶∞14πε2∫y∈Ω:y−x=εdsy,forx∈Σ,for almost x∈Σ, σx=0.5 [51].In equation (12), D is a hypersingular operator which is defined as(14)Dvx≔−∫Σ∂2U∗x,y∂nx∂nyvydsy.K
∗
:
H
−
1
/
2
Σ
⟶
H
−
1
/
2
Σ is the adjoint double-layer potential, which is defined as the conormal derivative of V:(15)K∗wx=∫Σ∂U∗x,y∂nxwydsy.The potential functions discussed above are integral operators, which are defined in both continuous space and infinite-dimensional function space. However, the root canal model is discretized, so these boundary operators should be further discretized into corresponding matrices which is discussed in Section4.


## 3.4. Constructing an Ideal Root Canal Model as a Benchmark

The ideal model serves as a benchmark for the evaluation of preparatory effect. For clearing the inner wall of the root canal, facilitating subsequent obturation and minimizing the transportation of centerline of the root canal, the prepared root canals should have the following criteria:(1)
The central axis should be consistent with the original one(2)
The new cross section should be enlarged evenly on the basis of the old one to ensure that all the inner walls are treated evenly(3)
The volume of the cut inner cavity should be as small as possible on the premise of satisfying Step (2)(4)
The transition between sections should be as smooth as possibleBased on the four criteria, the second result is better than the first result in Figures6(d) and 6(e).Figure 6
Comparison of two results before and after root canal preparation. (a) Before preparation. (b) After preparation. (c) Compound root canal. (d) Result I at A-A (5x). (e) Result II at A-A (5).
(a)
(b)
(c)
(d)
(e)Since most of the sections are nearly circular, we propose an ideal algorithm of root canal model generation based on the ideal inner diameter. The constructing steps are as follows:(1)
Assuming that the height of the root canal ish, a root canal is divided into three segments: coronal, middle, and apical from top to bottom with heights hc,hm,andha, respectively. For the coronal segment, three sections are intercepted at three heights of hc/3−a,hc/3 and hc/3+a (a=0.05hc). Also in the middle and apical segments, three sections each are intercepted in the same way as in coronal’s. A total of three groups of 9 sections Sici=13,Sjmj=13,andSkak=13 are obtained, and the set of boundary points of the sections are pici=13,pjmj=13,andpkak=13, respectively.(2)
To calculate the ideal inner diameter of the section, the following formulas are used:(16)rideal⋅≔A⋅π,where A⋅ represents the area of the given section.(3)
For each section, the average distanced¯⋅ between the actual boundary and the ideal boundary is calculated by using the following formula:(17)d¯⋅≔∑i=1Npm−c⋅−rideal⋅N,where pi denotes the ith boundary point of the current section and c⋅ is the barycentric coordinates of the boundary points of each section. N represents the number of boundary points of the current section. For each segment of the root canal, the cross section with the minimum d¯⋅ is selected as the representative section of the segment. Obviously, when the section is circular, its d¯⋅ equals zero exactly.(4)
In Step (3), the critical cross sectionsSc,Sm,andSa, of the three segments were obtained. The top section St and the bottom section Sb of the root canal can be achieved from the original model. The ideal root canal model can be generated by lofting in the order of St, Sc,Sm,Sa, and Sb.


## 4. Implementation of Algorithms

### 4.1. Spectrum Numerical Solution Algorithm

The potential operator is defined in infinite-dimensional function space, and the actual root canal model is a triangular mesh composed of finite points. In this paper, the Galerkin method is used to convert the continuous equation to weak form. Specifically, by defining the basis function, some constraints are imposed on the original function space [42]: finding ϕ∈ℍdΣ in function space ℍ, such as that for all u∈ℍdΣ, aϕ,u=fu∈ℝ, where a⋅,⋅ represents a bilinear form, and it is also a bounded linear function. Therefore, for the weak form of single-layer un, it is expressed by a set of basis function e1,e2,…,en:(18)un=∑j=1naei,ejuj.Thus, the original integral problem on the whole boundary is transformed into the issue of solving weak form linear equations.Let the discretized root canal model beM=V,E,F. Substitute equation (18) into integral equations (equations (10), (11), (14), and (15)). These integral equations will be discretized into summation form and can be expressed in matrix form eventually.Any piecewise linear hat functionϕi:Σ⟶ℝ on M equals one at their associated vertex and zero at all other vertices and is orthogonal at different subscripts. Therefore, any square integrable function u˜ on M can be expressed as(19)u˜=∑iciϕi,ci=∫Σϕi,udx,ci∈ℝ.Assume thatwx∈Σ is the solution of(20)Vwx=fx.Then, the weak form of equation (20) is(21)ϕ,VwΣ=ϕ,fΣ.AfterΣ is discretized into the mesh M, wandf are also discretized into vectors , so the following linear system can be obtained:(22)Vw=Gf,where V=vij,G=gij, and(23)vij=∑T1∈Ti∑T2∈Tj∬T1×T2U∗x,yϕixϕjxdsxdsy,gij=∑T1∈Ti∫Tiϕixϕjxdsx.Similarly, the potential integral equations in equation (12) can be listed as a series of linear systems, namely,(24)D⟶D=dij,K⟶K=kij,K∗⟶K∗=kij∗,where dij is defined as(25)dij=14π∑T1∈Ai∑T2∈Aj∬T1×T23y−x,nxy−x,nyx−y5−nx,nyx−y3ϕixϕjxdsy,and nxandny denote the normal directions.Also in equation (24), kij is defined as(26)kij=14π∑T1∈Ai∑T2∈Aj∬T1×T23x−y,nyx−y3ϕixϕjxdsy,where kij∗ is defined as(27)kij∗=14π∑T1∈Ai∑T2∈Aj∫T1×T23y−x,nxx−y3ϕixϕjxdsy.From the above discussion, we then obtain that the Steklov eigenvalue problem on discrete meshes is(28)D+0.5I+K∗V−10.5I+Ku=λGu.The matricesD,K,K∗,V,andG in equation (28) are symmetrical, and the summation and inverse matrices are symmetrical and positive definite matrices, which can be solved by using the conjugate gradient algorithm. Firstly, we rewrite equation (28) as follows:(29)A˜U¯=0.In this paper, the PCG method is used to solve equation (29). The main idea is to transform the linear system (29) into equivalent functional optimization problems. In view of the convergence of the steepest descent method and the insufficiency of the conjugate gradient method depending on the distribution of eigenvalues, the incomplete Cholesky decomposition conjugate gradient (ICCG) method not only keeps the matrix sparse but also makes the preoptimized matrix approximate identity matrix. Therefore, the method accelerates the iteration process [52].Assume thatLTL is an approximate decomposition of A:(30)A˜=LTL+R,where R is a residual matrix and can be used to control the sparsity of L. Since the elements of L should ensure that the preoptimal matrix is effective, LTL is an approximate decomposition of A. This paper adopts the following strategies.Referring to [53], if the element a˜ij of A˜ is not zero, then the element lij of L is calculated with the complete Cholesky decomposition formula. Otherwise, if a˜ij=0, then let lij=0. After finding the incomplete decomposition of A˜ (let L be nonsingular), A˜=LTL can be taken as the preoptimal matrix. Since L is not singular, A˜ is symmetrically positive definite. Thus, the process of solving the linear system A˜U¯=0 can be converted to solving the following two triangular matrix linear systems:(31)Ly˜=0,LTz˜=y˜.In this work, heuristic strategies were used to detecting the convergence rate at different resolutions to determine the optimal sparse degree of the direction. The specific algorithm is as follows.


### 4.2. Root Canal Shape Evaluation Based on the Spectral Difference

In many kinds of literature, shape difference evaluation based on landmark points is the primary method, including geodesic distance based on point pairs and spectral distance based on point pairs. These methods require that there are well-matched points between two models to be compared, and registration is needed before comparison. Most of the shapes before and after root canal treatment have nonrigid deformation. It is almost impossible to find the corresponding points in the shape of the operation, so the above methods are not suitable for the problems studied in this paper.In the related literature of shape differences, Ovsjanikov et al. [38] presented a novel representation of maps between pairs of shapes, and an energy measuring distortions induced by a map was defined on the function space to measure the difference of shapes. Huang et al. [54] proposed a pipeline for constructing an area-based shape difference operator on point clouds and showed numerically that the results are robust and informative.Suppose the spectrum extracted from the two root canal models isλ=λi and μ=μi. A naive difference evaluation method is to calculate the Euclidean distance between λ and μ as(32)dsM,M˜=λ−μ2.However, the influence of different eigenvalues on shape is different since the first dozens of eigenvalues determined the shape.Kac [17] proposed the spectral distance weighted by i−1 but did not consider the influence factors of the regional characteristics of a specific shape.Based on the spectral theory, we know that the low-frequency eigenvalue reflects the global geometry (that is, the high-frequency signal is cancelled in a flat place) and the high-frequency eigenvalue reflects the local geometry (resonance occurs at the end of the shape). This principle was used in document [19, 21, 55, 56] to realize the interchange of gestures and details of multiple animated characters. Therefore, for each eigenvalue, its contribution to the difference between two shapes is proportional to the size of the variation range of the shape details. Intuitively, the curvature is the amount by which a geometric object such as a surface deviates from being a flat plane. The intensity of detail changes on the mesh can be quantified by the Gaussian curvature. The Gaussian curvature of the triangular mesh surface can be calculated by using the following formula:(33)Kx=1Ax2π−∑xi∈N1xθi,where x∈V and Axis the area of the 1-ring region centered on x. θi is the angle of one of the adjacent triangles. Examples of the Gaussian curvature are shown in Figure 7.Figure 7
An example of calculating the Gauss curvature of the root canal model.Accordingly, we propose a spectral difference evaluation algorithm based on Gaussian curvature statistical distribution. Suppose the mesh model of a root canal isM=V,E,F, Nv represent the number vertex in V. We select the first l eigenvalues in the spectrum to evaluate the shape difference. The spectral difference evaluation result dM,N can be calculated by using Algorithm 1.Algorithm 1: Shape evaluation based on spectral difference and Gauss curvature.
(1)
F
or
all
v in V
(2)
CalculateGauss curvature kx
(3)
E
nd
for
(4)
S
ort
k
x
,
then
get
K
^
=
k
^
j
1
N
x
,
m
k
=
k
^
1
≤
k
^
2
≤
⋯
≤
k
^
N
x
=
M
K.
(5)
d
⟶
M
k
−
m
k
/
l
(6)
I
1
∪
I
2
∪
⋯
∪
I
i
∪
⋯
∪
I
l, whereIi≔mK+i−1d,mK+id
(7)
F
or
each
I
i
(8)
ri⟵Nvi/Nv, wherevi=v∣Kvi∈Ii
(9)
E
nd
for
(10)
d
⟵
0
(11)
F
or
each
λ
i in λi1l
(12)
d
+
=
r
i
λ
i
−
μ
i
2
(13)
E
nd
for
(14)
d
M
,
N
⟵
dThrough the above algorithm, we can quantify the weight of the difference between two spectrums in the overall shape difference comprehensively. Therefore, we can get better experimental results.


## 4.1. Spectrum Numerical Solution Algorithm

The potential operator is defined in infinite-dimensional function space, and the actual root canal model is a triangular mesh composed of finite points. In this paper, the Galerkin method is used to convert the continuous equation to weak form. Specifically, by defining the basis function, some constraints are imposed on the original function space [42]: finding ϕ∈ℍdΣ in function space ℍ, such as that for all u∈ℍdΣ, aϕ,u=fu∈ℝ, where a⋅,⋅ represents a bilinear form, and it is also a bounded linear function. Therefore, for the weak form of single-layer un, it is expressed by a set of basis function e1,e2,…,en:(18)un=∑j=1naei,ejuj.Thus, the original integral problem on the whole boundary is transformed into the issue of solving weak form linear equations.Let the discretized root canal model beM=V,E,F. Substitute equation (18) into integral equations (equations (10), (11), (14), and (15)). These integral equations will be discretized into summation form and can be expressed in matrix form eventually.Any piecewise linear hat functionϕi:Σ⟶ℝ on M equals one at their associated vertex and zero at all other vertices and is orthogonal at different subscripts. Therefore, any square integrable function u˜ on M can be expressed as(19)u˜=∑iciϕi,ci=∫Σϕi,udx,ci∈ℝ.Assume thatwx∈Σ is the solution of(20)Vwx=fx.Then, the weak form of equation (20) is(21)ϕ,VwΣ=ϕ,fΣ.AfterΣ is discretized into the mesh M, wandf are also discretized into vectors , so the following linear system can be obtained:(22)Vw=Gf,where V=vij,G=gij, and(23)vij=∑T1∈Ti∑T2∈Tj∬T1×T2U∗x,yϕixϕjxdsxdsy,gij=∑T1∈Ti∫Tiϕixϕjxdsx.Similarly, the potential integral equations in equation (12) can be listed as a series of linear systems, namely,(24)D⟶D=dij,K⟶K=kij,K∗⟶K∗=kij∗,where dij is defined as(25)dij=14π∑T1∈Ai∑T2∈Aj∬T1×T23y−x,nxy−x,nyx−y5−nx,nyx−y3ϕixϕjxdsy,and nxandny denote the normal directions.Also in equation (24), kij is defined as(26)kij=14π∑T1∈Ai∑T2∈Aj∬T1×T23x−y,nyx−y3ϕixϕjxdsy,where kij∗ is defined as(27)kij∗=14π∑T1∈Ai∑T2∈Aj∫T1×T23y−x,nxx−y3ϕixϕjxdsy.From the above discussion, we then obtain that the Steklov eigenvalue problem on discrete meshes is(28)D+0.5I+K∗V−10.5I+Ku=λGu.The matricesD,K,K∗,V,andG in equation (28) are symmetrical, and the summation and inverse matrices are symmetrical and positive definite matrices, which can be solved by using the conjugate gradient algorithm. Firstly, we rewrite equation (28) as follows:(29)A˜U¯=0.In this paper, the PCG method is used to solve equation (29). The main idea is to transform the linear system (29) into equivalent functional optimization problems. In view of the convergence of the steepest descent method and the insufficiency of the conjugate gradient method depending on the distribution of eigenvalues, the incomplete Cholesky decomposition conjugate gradient (ICCG) method not only keeps the matrix sparse but also makes the preoptimized matrix approximate identity matrix. Therefore, the method accelerates the iteration process [52].Assume thatLTL is an approximate decomposition of A:(30)A˜=LTL+R,where R is a residual matrix and can be used to control the sparsity of L. Since the elements of L should ensure that the preoptimal matrix is effective, LTL is an approximate decomposition of A. This paper adopts the following strategies.Referring to [53], if the element a˜ij of A˜ is not zero, then the element lij of L is calculated with the complete Cholesky decomposition formula. Otherwise, if a˜ij=0, then let lij=0. After finding the incomplete decomposition of A˜ (let L be nonsingular), A˜=LTL can be taken as the preoptimal matrix. Since L is not singular, A˜ is symmetrically positive definite. Thus, the process of solving the linear system A˜U¯=0 can be converted to solving the following two triangular matrix linear systems:(31)Ly˜=0,LTz˜=y˜.In this work, heuristic strategies were used to detecting the convergence rate at different resolutions to determine the optimal sparse degree of the direction. The specific algorithm is as follows.


## 4.2. Root Canal Shape Evaluation Based on the Spectral Difference

In many kinds of literature, shape difference evaluation based on landmark points is the primary method, including geodesic distance based on point pairs and spectral distance based on point pairs. These methods require that there are well-matched points between two models to be compared, and registration is needed before comparison. Most of the shapes before and after root canal treatment have nonrigid deformation. It is almost impossible to find the corresponding points in the shape of the operation, so the above methods are not suitable for the problems studied in this paper.In the related literature of shape differences, Ovsjanikov et al. [38] presented a novel representation of maps between pairs of shapes, and an energy measuring distortions induced by a map was defined on the function space to measure the difference of shapes. Huang et al. [54] proposed a pipeline for constructing an area-based shape difference operator on point clouds and showed numerically that the results are robust and informative.Suppose the spectrum extracted from the two root canal models isλ=λi and μ=μi. A naive difference evaluation method is to calculate the Euclidean distance between λ and μ as(32)dsM,M˜=λ−μ2.However, the influence of different eigenvalues on shape is different since the first dozens of eigenvalues determined the shape.Kac [17] proposed the spectral distance weighted by i−1 but did not consider the influence factors of the regional characteristics of a specific shape.Based on the spectral theory, we know that the low-frequency eigenvalue reflects the global geometry (that is, the high-frequency signal is cancelled in a flat place) and the high-frequency eigenvalue reflects the local geometry (resonance occurs at the end of the shape). This principle was used in document [19, 21, 55, 56] to realize the interchange of gestures and details of multiple animated characters. Therefore, for each eigenvalue, its contribution to the difference between two shapes is proportional to the size of the variation range of the shape details. Intuitively, the curvature is the amount by which a geometric object such as a surface deviates from being a flat plane. The intensity of detail changes on the mesh can be quantified by the Gaussian curvature. The Gaussian curvature of the triangular mesh surface can be calculated by using the following formula:(33)Kx=1Ax2π−∑xi∈N1xθi,where x∈V and Axis the area of the 1-ring region centered on x. θi is the angle of one of the adjacent triangles. Examples of the Gaussian curvature are shown in Figure 7.Figure 7
An example of calculating the Gauss curvature of the root canal model.Accordingly, we propose a spectral difference evaluation algorithm based on Gaussian curvature statistical distribution. Suppose the mesh model of a root canal isM=V,E,F, Nv represent the number vertex in V. We select the first l eigenvalues in the spectrum to evaluate the shape difference. The spectral difference evaluation result dM,N can be calculated by using Algorithm 1.Algorithm 1: Shape evaluation based on spectral difference and Gauss curvature.
(1)
F
or
all
v in V
(2)
CalculateGauss curvature kx
(3)
E
nd
for
(4)
S
ort
k
x
,
then
get
K
^
=
k
^
j
1
N
x
,
m
k
=
k
^
1
≤
k
^
2
≤
⋯
≤
k
^
N
x
=
M
K.
(5)
d
⟶
M
k
−
m
k
/
l
(6)
I
1
∪
I
2
∪
⋯
∪
I
i
∪
⋯
∪
I
l, whereIi≔mK+i−1d,mK+id
(7)
F
or
each
I
i
(8)
ri⟵Nvi/Nv, wherevi=v∣Kvi∈Ii
(9)
E
nd
for
(10)
d
⟵
0
(11)
F
or
each
λ
i in λi1l
(12)
d
+
=
r
i
λ
i
−
μ
i
2
(13)
E
nd
for
(14)
d
M
,
N
⟵
dThrough the above algorithm, we can quantify the weight of the difference between two spectrums in the overall shape difference comprehensively. Therefore, we can get better experimental results.


## 5. Experiment and Results

Traditional evaluation indicators for the effectiveness of systems include volume and area difference, as well as the transportation of the root canal. In this paper, the conicity and permeability of root canals are evaluated from the spectral distance between shapes. To access the superiority of the algorithm, the accuracy and the robustness of the algorithm are evaluated.


### 5.1. Programming

This paper mainly deals with the numerical solution of the integral equation. The most crucial problem is to find the eigenvalues of a matrix, i.e., to solve large-scale sparse linear system. Therefore, we choose the popular eigen algorithm library to calculate. It provides a series of built-in solver libraries and encapsulates some external solver libraries, including LLT and LDLT used in this paper. The algorithm proposed in this paper was implemented through C++ programming with the hardware and software shown in Table1.Table 1
Simulation environment.
Hardware types
Environment Info
CPU
Intel Core i7-3820QM CPU, 2.70 GHz
DRAM
DDR3 1600 MHz
Hard drive
2.5″ 7 mm 128 GB, 1TB, SCSI
OS
Microsoft Windows10 64 bit
IDE
Visual Studio 2012, C++
Interface design
QT4.8
Dependency library
VCGLIB/Meshlab 1.3/eigen3.0


### 5.2. Accuracy Evaluation

#### 5.2.1. Comparison with LBO-Based Method

The Laplace–Beltrami operator is an intrinsic invariant operator whose spectrum also reflects the geometric and topological characteristics of some shapes, so it is often used for shape difference comparison, such as [57]. Most of the models processed by previous works are mesh models with low genus or high regularity. However, the mesh generated from the point cloud of the root canal obtained from μCT is often of poor quality. Therefore, for some spectral methods with low robustness, there may be more problems in the aspect of accuracy.In this paper, the Steklov spectrum method and LBO spectrum method are used to calculate three groups of root canal models, respectively, and the results are shown in Figure8. Our algorithm has an apparent distinction between the apparent uneven root canal shape, the rough root canal, and the different topology, while the performance of the LBO method is not ideal. Therefore, we infer that our algorithm is more suitable for the root canal model with complex geometry and topology in these cases.Figure 8
Comparison of the Steklov spectrum and LBO spectrum of the root canal model. Yellow indicates the original root canal, and grey indicates the root canal after preparation in all graphs and charts.


#### 5.2.2. Comparison with Volume/Area Evaluation Method

Our method is to compare the prepared root canal with the idealized root canal and then calculate the overall geometric structure and topology. Shape difference based on the extrinsic spectrum accurately perceives the concavity and convexity of the shape. So our method is volume aware, area aware, and slice aware in real sense.To evaluate the comparative effect of SSD, we infer that the best root canal shape is the smallest distance from the ideal root canal shape in theory, and the best therapy result is obtained. In [32], spiral computerized tomography (SCT) was reported to be a useful tool in various in vivo and laboratory studies. It was concluded that with SCT, three-dimensional volume measurements are possible without sectioning specimens, thus avoiding loss of material.Specifically, the comparative method is as follows: A total of 60 extracted first molars were allocated into three groups randomly. Each group was treated with one instrument, respectively. The effect of root canal preparation was evaluated by the literature [2, 14, 58], which was based on the percentage of void (POV). The lower the POV, the better the obturation effect.The aim of this comparison is to assess the correctness of the volume and area evaluation method and our method. Statistical analysis was performed with ANOVA (Scheffé tests).The results of experiments summarized in Table2 show that there is no significant difference among the three groups of results obtained by the volume-based discriminant method, and the third group which is considered to be the best one is not consistent with the ground truth evaluation result. Therefore, the approach is considered to be unsuccessful.Table 2
The evaluation comparison among volume method, area method, and our method for three systems (means ± SD,n = 20).
System
A
B
C
Best performance system
Consistency with ground truth evaluation result
Volume method
Mean volume of root canal (mm3)
Before
3.68 ± 0.92
5.02 ± 1.42
4.18 ± 1.82
SystemC
F
After
4.12 ± 0.29
5.99 ± 2.01
5.29 ± 1.12
Diff.
0.44b
0.89a
1.11
a
∗
Area method
Mean area of root canal (mm2)
Before
39 ± 3.2
36.4 ± 3.3
33.6 ± 4.2
SystemB
F
After
37.2 ± 4.0
34.9 ± 3.6
32.2 ± 3.2
Diff.
−1.8a
−1.5a
−1.9a∗
Our method
Mean SSD of root canal
Before
5.9 ± 0.6
6.5 ± 0.7
7.8 ± 1.0
SystemA
T
After
1.1 ± 0.2a∗
1.7 ± 0.3b
2.3 ± 0.3b
Ground truth result
Mean POV
Coronal third
2.1 ± 0.4a∗
3.8 ± 0.6b
2.9 ± 0.6b
SystemA
N/A
Middle third
1.8 ± 0.2a∗
2.3 ± 0.5b
2.1 ± 0.5b
SystemA
Apical third
3.5 ± 0.5a∗
6.9 ± 0.8b
5.7 ± 0.9b
SystemA
∗The minimum shape difference in three systems. Cells with the same letter in the superscript indicate no statistical significance between different instruments (ANOVA, Scheffé tests, P>0.05).The results obtained by the area-based discriminant method did not express a significant difference between the second and the third groups, and the best (second) group is also inconsistent with the ground truth evaluation result either.The evaluation results of the SSD method have achieved good efficacy in two respects. First, only the SSD method is consistent with the ground truth evaluation result. Second, among the three methods, only SSD is capable of distinguishing the significant difference of SystemA from the other. Therefore, it can be concluded that the new method can evaluate the effect of root canal preparation comprehensively and accurately. Hence, the experiment achieves the desired purpose in the dataset.


### 5.3. Robustness Evaluation

In Figure9, it shows an example of the robustness of our method to different sampling density. We compare the eigenvalues of a series of the root canal model with different resolution. It is observed that some geometric details in some regions are lost during the downsampling and the spectrum remains almost unaffected. According to the theory of the Steklov spectrum, the first dozens of eigenvalues and eigenfunctions are robust to downsampling for a mesh with 10k to 100k vertices. Moreover, the number of vertices is generally sufficient for shape difference evaluation.Figure 9
Three root canal meshes represented with 10K (blue), 50K (red), and 100K (yellow) vertices share similar Steklov spectrums.


## 5.1. Programming

This paper mainly deals with the numerical solution of the integral equation. The most crucial problem is to find the eigenvalues of a matrix, i.e., to solve large-scale sparse linear system. Therefore, we choose the popular eigen algorithm library to calculate. It provides a series of built-in solver libraries and encapsulates some external solver libraries, including LLT and LDLT used in this paper. The algorithm proposed in this paper was implemented through C++ programming with the hardware and software shown in Table1.Table 1
Simulation environment.
Hardware types
Environment Info
CPU
Intel Core i7-3820QM CPU, 2.70 GHz
DRAM
DDR3 1600 MHz
Hard drive
2.5″ 7 mm 128 GB, 1TB, SCSI
OS
Microsoft Windows10 64 bit
IDE
Visual Studio 2012, C++
Interface design
QT4.8
Dependency library
VCGLIB/Meshlab 1.3/eigen3.0


## 5.2. Accuracy Evaluation

### 5.2.1. Comparison with LBO-Based Method

The Laplace–Beltrami operator is an intrinsic invariant operator whose spectrum also reflects the geometric and topological characteristics of some shapes, so it is often used for shape difference comparison, such as [57]. Most of the models processed by previous works are mesh models with low genus or high regularity. However, the mesh generated from the point cloud of the root canal obtained from μCT is often of poor quality. Therefore, for some spectral methods with low robustness, there may be more problems in the aspect of accuracy.In this paper, the Steklov spectrum method and LBO spectrum method are used to calculate three groups of root canal models, respectively, and the results are shown in Figure8. Our algorithm has an apparent distinction between the apparent uneven root canal shape, the rough root canal, and the different topology, while the performance of the LBO method is not ideal. Therefore, we infer that our algorithm is more suitable for the root canal model with complex geometry and topology in these cases.Figure 8
Comparison of the Steklov spectrum and LBO spectrum of the root canal model. Yellow indicates the original root canal, and grey indicates the root canal after preparation in all graphs and charts.


### 5.2.2. Comparison with Volume/Area Evaluation Method

Our method is to compare the prepared root canal with the idealized root canal and then calculate the overall geometric structure and topology. Shape difference based on the extrinsic spectrum accurately perceives the concavity and convexity of the shape. So our method is volume aware, area aware, and slice aware in real sense.To evaluate the comparative effect of SSD, we infer that the best root canal shape is the smallest distance from the ideal root canal shape in theory, and the best therapy result is obtained. In [32], spiral computerized tomography (SCT) was reported to be a useful tool in various in vivo and laboratory studies. It was concluded that with SCT, three-dimensional volume measurements are possible without sectioning specimens, thus avoiding loss of material.Specifically, the comparative method is as follows: A total of 60 extracted first molars were allocated into three groups randomly. Each group was treated with one instrument, respectively. The effect of root canal preparation was evaluated by the literature [2, 14, 58], which was based on the percentage of void (POV). The lower the POV, the better the obturation effect.The aim of this comparison is to assess the correctness of the volume and area evaluation method and our method. Statistical analysis was performed with ANOVA (Scheffé tests).The results of experiments summarized in Table2 show that there is no significant difference among the three groups of results obtained by the volume-based discriminant method, and the third group which is considered to be the best one is not consistent with the ground truth evaluation result. Therefore, the approach is considered to be unsuccessful.Table 2
The evaluation comparison among volume method, area method, and our method for three systems (means ± SD,n = 20).
System
A
B
C
Best performance system
Consistency with ground truth evaluation result
Volume method
Mean volume of root canal (mm3)
Before
3.68 ± 0.92
5.02 ± 1.42
4.18 ± 1.82
SystemC
F
After
4.12 ± 0.29
5.99 ± 2.01
5.29 ± 1.12
Diff.
0.44b
0.89a
1.11
a
∗
Area method
Mean area of root canal (mm2)
Before
39 ± 3.2
36.4 ± 3.3
33.6 ± 4.2
SystemB
F
After
37.2 ± 4.0
34.9 ± 3.6
32.2 ± 3.2
Diff.
−1.8a
−1.5a
−1.9a∗
Our method
Mean SSD of root canal
Before
5.9 ± 0.6
6.5 ± 0.7
7.8 ± 1.0
SystemA
T
After
1.1 ± 0.2a∗
1.7 ± 0.3b
2.3 ± 0.3b
Ground truth result
Mean POV
Coronal third
2.1 ± 0.4a∗
3.8 ± 0.6b
2.9 ± 0.6b
SystemA
N/A
Middle third
1.8 ± 0.2a∗
2.3 ± 0.5b
2.1 ± 0.5b
SystemA
Apical third
3.5 ± 0.5a∗
6.9 ± 0.8b
5.7 ± 0.9b
SystemA
∗The minimum shape difference in three systems. Cells with the same letter in the superscript indicate no statistical significance between different instruments (ANOVA, Scheffé tests, P>0.05).The results obtained by the area-based discriminant method did not express a significant difference between the second and the third groups, and the best (second) group is also inconsistent with the ground truth evaluation result either.The evaluation results of the SSD method have achieved good efficacy in two respects. First, only the SSD method is consistent with the ground truth evaluation result. Second, among the three methods, only SSD is capable of distinguishing the significant difference of SystemA from the other. Therefore, it can be concluded that the new method can evaluate the effect of root canal preparation comprehensively and accurately. Hence, the experiment achieves the desired purpose in the dataset.


## 5.2.1. Comparison with LBO-Based Method

The Laplace–Beltrami operator is an intrinsic invariant operator whose spectrum also reflects the geometric and topological characteristics of some shapes, so it is often used for shape difference comparison, such as [57]. Most of the models processed by previous works are mesh models with low genus or high regularity. However, the mesh generated from the point cloud of the root canal obtained from μCT is often of poor quality. Therefore, for some spectral methods with low robustness, there may be more problems in the aspect of accuracy.In this paper, the Steklov spectrum method and LBO spectrum method are used to calculate three groups of root canal models, respectively, and the results are shown in Figure8. Our algorithm has an apparent distinction between the apparent uneven root canal shape, the rough root canal, and the different topology, while the performance of the LBO method is not ideal. Therefore, we infer that our algorithm is more suitable for the root canal model with complex geometry and topology in these cases.Figure 8
Comparison of the Steklov spectrum and LBO spectrum of the root canal model. Yellow indicates the original root canal, and grey indicates the root canal after preparation in all graphs and charts.


## 5.2.2. Comparison with Volume/Area Evaluation Method

Our method is to compare the prepared root canal with the idealized root canal and then calculate the overall geometric structure and topology. Shape difference based on the extrinsic spectrum accurately perceives the concavity and convexity of the shape. So our method is volume aware, area aware, and slice aware in real sense.To evaluate the comparative effect of SSD, we infer that the best root canal shape is the smallest distance from the ideal root canal shape in theory, and the best therapy result is obtained. In [32], spiral computerized tomography (SCT) was reported to be a useful tool in various in vivo and laboratory studies. It was concluded that with SCT, three-dimensional volume measurements are possible without sectioning specimens, thus avoiding loss of material.Specifically, the comparative method is as follows: A total of 60 extracted first molars were allocated into three groups randomly. Each group was treated with one instrument, respectively. The effect of root canal preparation was evaluated by the literature [2, 14, 58], which was based on the percentage of void (POV). The lower the POV, the better the obturation effect.The aim of this comparison is to assess the correctness of the volume and area evaluation method and our method. Statistical analysis was performed with ANOVA (Scheffé tests).The results of experiments summarized in Table2 show that there is no significant difference among the three groups of results obtained by the volume-based discriminant method, and the third group which is considered to be the best one is not consistent with the ground truth evaluation result. Therefore, the approach is considered to be unsuccessful.Table 2
The evaluation comparison among volume method, area method, and our method for three systems (means ± SD,n = 20).
System
A
B
C
Best performance system
Consistency with ground truth evaluation result
Volume method
Mean volume of root canal (mm3)
Before
3.68 ± 0.92
5.02 ± 1.42
4.18 ± 1.82
SystemC
F
After
4.12 ± 0.29
5.99 ± 2.01
5.29 ± 1.12
Diff.
0.44b
0.89a
1.11
a
∗
Area method
Mean area of root canal (mm2)
Before
39 ± 3.2
36.4 ± 3.3
33.6 ± 4.2
SystemB
F
After
37.2 ± 4.0
34.9 ± 3.6
32.2 ± 3.2
Diff.
−1.8a
−1.5a
−1.9a∗
Our method
Mean SSD of root canal
Before
5.9 ± 0.6
6.5 ± 0.7
7.8 ± 1.0
SystemA
T
After
1.1 ± 0.2a∗
1.7 ± 0.3b
2.3 ± 0.3b
Ground truth result
Mean POV
Coronal third
2.1 ± 0.4a∗
3.8 ± 0.6b
2.9 ± 0.6b
SystemA
N/A
Middle third
1.8 ± 0.2a∗
2.3 ± 0.5b
2.1 ± 0.5b
SystemA
Apical third
3.5 ± 0.5a∗
6.9 ± 0.8b
5.7 ± 0.9b
SystemA
∗The minimum shape difference in three systems. Cells with the same letter in the superscript indicate no statistical significance between different instruments (ANOVA, Scheffé tests, P>0.05).The results obtained by the area-based discriminant method did not express a significant difference between the second and the third groups, and the best (second) group is also inconsistent with the ground truth evaluation result either.The evaluation results of the SSD method have achieved good efficacy in two respects. First, only the SSD method is consistent with the ground truth evaluation result. Second, among the three methods, only SSD is capable of distinguishing the significant difference of SystemA from the other. Therefore, it can be concluded that the new method can evaluate the effect of root canal preparation comprehensively and accurately. Hence, the experiment achieves the desired purpose in the dataset.


## 5.3. Robustness Evaluation

In Figure9, it shows an example of the robustness of our method to different sampling density. We compare the eigenvalues of a series of the root canal model with different resolution. It is observed that some geometric details in some regions are lost during the downsampling and the spectrum remains almost unaffected. According to the theory of the Steklov spectrum, the first dozens of eigenvalues and eigenfunctions are robust to downsampling for a mesh with 10k to 100k vertices. Moreover, the number of vertices is generally sufficient for shape difference evaluation.Figure 9
Three root canal meshes represented with 10K (blue), 50K (red), and 100K (yellow) vertices share similar Steklov spectrums.


## 6. Conclusions and Prospects

In this study, we addressed the problem of evaluating the operation result of a root canal by introducing a spectrum-based procedure called the Steklov spectrum distance, which measures the shape difference between the original shape and its ideal shape with its Steklov spectrum. We find it remarkable that the use of SSD simplifies the problem of evaluating the intrinsic and extrinsic difference between shapes, allowing to significantly improve the accuracy of standard pipelines comparing with the traditional volume-area literature. SSD appears to be a valuable tool to assess the efficacy of various root canal preparation systems.Within the limitations of this study, one is that the error of SSD was seen in all groups, and it even diverges in the calculation of certain root canals. To solve these problems, we think that it can be optimized in two ways: one is to design a more expressive distance norm according to the specific shape of root canal; the other is to improve the numerical solution method of partial differential equation, using better presolver and matrix block parallel algorithm to improve the accuracy and efficiency of the algorithm. Another shortage to be improved is that, limited by theμCT itself, this study can only be performed in vitro, rather than applied to direct clinical surgical evaluation. In the future, we will continue to consider the related research based on the data sources of high-speed acquisition of CT images in vivo.



---

*Source: 4830914-2019-11-03.xml*