Method and device for extracting root canal feature morphology based on CBCT data, equipment and medium

Abstract

Embodiments of the present application disclose a method and device for extracting root canal feature morphology based on CBCT data, equipment and medium, belong to the field of oral medicine and medical image technology, wherein the method comprises: filtering and denoising CBCT images and performing pulp segmentation processing, and performing surface rendering three-dimensional reconstruction on tooth root canals to generate a three-dimensional model; based on the regular point cloud data in the three-dimensional model, a point-surface method is used to realize surface reconstruction to generate a root canal entity model; the root canal entity model is subjected to coordinate alignment processing in a three-dimensional space, and the approximate center line of the fitted root canal is used to measure the root canal length, root canal opening direction, root canal curvature and pulp depth. The present application quickly reconstructs the corresponding three-dimensional model based on the CBCT tomographic images obtained by scanning the teeth of a patient, and obtains relevant data such as pulp depth and root canal opening size from the three-dimensional model of the teeth to assist the performance of oral surgery.

Description

Technical Field

[0001] This application relates to a method, device, equipment, and medium for extracting root canal features based on CBCT data, belonging to the fields of oral medicine and medical imaging technology. Background Technology

[0002] Accurate measurement of dental pulp morphological data is a crucial part of clinical oral surgery. However, tooth structures vary greatly among individuals due to differences in age, sex, and ethnicity, especially the highly complex structures of the pulp and root canals. This morphological data not only provides surgeons with preoperative risk assessment and efficient surgical planning—such as determining the optimal location and direction for access, assessing dentin fracture risk, and locating root canal orifices—but also provides intraoperative feedback on surgical progress and related information, including patient condition updates, visual information about the procedure, and biomechanical analysis of instruments operating within the pulp.

[0003] CBCT data is used to display the entire tooth structure and diseased tissue in three dimensions, helping doctors make accurate diagnoses and treatments. The process of acquiring CBCT data mainly includes the following steps:

[0004] Step 1, X-ray irradiation: CBCT uses KV-level X-rays for irradiation. The X-ray source is a three-dimensional cone beam X-ray, which is used to perform a ring-shaped digital irradiation around the object being irradiated.

[0005] Step 2, Data Acquisition: Collect two-dimensional projection data from multiple angles using a flat panel detector.

[0006] Step 3, Data Reconstruction: The collected two-dimensional projection data is reconstructed by computer, and cone-beam CT reconstruction algorithm (such as FDK algorithm) is used to reconstruct the two-dimensional projection data into a three-dimensional image.

[0007] Existing CBCT data cannot directly reflect the morphological features of tooth root canals. Therefore, how to extract morphological features of tooth root canals based on CBCT data to assist in oral surgery is a pressing technical problem that needs to be solved. Summary of the Invention

[0008] To address the aforementioned technical problems, embodiments of this application provide a method, apparatus, device, and medium for extracting root canal features based on CBCT data. This allows dentists to obtain relevant data such as pulp depth and root canal orifice size from a three-dimensional model of the tooth, which can help determine the pulpotomy plan.

[0009] Other features and advantages of this application will become apparent from the following detailed description, or may be learned in part from practice of this application.

[0010] According to one aspect of the embodiments of this application, a method for extracting root canal features based on CBCT data is provided, the method comprising:

[0011] The CBCT images were filtered and denoised, and the pulp was segmented. The root canals of the teeth were then rendered and reconstructed in three dimensions to generate a three-dimensional model.

[0012] Based on the regular point cloud data in the three-dimensional model, the point-surface method is used to reconstruct the surface and generate the root canal solid model.

[0013] The root canal physical model is aligned in three-dimensional space, and the morphological characteristics of the root canal are measured by fitting the approximate centerline of the root canal; wherein, the morphological characteristics of the root canal include at least one of root canal length, root canal orifice direction, root canal curvature, and pulp depth.

[0014] Furthermore, the CBCT images are filtered and denoised, and the pulp is segmented. The root canals are then reconstructed using surface mapping to generate a 3D model, including:

[0015] Based on the acquired CBCT images, the neighborhood edge model is used to improve the denoising effect of median filtering, and the optimal threshold for pulp segmentation is calculated based on the intercept histogram of the inverse cross-entropy, thus obtaining the preprocessed CBCT sequence images.

[0016] The surface model of the root canal is drawn by rendering the preprocessed CBCT sequence images using a surface rendering method.

[0017] Furthermore, based on the acquired CBCT images, a neighborhood edge model is used to improve the denoising effect of median filtering, including:

[0018] Within a neighborhood, eight edge models are established for the center pixel of the neighborhood, and the pixel difference between the center pixel and its neighboring pixels is calculated using the following formula:

[0019] |M i -M c |=d(i=1,2,…,7,8) (1)

[0020] In the formula, M i Let M represent the i-th neighboring pixel, where i is a positive integer not greater than 8. The eight neighboring pixels are located above, below, left, right, and at four 45° angles from the center pixel. c d represents the center pixel, and d represents the pixel difference between the center pixel and the i-th neighboring pixel.

[0021] Set the threshold size to T1, and record the number of adjacent points with gray levels similar to the central pixel as n. If the pixel difference d between the central pixel and the adjacent points is less than T1, then determine that the adjacent points are similar in gray level to the central pixel; let n = n + 1. When all adjacent points have been traversed, if n min ≤n≤n max , b min is the minimum value of the number of adjacent points, and n max is the maximum value of the number of adjacent points, then determine that the number of adjacent points conforms to the established edge model, and determine that the central pixel is an edge point and directly retain it;

[0022] For non-edge points, replace the non-edge points with the median value of all pixel points in the neighborhood. Calculate the median value of all pixel points through the following formula:

[0023] F med =medM (x,y) =med[M (x+m,y+n) ;(m,n)∈A] (2)

[0024] In the formula, F med is the median value of all pixel points, A is the neighborhood around the central pixel (x, y), m and n are the length and width of the neighborhood window, M (x,y) is the pixel value of the pixel point (x, y), med is the median operation, and M (x+m,y+n) is the pixel value of the pixel point (x + m, y + n);

[0025] Sort all the pixel values in the neighborhood to obtain the minimum value F min , median value F med , and maximum value F max ;

[0026] Execute the first step: If d1 = F med -F min > 0 and d2 = F max -F med > 0, then execute the third step; if not satisfied, then execute the second step;

[0027] The second step includes: expanding the neighborhood area by one circle. If the neighborhood area meets the second condition, then sort all the pixel values in the neighborhood again, update the minimum value F min , median value F med , and maximum value F max and repeat the first step. If the neighborhood area does not meet the second condition, then do not replace the gray level value of the central pixel; where the second condition is that the neighborhood area is less than the maximum neighborhood threshold S max ;

[0028] The third step includes: If g1 = Mc -F min >0 and g2=F max -M c If the grayscale value is greater than 0, then the grayscale value of the center pixel is not replaced; otherwise, the grayscale value of the center pixel is replaced with the median value F. med .

[0029] Furthermore, the optimal threshold for pulp segmentation is calculated based on the intercept histogram of the inverse cross-entropy, resulting in preprocessed CBCT sequence images, including:

[0030] Based on the denoised CBCT images, gamma transform is used to enhance image contrast, and the neighborhood average image and gradient image are calculated respectively.

[0031] The neighborhood average image and gradient image are combined with the grayscale image and added pixel by pixel to form a hybrid image;

[0032] The optimal threshold for pulp segmentation is calculated using the intercept histogram of the mixed image;

[0033] Each pixel value in the blended image is compared with the optimal threshold and converted into a binary image;

[0034] Remove all unnecessary connected tissues on the image boundaries of the binary image, retaining only the pulp region, to obtain a preprocessed CBCT sequence image;

[0035] The step of calculating the optimal threshold for pulp segmentation using the intercept histogram of the mixed image includes:

[0036] A three-dimensional space is generated by grayscale information, neighborhood average grayscale information and gradient compensation. The three-dimensional space is divided into target region O and background region B by a plane.

[0037] Based on the aforementioned three-dimensional space, the final gradient is calculated using the following formula:

[0038]

[0039] In the formula, G H (x,y) and G V (x,y) represents the gradient values ​​of a point (x,y) in the horizontal and vertical directions in a given image, where x represents the x-coordinate of the point in the given image, y represents the y-coordinate of the point in the given image, G(x,y) represents the final gradient, g(x,y+1) represents the gray value of point (x,y+1), g(x,y-1) represents the gray value of point (x,y-1), g(x+1,y) represents the gray value of point (x+1,y), and g(x-1,y) represents the gray value of point (x+1,y).

[0040] The hybrid image is represented as:

[0041] F(x,y)=f(x,y)+g(x,y)-G(x,y) (5)

[0042] In the formula, F(x,y) represents the mixed image, and g(x,y) and f(x,y) represent the gray value of point (x,y) and the average gray value of the neighborhood of point (x,y).

[0043] Based on the segmentation threshold T, the hybrid image is divided into the target portion Ω. O ∈{(x,y)|F(x,y)=0,1,…,T} and background part Ω B ∈{(x,y)|F(x,y)=T+1,T+2,…,2T-2}, pixel values ​​less than the segmentation threshold T belong to the target part, and pixel values ​​greater than the segmentation threshold T belong to the background part;

[0044] The inverse cross-entropy between the target region and the background region is calculated using the following formula:

[0045]

[0046] In the formula, E(O,B) represents the reciprocal cross-entropy; P(k) represents the prior probability that the average gray value of the neighborhood of the mixed image is k. This represents the average gray level of the target region. This represents the average grayscale value of the background area;

[0047] Based on the calculated inverse cross-entropy E(O,B), the optimal threshold for pulp segmentation is determined using the following formula:

[0048]

[0049] In the formula, The optimal threshold for pulp segmentation is given by arg min, which is the index operation for finding the minimum element in the array.

[0050] Furthermore, the surface model of the root canal is drawn using surface rendering of the preprocessed CBCT sequence images, including:

[0051] In three-dimensional space, uniform sampling is performed on the upper and lower jaw regions in the x, y, and z directions with sampling intervals of Δx, Δy, and Δz. The volume data is represented by a ternary function; where a cubic region consisting of eight adjacent sampling points is a voxel.

[0052] For each voxel, an isosurface is extracted; wherein, the isosurface is a surface composed of points with the same properties in three-dimensional space;

[0053] Traverse the entire volume data to find the voxels containing isosurfaces, and determine the intersection points P(x,y,z) of the isosurfaces and the edges of the voxels using the following formula:

[0054]

[0055] In the formula, T is the isosurface threshold, M1 and M2 are the gray values ​​of the first and second vertices on the edge where the voxel intersects with the isosurface, x1, y1, z1 are the position coordinates of the first vertex on the edge where the voxel intersects with the isosurface, x2, y2, z2 are the position coordinates of the second vertex on the edge where the voxel intersects with the isosurface, and x, y, z are the position coordinates of the intersection points of the isosurface and each edge of the voxel.

[0056] Using the intersection points of the isosurface and each edge of the voxel, P(x,y,z), and the triangular facets formed by the first vertex P1(x1,y1,z1) and the second vertex P2(x2,y2,z2) on the edge where the voxel intersects with the isosurface as isosurfaces, the isosurfaces are combined to form an isosurface triangular network.

[0057] To render the equivalent triangular network with better 3D effects, a lighting model needs to be selected based on the normal vectors of each triangular facet to render the equivalent triangular network and obtain the surface model of the root canal; wherein, the formula for calculating the normal vectors of each triangular facet is:

[0058]

[0059] In the formula, G x G y G z M represents the gradient at the vertex of the isosurface along the X, Y, and Z axes, respectively. (x+a,y,z) M (x-a,y,z) M represents the grayscale value of two points (x+a,y,z) and (xa,y,z) on the edge where the isosurfaces intersect. (x,y+b,z) M (x,y-b,z) M represents the grayscale values ​​of (x,y+b,z) and (x,yb,z). (x,y,z+c) M (x,y,z-c) V represents the grayscale values ​​of (x,y,z+x) and (x,y,zc), where a, b, and c represent the distances between the vertices of the two isosurfaces along the X, Y, and Z axes, respectively; x V y V z The normal vectors G of the triangular facets along the X, Y, and Z axes, respectively. x1 G y1 G z1 Let G be the gradient components of point P1(x1,y1,z1) along the X, Y, and Z axes, respectively. x2 G y2 Cz2 These are the gradient components of point P2(x2,y2,z2) along the X, Y, and Z axes, respectively.

[0060] Furthermore, based on the regular point cloud data in the aforementioned 3D model, a point-surface method is used to reconstruct the surface and generate a root canal solid model, including:

[0061] The n-order triangular Bezier surface is controlled by (n+1)(n+2) / 2 vertices p. i,j,k Definition of a triangular array consisting of (i,j,k≥0,i+j+k=n):

[0062]

[0063] In the formula, P(u,v,w) is any point on the Bezier surface of the triangular domain, u,v,w∈[0,1] represents the coordinates of the centroid inside the triangle, i,j,k represent the sequence numbers of any three different control vertices, and B i,j,k n (u,v,w) denotes the nth-order Bezier basis functions.

[0064] The coordinates of the centroid of the triangle can be calculated using the following formula:

[0065]

[0066] In the formula, area represents the area of ​​the triangle;

[0067] The nth-order control vertex P i,j,k Convert to an (n-1)th order control vertex Q i,j,k =uP i+1,j,k +vP i,j+1,k +wP i,j,k+1 By recursively applying this process, the remaining points are eventually points on the Bezier surface of the triangular domain, thus reconstructing the root canal solid model.

[0068] Furthermore, the root canal physical model is aligned in three-dimensional space, and the morphological characteristics of the root canal are measured using an approximate centerline fitted to the root canal, including:

[0069] x i =(x i ,y i ,z i (i = 1, 2, ..., n) can be represented as matrix X:

[0070]

[0071] In the formula, (x i y i , z i() represents the coordinates of the i-th point cloud, and n represents the amount of data in the point cloud;

[0072] The matrix can be obtained by decentralizing the point cloud data of a single tooth using the following formula.

[0073]

[0074] in,

[0075] Calculated using the following formula Covariance matrix:

[0076]

[0077] In the formula, V represents The covariance matrix, S represents the point cloud stretching matrix, satisfying S = S T R represents the point cloud rotation matrix, L represents the covariance in the new coordinate system after the point cloud is dimensionality reduced and transformed, and D represents the covariance matrix of the white data.

[0078] If the eigenvalues ​​of V are λ1, λ2, ..., λ n The corresponding feature vectors are p1, p2, ..., p n Then we have:

[0079] V·p i =λ i ·V (17)

[0080] It can be seen that the eigenvectors correspond to the rotation matrix R, representing the coordinate axis direction of each component; the eigenvalues ​​correspond to the stretching matrix S, representing the variance of the data in the corresponding component's coordinate axis direction.

[0081] Arrange the eigenvectors into a matrix from top to bottom according to the magnitude of their corresponding eigenvalues, and take the first m principal components to form a matrix W. m And determine the Z-axis vector of the principal components of the original point cloud after it has been reduced to m dimensions as follows:

[0082]

[0083] Let the Z-axis vector of the principal components obtained according to equation (18) be k = (z x ,z y ,z z ), with i = (1, -z) x / z y The X-axis vector is used as the X-axis vector, and the Y-axis vector is obtained by cross product of the X-axis and Z-axis vectors, which determines the principal component directions of the three-dimensional point cloud data of teeth.

[0084] The root canal length, orifice orientation, root canal curvature, and pulp depth are measured using a centerline approximation method in a coordinate system aligned with the principal component directions.

[0085] According to one aspect of the embodiments of this application, a device for extracting root canal features based on CBCT data is provided, comprising:

[0086] The 3D reconstruction module is configured to filter and denoise CBCT images and segment the pulp, and to perform surface rendering and 3D reconstruction of the root canals to generate a 3D model.

[0087] The surface reconstruction module is configured to use the point-surface method to reconstruct the surface based on the regular point cloud data in the three-dimensional model, and generate the root canal solid model.

[0088] The parameter estimation module is configured to perform coordinate alignment processing on the root canal entity model in three-dimensional space and measure the root canal morphological features by fitting the approximate centerline of the root canal; wherein the root canal morphological features include at least one of root canal length, root canal orifice direction, root canal curvature, and pulp depth.

[0089] According to one aspect of the embodiments of this application, an electronic device is provided, including: a controller; and a memory for storing one or more programs, which, when executed by the controller, cause the controller to implement the root canal feature morphology extraction method based on CBCT data described above.

[0090] According to one aspect of the embodiments of this application, a computer-readable storage medium is also provided, on which computer-readable instructions are stored, which, when executed by a computer processor, cause the computer to perform the above-described method for extracting root canal features based on CBCT data.

[0091] According to one aspect of the embodiments of this application, a computer program product or computer program is also provided, which includes computer instructions stored in a computer-readable storage medium. A processor of a computer device reads the computer instructions from the computer-readable storage medium and executes the computer instructions, causing the computer device to perform the aforementioned method for extracting root canal features based on CBCT data.

[0092] The technical solutions provided in the embodiments of this application have at least the following advantages:

[0093] This application utilizes CBCT tomographic images obtained from scanning a patient's teeth to rapidly reconstruct a corresponding three-dimensional model, containing a wealth of morphological and locational information about the tooth structure, pulp, and root canals. Since access cavity preparation is required before root canal preparation, obtaining data such as pulp depth and root canal orifice size from the three-dimensional model of the tooth can help dentists determine the access cavity preparation plan, which has practical significance and application value for current root canal preparation surgery.

[0094] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit this application. Attached Figure Description

[0095] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application. It is obvious that the drawings described below are merely some embodiments of this application, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort. In the drawings:

[0096] Figure 1 This is an overall flowchart illustrating an exemplary embodiment of the present application of a method for extracting root canal features based on CBCT data;

[0097] Figure 2 This is a flowchart illustrating the root canal surface reconstruction process in an exemplary embodiment of this application;

[0098] Figure 3 This is a schematic diagram of an edge model shown in an exemplary embodiment of this application;

[0099] Figure 4 This is a comparison diagram of filtering effects shown in an exemplary embodiment of this application, wherein (a) means filtering; and (b) median filtering with edge preservation.

[0100] Figure 5 This is an exemplary embodiment of the present application illustrating the overall flowchart of pulp segmentation;

[0101] Figure 6 This is a schematic diagram of a three-dimensional segmented space shown in an exemplary embodiment of this application;

[0102] Figure 7 This is a schematic diagram of a voxel model shown in an exemplary embodiment of this application;

[0103] Figure 8 This is a schematic diagram of the centroid coordinates of a triangle shown in an exemplary embodiment of this application;

[0104] Figure 9This is a schematic diagram of a third-order triangular Bezier surface control network illustrated in an exemplary embodiment of this application;

[0105] Figure 10 This is a schematic diagram illustrating a recursive calculation process in an exemplary embodiment of this application;

[0106] Figure 11 This is a flowchart illustrating the measurement of root canal morphology features according to an exemplary embodiment of this application;

[0107] Figure 12 This is a schematic diagram of a cross-sectional curve shown in an exemplary embodiment of this application;

[0108] Figure 13 This is a schematic diagram of a fitted projection curve shown in an exemplary embodiment of this application;

[0109] Figure 14 This is a schematic diagram of the root canal centerline shown in an exemplary embodiment of this application;

[0110] Figure 15 This is a structural diagram of a root canal feature morphology extraction device based on CBCT data, as illustrated in an exemplary embodiment of this application. Detailed Implementation

[0111] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numbers in different drawings denote the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this application as detailed in the appended claims.

[0112] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0113] The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily need to be performed in the described order. For example, some operations / steps can be broken down, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0114] In this application, "multiple" refers to two or more. "And / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A alone, A and B simultaneously, or B alone. The character " / " generally indicates that the preceding and following related objects have an "or" relationship.

[0115] Please see Figure 1 , Figure 1 This is an overall flowchart illustrating an exemplary embodiment of the present application of a method for extracting root canal features based on CBCT data. One aspect of this application provides a method for extracting root canal features based on CBCT data, such as... Figure 1 As shown, the method includes steps S100 to S300. Details are as follows.

[0116] S100 performs filtering and noise reduction on CBCT images and pulp segmentation, and performs surface rendering and 3D reconstruction of root canals to generate 3D models.

[0117] In this embodiment, to reduce interference from irrelevant tissues such as gums and muscles in the image, a preprocessing operation of filtering and denoising is required on the CBCT data before tooth surface reconstruction. Then, the preprocessed CBCT image is used for surface rendering and reconstruction to generate three-dimensional point cloud data of the target tooth.

[0118] In some embodiments, please refer to Figure 2 This is a flowchart of root canal surface reconstruction. Step S100 can be implemented in practice through the following steps S110 and S120.

[0119] S110: Based on the acquired CBCT images, the neighborhood edge model is used to improve the denoising effect of median filtering, and the optimal threshold for pulp segmentation is calculated based on the intercept histogram of the inverse cross-entropy, thus obtaining the preprocessed CBCT sequence images.

[0120] In this embodiment, step S110 can be implemented in two steps. One step is CBCT image preprocessing, which uses a neighborhood edge model to improve the denoising effect of median filtering based on the acquired CBCT image. The other step is pulp segmentation, which calculates the optimal threshold for pulp segmentation based on the intercept histogram of the inverse cross-entropy to obtain the preprocessed CBCT sequence image.

[0121] Specifically, during CBCT image preprocessing, the image edges are first determined and preserved. For example, within a 3×3 neighborhood, the center pixel M of that neighborhood is selected for preservation. c Establish 8 edge models such as Figure 3 As shown, calculate the pixel difference between this pixel and its eight nearest neighbors (top, bottom, left, right, and 45° angle):

[0122] |M i -M c | = d(i = 1, 2, …, 7, 8) (1)

[0123] In the formula, M i represents the i-th neighboring point, where i is a positive integer not greater than 8. The 8 neighboring points are respectively located above, below, left, right of the central pixel point and in 4 45° angular directions. M c represents the central pixel point, and d represents the pixel difference between the central pixel point and the i-th neighboring point.

[0124] Set the threshold size as T1, and record the number of neighboring points with gray levels similar to that of the central pixel point as n. When the pixel difference d < T1, it is considered that the neighboring point has a gray level similar to that of the central pixel point, and let n = n + 1. After traversing all neighboring points, if n min ≤ n ≤ n max That is, the number of neighboring points conforms to the edge model established above, then it is judged that the central pixel point is an edge point and should be directly retained. n min is the minimum value of the number of neighboring points, and n max is the maximum value of the number of neighboring points.

[0125] For non-edge points, replace the non-edge points with the median value of all pixel points in the neighborhood. Calculate the median value of all pixel points through the following formula:

[0126] F med = medM (x,y) = med[M (x+m,y+n) ; (m, n) ∈ A] (2)

[0127] In the formula, F med is the median value of all pixel points, A is the neighborhood around the central pixel point (x, y), m and n are the length and width of the neighborhood window, M (x,y) is the pixel value of the pixel point (x, y), med is the median operation, and M (x+m,y+n) is the pixel value of the pixel point (x + m, y + n);

[0128] Execute the first step: If d1 = F med - F min > 0 and d2 = F max - F med > 0, then execute the third step. If not satisfied, execute the second step;

[0129] The second step of median filtering includes: expanding the neighborhood area by one circle. If the neighborhood area meets the second condition, then sort all pixel values in the neighborhood again and update the minimum value F min and the median value F medMaximum value F max Then repeat step one; if the neighborhood area does not meet the second condition, then do not replace the grayscale value of the center pixel; wherein, the second condition is that the neighborhood area is less than the maximum neighborhood threshold S. max ;

[0130] The third step includes: if g1 = M c -F min >0 and g2=F max -M c If the grayscale value is greater than 0, then the grayscale value of the center pixel is not replaced; otherwise, the grayscale value of the center pixel is replaced with the median value F. med .

[0131] The effect of using mean filtering and edge-preserving median filtering is, for example... Figure 4 As shown in (a) and (b).

[0132] The overall flowchart of pulp segmentation is as follows: Figure 5 As shown, during pulp segmentation, for a specific layer of CBCT image, gamma transform is first used to enhance image contrast. Then, the neighborhood average image and gradient image are calculated separately. These two images are combined with a grayscale image and added pixel by pixel to form a hybrid image. Based on the hybrid image, its histogram can be used to find the optimal threshold. Each pixel value in the hybrid image is then compared with the optimal threshold and converted into a binary image. Finally, all unwanted connected tissues at the image boundaries are removed, retaining only the pulp region.

[0133] A three-dimensional space is generated by combining grayscale information, neighborhood average grayscale information, and gradient compensation, such as... Figure 6 As shown, the entire space is divided into the target region O and the background region B by a plane α.

[0134] The gray value of point (x,y) and the average gray value of its neighborhood are denoted as g(x,y) and f(x,y), respectively, and the gradient compensation is denoted as G(x,y). T is the segmentation threshold of the two-dimensional plane determined by g(x,y) and f(x,y), and G... T This refers to gradient compensation at a threshold of T. The gradient compensation for each threshold is a corresponding fixed value. Therefore, the optimal segmentation plane α is determined only by the threshold T, and the value of T is determined by the line k in the two-dimensional plane. Each line k intersects the dashed line v, so the number of such lines is 2L-1, meaning the threshold range is T∈1,2,…,2L-1.

[0135] Calculate the gradient values ​​of each point in the horizontal and vertical directions. Given a point (x, y) in an image with height N and width M, where x∈1,2,…,M-1,M and y∈1,2,…,N-1,N, the gradient values ​​in the horizontal and vertical directions are denoted as G.H (x,y) and G V (x,y). The final gradient G(x,y) is the combination of the gradients in the two directions.

[0136]

[0137] In the formula, G H (x,y) and G V (x,y) represents the gradient values ​​of a point (x,y) in the horizontal and vertical directions in a given image, where x represents the x-coordinate of the point in the given image, y represents the y-coordinate of the point in the given image, G(x,y) represents the final gradient, g(x,y+1) represents the gray value of point (x,y+1), g(x,y-1) represents the gray value of point (x,y-1), g(x+1,y) represents the gray value of point (x+1,y), and g(x-1,y) represents the gray value of point (x+1,y).

[0138] Since the threshold is the same, and both the grayscale axis and the neighborhood average grayscale axis are T, the grayscale image, the neighborhood average image, and the gradient image are added together to generate a mixed image, represented as:

[0139] F(x,y)=f(x,y)+g(x,y)-G(x,y) (5)

[0140] In the formula, F(x,y) represents the mixed image.

[0141] For blended images, the optimal threshold can be found by calculating the corresponding histogram. In the histogram, the blended image with a threshold of T is divided into the target portion Ω. O ∈{(x,y)|F(x,y)=0,1,…,T} and background part Ω B ∈{(x,y)|F(x,y)=T+1,T+2,…,2T-2}, pixel values ​​less than T belong to Ω. O Pixel values ​​greater than T belong to Ω. B .

[0142] Then, the inverse cross-entropy between the target region and the background region is calculated using the following formula:

[0143]

[0144] In the formula, E(O,B) represents the reciprocal cross-entropy; P(k) represents the prior probability that the average gray value of the neighborhood of the mixed image is k. This represents the average gray level of the target region. This represents the average grayscale value of the background area;

[0145] The inverse cross-entropy E(O,B) reflects the deviation before and after image segmentation. Therefore, when E(O,B) reaches its minimum value, we have:

[0146]

[0147] In the formula, This is the optimal threshold for pulp segmentation, and arg min is the operation to find the index of the smallest element in the array.

[0148] S120 uses surface rendering to draw the surface model of the root canal from the preprocessed CBCT sequence image.

[0149] The purpose of step S120 is to achieve three-dimensional reconstruction, thereby obtaining the surface model of the root canal, i.e., the three-dimensional model of the tooth.

[0150] For example, if uniform sampling is performed in the x, y, z directions of the patient's upper and lower jaw regions in three-dimensional space, with sampling intervals of Δx, Δy, and Δz, then the volume data can be represented by a ternary function p(i,j,k). The cubic region formed by eight adjacent sampling points constitutes one voxel, such as... Figure 7 As shown.

[0151] The essence of the McLeod algorithm is to extract isosurfaces for each voxel. An isosurface is a surface composed of points with the same properties in 3D space. Reconstructing different tissues requires setting different isosurface thresholds; for teeth, the threshold is set to around 750. Comparing the vertices of each voxel with the set isosurface threshold reveals that these vertices have two different states. When the states of two vertices on one side of a voxel differ, that side will inevitably intersect with an isosurface. Due to the complementary and rotational symmetries of space, the configurations of voxels intersecting with isosurfaces can be simplified to 15.

[0152] Traverse the entire volume data, find the voxels containing isosurfaces, and use linear interpolation to determine the intersection points P(x,y,z) between the isosurfaces and the edges of the voxels:

[0153]

[0154] In the formula, T is the isosurface threshold, M1 and M2 are the gray values ​​of the first and second vertices on the edge where the voxel intersects the isosurface, x1, y1, z1 are the position coordinates of the first vertex on the edge where the voxel intersects the isosurface, x2, y2, z2 are the position coordinates of the second vertex on the edge where the voxel intersects the isosurface, and x, y, z are the position coordinates of the intersection points of the isosurface and each edge of the voxel.

[0155] Where T is the isosurface threshold, M1 and M2 are the gray values ​​of two vertices on the edge where the voxel intersects the isosurface, and P1(x1,y1,z1) and P2(x2,y2,z2) are the positions of these two vertices. The triangular patches formed by these intersection points are the extracted isosurfaces. To render the isosurface triangulation network with better 3D effects, a suitable lighting model needs to be selected based on the normal vectors of each triangular patch. First, calculate the gradients in the 3D directions at the vertices of the isosurface:

[0156]

[0157] In the formula, G x G y G z M represents the gradient at the vertex of the isosurface along the X, Y, and Z axes, respectively. (x+a,y,z) M (x-a,y,z) M represents the gray values ​​of two vertices (x+z,y,z) and (xa,y,z) on the edge where the isosurfaces intersect. (x,y+b,z) M (x,y-b,z) M represents the grayscale values ​​of (x,y+b,z) and (x,yb,z). (x,y,z+c) M (x,y,z-c) V represents the grayscale values ​​of (x,y,z+c) and (x,y,zc), where a, b, and c represent the distances between the vertices of the two isosurfaces along the X, Y, and Z axes, respectively; x V y V z The normal vectors G of the triangular facets along the X, Y, and Z axes, respectively. x1 G y1 G z1 Let G be the gradient components of point P1(x1,y1,z1) along the X, Y, and Z axes, respectively. x2 G y2 G z2 These are the gradient components of point P2(x2,y2,z2) along the X, Y, and Z axes, respectively.

[0158] S200, based on regular point cloud data in a 3D model, uses the point-surface method to reconstruct the surface and generate a root canal solid model.

[0159] In this embodiment, a method based on triangular domain Bezier surface reconstruction is used to recursively refine a large number of triangular facets on the root canal surface using a mesh, thereby realizing the surface reconstruction of the root canal model.

[0160] For example, a triangular Bezier surface is a control vertex connected by elementary functions, that is, an n-order triangular Bezier surface consists of (n+1)(n+2) / 2 control vertices p. i,j,kDefinition of a triangular array consisting of (i,j,k≥0,i+j+k=n):

[0161]

[0162] In the formula, P(u,v,w) is any point on the Bezier surface of the triangular domain, u,v,w∈[0,1] represents the coordinates of the centroid inside the triangle, i,j,k represent any three distinct control vertices, and B i,j,k n (u,v,w) denotes the nth-order Bezier basis functions.

[0163] Equation (11) is the Bezier basis function, and the coordinates of the centroid of the triangle in the equation are as follows: Figure 8 As shown.

[0164] The coordinates of the centroid of the triangle can be calculated using equation (13):

[0165]

[0166] In the formula, area represents the area of ​​the triangle.

[0167] In this embodiment, the parameters u = v = w = 1 / 3 can be set.

[0168] In the calculation of the Bezier surface in the third-order triangular domain, 10 control vertices, denoted as p, will be used. i,j,k (i,j,k≥0,i+j+k=10), these 10 control vertices are connected by straight lines following the index order, generating a surface control network composed of triangles. The control vertices of the network correspond one-to-one with the nodes of the triangular regions, as shown below. Figure 9 As shown.

[0169] nth-order control vertex P i,j,k It will be converted into an n-1 order control vertex Q i,j,k =uP i+1,j,k +vP i,j+1,k +wP i,j,j+1 This process is repeated recursively until the last remaining point is a point on the Bezier surface of the triangular region. The recursive calculation process is as follows: Figure 10 As shown.

[0170] S300 aligns the root canal physical model in three-dimensional space and uses the approximate centerline of the root canal to measure the morphological characteristics of the root canal; wherein, the morphological characteristics of the root canal include at least one of root canal length, root canal orifice direction, root canal curvature and pulp depth.

[0171] In this embodiment, the purpose of step S300 is to measure the morphological characteristics of the tooth root canal. Specifically, principal component analysis is used to determine the direction of the principal components of the tooth and use it as the Z-axis of the reference coordinate system. Then, a centerline approximation method is used to detect the root canal length and the direction of the root canal orifice.

[0172] For example, the procedure for measuring root canal morphology is as follows: Figure 11 As shown, it includes the following steps S301 and S302.

[0173] Step S301, coordinate alignment.

[0174] x i =(x i ,y i ,z i (i = 1, 2, ..., n) can be represented as matrix X:

[0175]

[0176] The matrix can be obtained by decentralizing the point cloud data of a single tooth using the following formula. :

[0177]

[0178] in,

[0179] Calculated using the following formula Covariance matrix:

[0180]

[0181] In the formula, V represents The covariance matrix, S represents the point cloud stretching matrix, satisfying S = S T R represents the point cloud rotation matrix, L represents the covariance in the new coordinate system after the point cloud is dimensionality reduced and transformed, and D represents the covariance matrix of the white data.

[0182] If the eigenvalues ​​of V are λ1, λ2, ..., λ n The corresponding feature vectors are p1, p2, ..., p n Then we have:

[0183] V·p i =λ i ·V (17)

[0184] It can be seen that the eigenvectors correspond to the rotation matrix R, representing the coordinate axis direction of each component; the eigenvalues ​​correspond to the stretching matrix S, representing the variance of the data in the corresponding component's coordinate axis direction.

[0185] Arrange the eigenvectors into a matrix from top to bottom according to the magnitude of their corresponding eigenvalues, and take the first m principal components to form a matrix W. m And determine the data after reducing the original point cloud to m dimensions as follows:

[0186]

[0187] In the formula, F represents the Z-axis vector of the principal component;

[0188] Let the Z-axis vector of the principal components obtained according to equation (18) be k = (z x ,z y ,z z ), with i = (1, -z) x / z y The X-axis vector is used as the X-axis vector, and the Y-axis vector is obtained by cross product of the X-axis and Z-axis vectors, which determines the principal component directions of the three-dimensional point cloud data of the teeth.

[0189] Step S302, parameter calculation.

[0190] The centerline inside a tooth root canal is a spatial curve. The projections of this centerline onto the buccal-lingual and mesiodistal directions of the root canal are two planar curves. Therefore, the constructed root canal model is first sectioned in both the buccal-lingual and mesiodistal planes, as shown below. Figure 12 As shown, two sets of cross-sectional curves of the root canal profile can be obtained.

[0191] By performing spline fitting on the transverse midpoints of these two sets of cross-sectional curves, the projection curves of the root canal centerline in the buccal-lingual and mesiodistal directions can be obtained. Here, the XOZ plane represents the buccal-lingual direction, and the YOZ plane represents the mesiodistal direction. Figure 13 As shown.

[0192] Select two surfaces containing these two projection curves and perpendicular to their corresponding buccal-lingual plane or mesiodistal plane, and let them intersect. Figure 14 As shown, the spatial curve formed after the intersection is the centerline of the tooth root canal.

[0193] Therefore, the length of the entire root canal and the opening direction of the root canal orifice can be calculated based on the approximate root canal center curve. The angle between the tangent direction of the center line at a point on the root canal orifice and the direction of the principal components of the tooth is the opening direction of the root canal orifice.

[0194] Please see Figure 15 , Figure 15 This is a structural diagram of a root canal feature morphology extraction device based on CBCT data, illustrating an exemplary embodiment of this application. Another aspect of this application provides a root canal feature morphology extraction device based on CBCT data, the device comprising:

[0195] The three-dimensional reconstruction module 1501 is configured to perform filtering denoising and pulp segmentation processing on the CBCT image, and perform surface rendering three-dimensional reconstruction on the tooth root canal to generate a three-dimensional model;

[0196] The surface reconstruction module 1502 is configured to perform surface reconstruction by using the point-plane method based on the regular point cloud data in the three-dimensional model to generate a root canal solid model;

[0197] The parameter estimation module 1503 is configured to perform coordinate alignment processing on the root canal solid model in the three-dimensional space, and measure the morphological characteristics of the tooth root canal by using the approximate center line fitting the root canal; wherein, the morphological characteristics of the tooth root canal include at least one of the root canal length, the root canal opening direction, the root canal curvature, and the pulp depth.

[0198] In some embodiments, the three-dimensional reconstruction module is further configured to:

[0199] Based on the acquired CBCT image, use the neighborhood edge model to improve the denoising effect of the median filter, and calculate the optimal threshold for segmenting the pulp based on the intercept histogram of the reciprocal cross entropy to obtain the preprocessed CBCT sequence image;

[0200] Draw the surface model of the root canal by surface rendering the preprocessed CBCT sequence image.

[0201] In some embodiments, the three-dimensional reconstruction module is further configured to:

[0202] In a neighborhood, establish 8 edge models for the central pixel point of the neighborhood, and calculate the pixel difference between the central pixel point and the neighboring points through the following formula: <o:p>< / o:p> <o:p>< / o:p><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:mrow><m:m​​​​​​​​​​​​​​​​​​​​The minimum number of neighboring points, n max If the number of neighboring points is the maximum value, then the number of neighboring points is determined to conform to the established edge model, and the center pixel is determined to be an edge point and directly retained;

[0206] For non-edge points, the non-edge point is replaced with the median of all pixels in its neighborhood, and the median of all pixels is calculated using the following formula:

[0207] F med =medM (x,y) =med[M (x+m,y+n) ;(m,n)∈A] (2)

[0208] In the formula, F med Let A be the median of all pixels, A be the neighborhood surrounding the center pixel (x, y), m and n be the length and width of the neighborhood window, and M be the median of all pixels. (x,y) Here, M represents the pixel value at (x, y), med is the median value, and M is the median value. (x+m,y+n) The pixel value at pixel (x+m, y+n);

[0209] Sort all pixel values ​​within the neighborhood and obtain the minimum value F of all pixels. min Median F med Maximum value F max ;

[0210] Execute the first step: If d1 = F med -F min >0 and d2=F max -F med If the value is greater than 0, proceed to step three; otherwise, proceed to step two.

[0211] The second step includes: expanding the neighborhood area by one ring; if the neighborhood area satisfies the second condition, then sorting all pixel values ​​within the neighborhood again and updating the minimum value F of all pixels. min Median F med Maximum value F max Then repeat step one; if the neighborhood area does not meet the second condition, then do not replace the grayscale value of the center pixel; wherein, the second condition is that the neighborhood area is less than the maximum neighborhood threshold S. max ;

[0212] The third step includes: if g1 = M c -F min >0 and g2=F max -M c If the grayscale value is greater than 0, then the grayscale value of the center pixel is not replaced; otherwise, the grayscale value of the center pixel is replaced with the median value f. med .

[0213] In some embodiments, the three-dimensional reconstruction module is further configured to:

[0214] Based on the denoised CBCT images, gamma transform is used to enhance image contrast, and the neighborhood average image and gradient image are calculated respectively.

[0215] The neighborhood average image and gradient image are combined with the grayscale image and added pixel by pixel to form a hybrid image;

[0216] The optimal threshold for pulp segmentation is calculated using the intercept histogram of the mixed image;

[0217] Each pixel value in the blended image is compared with the optimal threshold and converted into a binary image;

[0218] Remove all unnecessary connected tissues on the image boundaries of the binary image, retaining only the pulp region, to obtain a preprocessed CBCT sequence image;

[0219] The step of calculating the optimal threshold for pulp segmentation using the intercept histogram of the mixed image includes:

[0220] A three-dimensional space is generated by grayscale information, neighborhood average grayscale information and gradient compensation. The three-dimensional space is divided into target region O and background region B by a plane.

[0221] Based on the aforementioned three-dimensional space, the final gradient is calculated using the following formula:

[0222]

[0223] In the formula, G H (x,y) and G V (x,y) represents the gradient values ​​of a point (x,y) in the horizontal and vertical directions in a given image, where x represents the x-coordinate of the point in the given image, y represents the y-coordinate of the point in the given image, G(x,y) represents the final gradient, g(x,y+1) represents the gray value of point (x,y+1), g(x,y-1) represents the gray value of point (x,y-1), g(x+1,y) represents the gray value of point (x+1,y), and g(x-1,y) represents the gray value of point (x+1,y).

[0224] The hybrid image is represented as:

[0225] F(x,y)=f(x,y)+g(x,y)-G(x,y) (5)

[0226] In the formula, F(x,y) represents the mixed image, and g(x,y) and f(x,y) represent the gray value of point (x,y) and the average gray value of the neighborhood of point (x,y).

[0227] Based on the segmentation threshold T, the hybrid image is divided into the target portion Ω. O ∈{(x,y)|F(x,y)=0,1,…,T} and background part Ω B ∈{(x,y)|F(x,y)=T+1,T+2,…,2T-2}, pixel values ​​less than the segmentation threshold T belong to the target part, and pixel values ​​greater than the segmentation threshold T belong to the background part;

[0228] The inverse cross-entropy between the target region and the background region is calculated using the following formula:

[0229]

[0230] In the formula, E(O,B) represents the reciprocal cross-entropy; P(k) represents the prior probability that the average gray value of the neighborhood of the mixed image is k. This represents the average gray level of the target region. This represents the average grayscale value of the background area;

[0231] Based on the calculated inverse cross-entropy E(O,B), the optimal threshold for pulp segmentation is determined using the following formula:

[0232]

[0233] In the formula, The optimal threshold for pulp segmentation is arg min.

[0234] In some embodiments, the three-dimensional reconstruction module is further configured to:

[0235] In three-dimensional space, uniform sampling is performed on the upper and lower jaw regions in the x, y, and z directions with sampling intervals of Δx, Δy, and Δz. The volume data is represented by a ternary function; where a cubic region consisting of eight adjacent sampling points is a voxel.

[0236] For each voxel, an isosurface is extracted; wherein, the isosurface is a surface composed of points with the same properties in three-dimensional space;

[0237] Traverse the entire volume data to find the voxels containing isosurfaces, and determine the intersection points P(x,y,z) of the isosurfaces and the edges of the voxels using the following formula:

[0238]

[0239] In the formula, T is the isosurface threshold, M1 and M2 are the gray values ​​of the first and second vertices on the edge where the voxel intersects with the isosurface, x1, y1, z1 are the position coordinates of the first vertex on the edge where the voxel intersects with the isosurface, x2, y2, z2 are the position coordinates of the second vertex on the edge where the voxel intersects with the isosurface, and x, y, z are the position coordinates of the intersection points of the isosurface and each edge of the voxel.

[0240] Using the intersection points of the isosurface and each edge of the voxel, P(x,y,z), and the triangular facets formed by the first vertex P1(x1,y1,z1) and the second vertex P2(x2,y2,z2) on the edge where the voxel intersects with the isosurface as isosurfaces, the isosurfaces are combined to form an isosurface triangular network.

[0241] To render the equivalent triangular network with better 3D effects, a lighting model needs to be selected based on the normal vectors of each triangular facet to render the equivalent triangular network and obtain the surface model of the root canal; wherein, the formula for calculating the normal vectors of each triangular facet is:

[0242]

[0243] In the formula, G x G y G z M represents the gradient at the vertex of the isosurface along the X, Y, and Z axes, respectively. (x+a,y,z) M (x-a,y,z) M represents the grayscale value of two points (x+a,y,z) and (xa,y,z) on the edge where the isosurfaces intersect. (x,y+b,z) M (x,y-b,z) M represents the grayscale values ​​of (x,y+b,z) and (x,yb,z). (x,y,z+c) M (x,y,z-c) V represents the grayscale values ​​of (x,y,z+c) and (x,y,zc), where a, b, and c represent the distances between the vertices of the two isosurfaces along the X, Y, and Z axes, respectively; x V y V z The normal vectors G of the triangular facets along the X, Y, and Z axes, respectively. x1 G y1 G z1 Let G be the gradient components of point P1(x1,y1,z1) along the X, Y, and Z axes, respectively. x2 G y2 G z2 These are the gradient components of point P2(x2,y2,z2) along the X, Y, and Z axes, respectively.

[0244] In some embodiments, the surface reconstruction module is further configured to:

[0245] The n-order triangular Bezier surface is controlled by (n+1)(n+2) / 2 vertices p. i,j,k Definition of a triangular array consisting of (i,j,k≥0,i+j+k=n):

[0246]

[0247] In the formula, P(u,v,w) is any point on the Bezier surface of the triangular domain, u,v,w∈[0,1] represents the coordinates of the centroid inside the triangle, i,j,k represent any three distinct control vertices, and B i,j,k n (u,v,w) denotes the nth-order Bezier basis functions.

[0248] The coordinates of the centroid of the triangle can be calculated using the following formula:

[0249]

[0250] In the formula, area represents the area of ​​the triangle;

[0251] The nth-order control vertex P i,j,k Convert to an (n-1)th order control vertex Q i,j,k =uP i+1,j,k +vP i,j+1,k +wP i,j,k+1 By recursively applying this process, the remaining points are eventually points on the Bezier surface of the triangular domain, thus reconstructing the root canal solid model.

[0252] In some embodiments, the parameter estimation module is further configured to:

[0253] x i =(x i ,y i ,z i (i = 1, 2, ..., n) can be represented as matrix X:

[0254]

[0255] The matrix can be obtained by decentralizing the point cloud data of a single tooth using the following formula.

[0256]

[0257] in,

[0258] Calculated using the following formula Covariance matrix:

[0259]

[0260] In the formula, V represents The covariance matrix, S represents the point cloud stretching matrix, satisfying S = S T R represents the point cloud rotation matrix, L represents the covariance in the new coordinate system after the point cloud is dimensionality reduced and transformed, and D represents the covariance matrix of the white data.

[0261] If the eigenvalues ​​of V are λ1, λ2, ..., λ n The corresponding feature vectors are p1, p2, ..., p n Then we have:

[0262] V·p i =λ i ·V (17)

[0263] It can be seen that the eigenvectors correspond to the rotation matrix R, representing the coordinate axis direction of each component; the eigenvalues ​​correspond to the stretching matrix S, representing the variance of the data in the corresponding component's coordinate axis direction.

[0264] Arrange the eigenvectors into a matrix from top to bottom according to the magnitude of their corresponding eigenvalues, and take the first m principal components to form a matrix W. m And determine the Z-axis vector of the principal components of the data after the original point cloud is reduced to m dimensions as follows:

[0265]

[0266] Let the Z-axis vector of the principal components obtained according to equation (18) be k = (z x ,z y ,z z ), with i = (1, -z) x / z y The X-axis vector is used as the X-axis vector, and the Y-axis vector is obtained by cross product of the X-axis and Z-axis vectors, which determines the principal component directions of the three-dimensional point cloud data of the teeth.

[0267] The root canal length, orifice orientation, root canal curvature, and pulp depth are measured using a centerline approximation method in a coordinate system aligned with the principal component directions.

[0268] It should be noted that the root canal feature morphology extraction device based on CBCT data provided in the above embodiments and the root canal feature morphology extraction method based on CBCT data provided in the aforementioned embodiments belong to the same concept. The specific way of performing the steps has been described in detail in the method embodiments, and will not be repeated here.

[0269] Another aspect of this application provides an electronic device, including: a controller; and a memory for storing one or more programs, which, when executed by the controller, perform the root canal feature morphology extraction method based on CBCT data in the various embodiments described above.

[0270] Specifically, according to embodiments of this application, the processes described above with reference to the flowcharts can be implemented as computer software programs. For example, embodiments of this application include a computer program product comprising a computer program carried on a computer-readable medium, the computer program including a computer program for performing the methods shown in the flowcharts. In such embodiments, the computer program can be downloaded and installed from a network via a communication component, and / or installed from a removable medium. When the computer program is executed by the central processing unit (CPU) 701, it performs various functions defined in the system of this application.

[0271] It should be noted that the computer-readable medium shown in the embodiments of this application can be a computer-readable signal medium or a computer-readable storage medium, or any combination of the two. A computer-readable storage medium can be, for example, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of a computer-readable storage medium may include, but are not limited to: an electrical connection having one or more wires, a portable computer disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM), flash memory, optical fiber, portable compact disc read-only memory (CD-ROM), optical storage device, magnetic storage device, or any suitable combination thereof. In this application, a computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device. In this application, a computer-readable signal medium can include a data signal propagated in baseband or as part of a carrier wave, carrying a computer-readable computer program. Such propagated data signals can take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. Computer-readable signal media can also be any computer-readable medium other than computer-readable storage media, which can send, propagate, or transmit a program for use by or in connection with an instruction execution system, apparatus, or device. The computer program contained on the computer-readable medium can be transmitted using any suitable medium, including but not limited to wireless, wired, etc., or any suitable combination thereof.

[0272] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of this application. Each block in a flowchart or block diagram may represent a module, segment, or portion of code, which contains one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those indicated in the drawings. For example, two consecutively indicated blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram or flowchart, and combinations of blocks in a block diagram or flowchart, may be implemented using a dedicated hardware-based system that performs the specified function or operation, or using a combination of dedicated hardware and computer instructions.

[0273] The modules / units described in the embodiments of this application can be implemented in software or hardware, and the described units can also be located in a processor. The names of these modules / units do not necessarily constitute a limitation on the module / unit itself.

[0274] Another aspect of this application provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the root canal feature morphology extraction method based on CBCT data as described above. This computer-readable storage medium may be included in the electronic device described in the above embodiments, or it may exist independently without being assembled into the electronic device.

[0275] Another aspect of this application provides a computer program product or computer program that includes computer instructions stored in a computer-readable storage medium. A processor of a computer device reads the computer instructions from the computer-readable storage medium and executes the computer instructions, causing the computer device to perform the root canal feature morphology extraction method based on CBCT data provided in the various embodiments described above.

[0276] According to one aspect of the embodiments of this application, a computer system is also provided, including a Central Processing Unit (CPU), which can perform various appropriate actions and processes based on a program stored in read-only memory (ROM) or a program loaded from storage into random access memory (RAM), such as performing the methods described above. Various programs and data required for system operation are also stored in the RAM. The CPU, ROM, and RAM are interconnected via a bus. Input / output (I / O) interfaces are also connected to the bus.

[0277] The following components are connected to the I / O interface: input sections including keyboards, mice, etc.; output sections including cathode ray tubes (CRTs), liquid crystal displays (LCDs), and speakers; storage sections including hard drives; and communication sections including network interface cards such as LAN (Local Area Network) cards and modems. The communication sections perform communication processing via networks such as the Internet. Drives are also connected to the I / O interface as needed. Removable media, such as disks, optical discs, magneto-optical discs, semiconductor memories, etc., are installed on the drive as needed so that computer programs read from them can be installed into the storage section as required.

[0278] The above description is merely a preferred exemplary embodiment of this application and is not intended to limit the implementation of this application. Those skilled in the art can easily make corresponding modifications or alterations based on the main concept and spirit of this application. Therefore, the scope of protection of this application should be determined by the scope of protection claimed in the claims.

Claims

1. A method for extracting root canal features based on CBCT data, characterized in that, The method includes: The CBCT images were filtered and denoised, and the pulp was segmented. The root canals of the teeth were then rendered and reconstructed in three dimensions to generate a three-dimensional model. Based on the regular point cloud data in the three-dimensional model, the point-surface method is used to reconstruct the surface and generate the root canal solid model. The root canal solid model is aligned in three-dimensional space, and the morphological characteristics of the root canal are measured by fitting the approximate centerline of the root canal; wherein, the morphological characteristics of the root canal include at least one of root canal length, root canal orifice direction, root canal curvature and pulp depth. Based on the regular point cloud data in the aforementioned 3D model, a point-surface method is used to reconstruct the surface and generate a root canal solid model, including: The n-order triangular Bezier surface is controlled by (n+1)(n+2) / 2 vertices p. i,j,k Definition of a triangular array consisting of (i,j,k≥0,i+j+k=n): In the formula, P(u,v,w) is any point on the Bezier surface of the triangular domain, u,v,w∈[0,1] represents the coordinates of the centroid inside the triangle, i,j,k represent the sequence numbers of any three different control vertices, and B i,j,k n (u,v,w) denotes the nth-order Bezier basis functions; The coordinates of the centroid of the triangle can be calculated using the following formula: In the formula, area represents the area of ​​the triangle; The nth-order control vertex P i,j,k Convert to an (n-1)th order control vertex Q i,j,k =uP i+1,j,k +vP i,j+1,k +wP i,j,k+1 By recursively applying this process, the remaining points are eventually points on the Bezier surface of the triangular domain, thus reconstructing the root canal solid model.

2. The method for extracting root canal features based on CBCT data according to claim 1, characterized in that, The CBCT images were filtered and denoised, and the pulp was segmented. The root canals were then rendered and reconstructed in 3D to generate a 3D model, including: Based on the acquired CBCT images, the neighborhood edge model is used to improve the denoising effect of median filtering, and the optimal threshold for pulp segmentation is calculated based on the intercept histogram of the inverse cross-entropy, thus obtaining the preprocessed CBCT sequence images. The surface model of the root canal is drawn by rendering the preprocessed CBCT sequence images using a surface rendering method.

3. The method for extracting root canal features based on CBCT data according to claim 2, characterized in that, Based on the acquired CBCT images, a neighborhood edge model is used to improve the denoising effect of median filtering, including: Within a neighborhood, eight edge models are established for the center pixel of the neighborhood, and the pixel difference between the center pixel and its neighboring pixels is calculated using the following formula: |M i -M c |=d(i=1,2,…,7,8) (1) In the formula, M i Let M represent the i-th neighboring pixel, where i is a positive integer not greater than 8. The eight neighboring pixels are located above, below, left, right, and at four 45° angles from the center pixel. c d represents the center pixel, and d represents the pixel difference between the center pixel and the i-th neighboring pixel. Set the threshold size to T1, and record the number of adjacent points with gray levels similar to the central pixel as n. If the pixel difference d between the central pixel and the adjacent point is less than T1, then determine that the adjacent point has a gray level similar to the central pixel; let n = n + 1. When all adjacent points have been traversed, if n min ≤ n ≤ n max , n min is the minimum value of the number of adjacent points, and n max is the maximum value of the number of adjacent points, then determine that the number of adjacent points conforms to the established edge model, and determine that the central pixel is an edge point and directly retain it; For non-edge points, the non-edge point is replaced with the median of all pixels in its neighborhood, and the median of all pixels is calculated using the following formula: F med =medM (x,y) =med[M (x+m,y+n) ;(m,n)∈A] (2) In the formula, F med A is the median of all pixels, and A is the center pixel. x The neighborhood of y), where m and n are the length and width of the neighborhood window, M (x,y) Here, M represents the pixel value at (x, y), med is the median value, and M is the median value. (x+m,y+n) The pixel value at pixel (x+m, y+n); Sort all pixel values ​​within the neighborhood and obtain the minimum value F of all pixels. min Median F med Maximum value F max ; Execute the first step: If d1 = F med -F min >0 and d2=F max -F med If the value is greater than 0, proceed to step three; otherwise, proceed to step two. The second step includes: expanding the neighborhood area by one ring; if the neighborhood area satisfies the second condition, then sorting all pixel values ​​within the neighborhood again and updating the minimum value F of all pixels. min Median F med Maximum value F max Then repeat step one; if the neighborhood area does not meet the second condition, then do not replace the grayscale value of the center pixel; wherein, the second condition is that the neighborhood area is less than the maximum neighborhood threshold S. max ; The third step includes: if g1 = M c -F min >0 and g2=F max -M c If the grayscale value is greater than 0, then the grayscale value of the center pixel is not replaced; otherwise, the grayscale value of the center pixel is replaced with the median value F. med .

4. The method for extracting root canal features based on CBCT data according to claim 2, characterized in that, The optimal threshold for pulp segmentation is calculated based on the intercept histogram of inverse cross-entropy, resulting in preprocessed CBCT sequence images, including: Based on the denoised CBCT images, gamma transform is used to enhance image contrast, and the neighborhood average image and gradient image are calculated respectively. The neighborhood average image and gradient image are combined with the grayscale image and added pixel by pixel to form a hybrid image; The optimal threshold for pulp segmentation is calculated using the intercept histogram of the mixed image; Each pixel value in the blended image is compared with the optimal threshold and converted into a binary image; Remove all unnecessary connected tissues on the image boundaries of the binary image, retaining only the pulp region, to obtain a preprocessed CBCT sequence image; The step of calculating the optimal threshold for pulp segmentation using the intercept histogram of the mixed image includes: A three-dimensional space is generated by grayscale information, neighborhood average grayscale information and gradient compensation. The three-dimensional space is divided into target region O and background region B by a plane. Based on the aforementioned three-dimensional space, the final gradient is calculated using the following formula: In the formula, G H (x,y) and G V (x,y) represents the gradient values ​​of a point (x,y) in the horizontal and vertical directions in a given image, where x represents the x-coordinate of the point in the given image, y represents the y-coordinate of the point in the given image, G(x,y) represents the final gradient, g(x,y+1) represents the gray value of point (x,y+1), g(x,y-1) represents the gray value of point (x,y-1), g(x+1,y) represents the gray value of point (x+1,y), and g(x-1,y) represents the gray value of point (x+1,y). The hybrid image is represented as: F(x,y)=f(x,y)+g(x,y)-G(x,y) (5) In the formula, F(x,y) represents the mixed image, and g(x,y) and f(x,y) represent the gray value of point (x,y) and the average gray value of the neighborhood of point (x,y). Based on the segmentation threshold T, the hybrid image is divided into the target portion Ω. O ∈{(x,y)|F(x,y)=0,1,…,T} and background part Ω B ∈{(x,y)|F(x,y)=T+1,T+2,…,2T-2}, pixel values ​​less than the segmentation threshold T belong to the target part, and pixel values ​​greater than the segmentation threshold T belong to the background part; The inverse cross-entropy between the target region and the background region is calculated using the following formula: In the formula, E(O,B) represents the reciprocal cross-entropy; P(k) represents the prior probability that the average gray value of the neighborhood of the mixed image is k. This represents the average gray level of the target region. This represents the average grayscale value of the background area; Based on the calculated cross-entropy E(O,B), the optimal threshold for pulp segmentation is determined using the following formula: In the formula, The optimal threshold for pulp segmentation is given by arg min, which is the index operation for finding the minimum element in the array.

5. The method for extracting root canal features based on CBCT data according to claim 4, characterized in that, The surface model of the root canal is drawn using surface rendering of the preprocessed CBCT sequence images, including: In three-dimensional space, uniform sampling is performed on the upper and lower jaw regions in the x, y, and z directions with sampling intervals of Δx, Δy, and Δz. The volume data is represented by a ternary function; where a cubic region consisting of eight adjacent sampling points is a voxel. For each voxel, an isosurface is extracted; wherein, the isosurface is a surface composed of points with the same properties in three-dimensional space; Traverse the entire volume data to find the voxels containing isosurfaces, and determine the intersection points P(x,y,z) of the isosurfaces and the edges of the voxels using the following formula: In the formula, T is the isosurface threshold, M1 and M2 are the gray values ​​of the first and second vertices on the edge where the voxel intersects with the isosurface, x1, y1, z1 are the position coordinates of the first vertex on the edge where the voxel intersects with the isosurface, x2, y2, z2 are the position coordinates of the second vertex on the edge where the voxel intersects with the isosurface, and x, y, z are the position coordinates of the intersection points of the isosurface and each edge of the voxel. Using the intersection point P(x,y,z) of the isosurface and each edge of the voxel, and the triangular facets formed by the first vertex P1(x1,y1,z1) and the second vertex P2(x2,y2,z2) on the edge where the voxel intersects the isosurface as isosurfaces, the isosurfaces are combined to form an isosurface triangular network. To render the equivalent triangular network with better 3D effects, a lighting model needs to be selected based on the normal vectors of each triangular facet to render the equivalent triangular network and obtain the surface model of the root canal; wherein, the formula for calculating the normal vectors of each triangular facet is: In the formula, G x G y G z M represents the gradient at the vertex of the isosurface along the X, Y, and Z axes, respectively. (x+a,y,z) M (x-a,y,z) M represents the grayscale value of two points (x+a,y,z) and (xa,y,z) on the edge where the isosurfaces intersect. (x,y+b,z) M (x,y-b,z) M represents the grayscale values ​​of (x,y+b,z) and (x,yb,z). (x,y,z+c) M (x,y,z-c) V represents the grayscale values ​​of (x,y,z+c) and (x,y,zc), where a, b, and c represent the distances between the vertices of the two isosurfaces along the X, Y, and Z axes, respectively; x V y V z The normal vectors G of the triangular facets along the X, Y, and Z axes, respectively. x1 G y1 G z1 Let G be the gradient components of point P1(x1,y1,z1) along the X, Y, and Z axes, respectively. x2 G y2 G z2 These are the gradient components of point P2(x2,y2,z2) along the X, Y, and Z axes, respectively.

6. The method for extracting root canal features based on CBCT data according to claim 1, characterized in that, The root canal physical model is aligned in three-dimensional space, and the morphological characteristics of the root canal are measured by fitting an approximate centerline. This includes: x i =(x i ,y i ,z i (i = 1, 2, ..., n) can be represented as matrix X: Calculated using the following formula Covariance matrix: In the formula, V represents The covariance matrix, S represents the point cloud stretching matrix, satisfying S = S T R represents the point cloud rotation matrix, L represents the covariance in the new coordinate system after the point cloud is dimensionality reduced and transformed, and D represents the covariance matrix of the white data. If the eigenvalues ​​of V are λ1, λ2, ..., λ n The corresponding feature vectors are p1, p2, ..., p n Then we have: V·p i =λ i ·V (17) It can be seen that the eigenvectors correspond to the rotation matrix R, representing the coordinate axis direction of each component; the eigenvalues ​​correspond to the stretching matrix S, representing the variance of the data in the corresponding component's coordinate axis direction. Arrange the eigenvectors into a matrix from top to bottom according to the magnitude of their corresponding eigenvalues, and take the first m principal components to form a matrix W. m And determine the Z-axis vector of the principal components of the data after the original point cloud is reduced to m dimensions as follows: Let the Z-axis vector of the principal components obtained according to equation (18) be k = (z x ,z y ,z z ), with i = (1, -z) x / z y The X-axis vector is used as the X-axis vector, and the Y-axis vector is obtained by cross product of the X-axis and Z-axis vectors, which determines the principal component directions of the three-dimensional point cloud data of teeth. The root canal length, orifice orientation, root canal curvature, and pulp depth are measured using a centerline approximation method in a coordinate system aligned with the principal component directions.

7. A device for extracting root canal features based on CBCT data, used to implement the method as described in any one of claims 1 to 6, characterized in that, The device includes: The 3D reconstruction module is configured to filter and denoise CBCT images and segment the pulp, and to perform surface rendering and 3D reconstruction of the root canals to generate a 3D model. The surface reconstruction module is configured to use the point-surface method to reconstruct the surface based on the regular point cloud data in the three-dimensional model, and generate the root canal solid model. The parameter estimation module is configured to perform coordinate alignment processing on the root canal entity model in three-dimensional space and measure the root canal morphological features by fitting the approximate centerline of the root canal; wherein the root canal morphological features include at least one of root canal length, root canal orifice direction, root canal curvature, and pulp depth.

8. An electronic device, characterized in that, include: Controller; A memory for storing one or more programs, which, when executed by the controller, cause the controller to implement the root canal feature morphology extraction method based on CBCT data as described in any one of claims 1 to 6.

9. A computer-readable storage medium, characterized in that, It stores computer-readable instructions, which, when executed by the computer's processor, cause the computer to perform the root canal feature morphology extraction method based on CBCT data as described in any one of claims 1 to 6.

Images